c Add a call to NINT to resolve a warning issued by the g95 Fortran compiler. c SJP 2009/04/15 c c Minor fixes to resolve warnings issued by the g95 Fortran compiler. c SJP 2009/04/14 c c Modified for v6.2-sjp1.4. c SJP 2004/08/09 c c Inserted code to detect floating point exceptions before they occur. If c division by zero is about to occur (loop 481) or the logarithm of a negative c number is about to be taken (loop 619), the run is aborted and the values of c key variables are written to standard output. c Note that this is only done for these two particular calculations, as these c are the only points at which I have ever encountered FPEs. c SJP 2004/01/05 c c Loops re-rolled so that code passes array-bounds checking. c SJP 2003/03/08 c c Parallel compiler directives replaced with equivalent OpenMP instructions. c SJP 2001/12/12 c c $Log: fst88.f,v $ c Revision 1.21 2000/11/16 01:16:47 rot032 c Fix from Hal for R21 on NEC. c c Revision 1.20 1998/07/17 01:07:24 mrd c Fix bug that affected lowest level heating rates. c c Revision 1.19 1997/12/17 23:22:59 ldr c Changes from MRD for parallel execution on NEC. c c Revision 1.18 1997/10/06 07:57:30 ldr c Final corrections to V5-1. c c Revision 1.17 1997/10/03 06:24:20 ldr c Put continuation symbols in column 6 for f77 compilers. c c Revision 1.16 1997/10/03 05:32:08 ldr c Changes for NEC from MRD and HBG c c Revision 1.15 1995/07/05 06:14:32 ldr c Changes from Dave Micklethwaite for Cray parallel execution tidied by LDR c and merged into V4-7-2l. c c Revision 1.13.1.1 1995/07/03 05:48:41 ldr c Mickles' changes to V4-5-30l for Cray multiprocessing, tidied by LDR. c c Revision 1.14 1994/08/08 17:21:20 ldr c Strip off excessive RCS comments at top of file. c c Revision 1.13 94/05/03 11:57:02 ldr c Simplify conversion of fluxes to heating rates (from MRD). c c Revision 1.12 93/12/17 15:32:37 ldr c Hack V4-4-52l to change all continuation chars to & c c Revision 1.11 93/08/31 14:19:00 ldr c Changes to V4-4 to get model running correctly on the VP. c c Revision 1.10 93/08/19 15:08:06 ldr c Minor cosmetic changes. c c Revision 1.9 93/07/27 14:58:46 ldr c Merge contents of common blocks clrtemp and lwoutclr into lwout. c c Revision 1.8 93/07/22 10:17:27 mrd c Add diagnostics for clear sky net surface radiation and downward radiation c at surface. c c Revision 1.7 93/06/23 14:30:32 mrd c Made padding of vtemp common block resolution independent c c Revision 1.6 93/06/16 14:12:38 ldr c Bugfix from Lawrie Rikus (via MRD). c c ***************************************************************** c subroutine fst88 is the main computation module of the c long-wave radiation code. in it all "emissivity" calculations, c including calls to table lookup subroutines. also,after calling c subroutine "spa88", final combined heating rates and ground c flux are obtained. c ***************************************************************** c inputs: c betinw,betawd,ab15wd bdwide c betad,bo3rnd,ao3rnd bandta c qh2o,p,delp2,delp,t,var1,var2, kdacom c var3,var4,cntval kdacom c temp,press radisw c ind,indx2,kmaxv,source,dsrce tabcom c skc1r,skc3r,kmaxvm,nrep1,nrep2 tabcom c nst1,nst2,nrp1,nrp2 tabcom c co2nbl,co2sp,co21 tfcom c outputs: c heatra,grnflx,topflx lwout c flx1e1 rdflux c c called by : radmn or main pgm c calls : clo88,e1e288,e3v88,spa88,nlte c c passed variables: c in e3v88: c emd = e3 function for h2o lines (0-560,1200-2200 cm-1) c computed in e3v88 c tpl = temperature input for e3 calculation in e3v88 c empl = h2o amount,input for e3 calculation in e3v88 c in e1e288: c e1cts1 = e1 function for the (i+1)th level using the c temperature of the ith data level,computed over c the frequency range 0-560,1200-2200 cm-1. (e1cts1- c e1ctw1) is used in obtaining the flux at the top c in the 0-160,1200-2200 cm-1 range (flx1e1). c e1cts2 = e1 function for the ith level, using the temp. of c the ith data level,computed over the frequency range c 0-560,1200-2200 cm-1. (e1cts2-e1ctw2) is also used c in obtaining the flux at the top in the 0-160,. c 1200-2200 cm-1 range. c e1flx = e1 fctn. for the ith level,using the temperature at c the top of the atmosphere. computed over the freq. c range 0-560,1200-2200 cm-1. used for q(approx) term. c (in common block tfcom) c e1ctw1 = like e1cts1,but computed over the 160-560 cm-1 range c and used for q(approx,cts) calculation c e1ctw2 = like e1cts2,but computed over the 160-560 cm-1 range c and used for q(approx,cts) calculation c fxo = temperature index used for e1 function and also c used for source function calc. in fst88. c dt = temp. diff.between model temps. and temps. at c tabular values of e1 and source fctns. used in c fst88 and in e1 function calc. *VOCL TOTAL,REPEAT(999999) subroutine fst88 !$OMP THREADPRIVATE ( /CLDCOM/ ) !$OMP THREADPRIVATE ( /CLRFLX/ ) !$OMP THREADPRIVATE ( /KDACOM/ ) !$OMP THREADPRIVATE ( /LWOUT/ ) !$OMP THREADPRIVATE ( /RADISW/ ) !$OMP THREADPRIVATE ( /RDFLUX/ ) !$OMP THREADPRIVATE ( /SRCCOM/ ) !$OMP THREADPRIVATE ( /TFCOM/ ) !$OMP THREADPRIVATE ( /VTEMP/ ) include 'HCON.f' include 'PARAMS.f' include 'RDPARM.f' include 'RADISW.f' include 'RNDDTA.f' include 'TABCOM.f' c include 'RDFLUX.f' include 'LWOUT.f' include 'KDACOM.f' include 'SRCCOM.f' include 'TFCOM.f' include 'CLDCOM.f' common / vtemp / ixo(imax,lp1),itop(imax),ibot(imax),indtc(imax) common / vtemp / vtmp1(imax,lp1),vtmp2(imax,lp1), & vtmp3(imax,lp1),c(imax,llp1),alp(imax,llp1),dsorc(imax,lp1), & totphi(imax,lp1),toto3(imax,lp1),tphio3(imax,lp1),rlog(imax,l), & delptc(imax),ptop(imax),pbot(imax),ftop(imax), & fbot(imax) common /vtemp/ over1d !Need over1d in common for SGI common /vtemp/ dummy(imax*(l*l+2*l)) dimension vtmp3x(imax+1,lp1),dsorcx(imax+1,lp1) c---dimension of variables equivalenced to those in vtemp--- dimension tval(imax,lp1),vsum1(imax,lp1),heatem(imax,lp1) dimension emxx(imax,l) dimension csub(imax,llp1) ,csub2(imax,llp1),c2(imax,llp1) dimension flx(imax,lp1),sum(imax,lp1) dimension oss(imax,lp1),css(imax,lp1),ss2(imax,lp1),tc(imax,lp1), & dtc(imax,lp1) dimension alpsq1(imax,lp1),alpsq2(imax,lp1) dimension delpr1(imax,lp1),delpr2(imax,lp1) dimension flxnet(imax,lp1),flxthk(imax,lp1) dimension z1(imax,lp1),ceval(imax,lp1) dimension cevalclr(imax,lp1) dimension totevv(imax,lp1) dimension avmo3(imax,lp1),avpho3(imax,lp1),fac1(imax,lp1) dimension avvo2(imax,lp1),avephj(imax,lp1) dimension over(imax,lp1,lp1),over1d(imax,lp1m) dimension emisst(imax,lp1,lp1) c---dimension of variables equivalenced to those in other common blks-- dimension to31d(imax,lp1m),emi21d(imax,lp1m) dimension co21d(imax,lp1m),emis1d(imax,lp1m) dimension avep1d(imax,lp1m),avep1(imax*lp1m) c---dimension of variables passed to other subroutines--- dimension e1cts1(imax,lp1),e1cts2(imax,l) dimension e1ctw1(imax,lp1),e1ctw2(imax,l) dimension fxo(imax,lp1),dt(imax,lp1) dimension emd(imax,llp1),tpl(imax,llp1),empl(imax,llp1) c---emx1 is a local variable used as input and output to e1e288-- dimension emx1(imax) equivalence (vtmp3,tval,vsum1,emxx,tc,heatem) equivalence (dsorc,totevv,delpr1,oss,sum,flxnet) equivalence (totphi,delpr2,css,flx,z1) equivalence (toto3,alpsq1,dtc,flxthk) equivalence (tphio3,alpsq2,ss2,ceval) equivalence (avephi,avep1d,avep1),(emi21d,emiss2) equivalence (over1d,over),(to31d,to3),(emis1d,emiss) equivalence (co21d,co21) equivalence (alp,csub) real flxclr(imax,lp1),vsum1clr(imax,lp1),vtmp1clr(imax,lp1) real ctsclr(imax,l) real heatemclr(imax,lp1),heatraclr(imax,lp1) integer i, j, k real temp_sjp c c c first section is table lookup for source function and c derivative (b and db/dt).also,the nlte co2 source function c is obtained c c---decrementing the index by 9 accounts for the tables beginning c at t=100k. CSJP do 101 i=1,imax*lp1 CSJP vtmp2(i,1)=aint(temp(i,1)*hp1) CSJP fxo(i,1)=vtmp2(i,1)-9. CSJP dt(i,1)=temp(i,1)-ten*vtmp2(i,1) CSJP ixo(i,1)=fxo(i,1) do 101 j = 1, lp1 do 101 i = 1, imax vtmp2(i, j) = aint(temp(i, j)*hp1) fxo(i, j) = vtmp2(i, j) - 9.0 dt(i, j) = temp(i, j) - ten*vtmp2(i, j) ixo(i, j) = nint(fxo(i, j)) 101 continue c c---source function for combined band 1 do 4114 i=1,imax do 4114 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),1) dsorcx(i,k)=dsrce(ixo(i,k),1) 4114 continue do 4112 k=1,lp1 do 4112 i=1,imax sorc(i,k,1)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 4112 continue c---source function for combined band 2 do 4214 i=1,imax do 4214 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),2) dsorcx(i,k)=dsrce(ixo(i,k),2) 4214 continue do 4212 k=1,lp1 do 4212 i=1,imax sorc(i,k,2)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 4212 continue c---source function for combined band 3 do 4314 i=1,imax do 4314 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),3) dsorcx(i,k)=dsrce(ixo(i,k),3) 4314 continue do 4312 k=1,lp1 do 4312 i=1,imax sorc(i,k,3)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 4312 continue c---source function for combined band 4 do 4414 i=1,imax do 4414 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),4) dsorcx(i,k)=dsrce(ixo(i,k),4) 4414 continue do 4412 k=1,lp1 do 4412 i=1,imax sorc(i,k,4)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 4412 continue c---source function for combined band 5 do 4514 i=1,imax do 4514 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),5) dsorcx(i,k)=dsrce(ixo(i,k),5) 4514 continue do 4512 k=1,lp1 do 4512 i=1,imax sorc(i,k,5)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 4512 continue c---source function for combined band 6 do 4614 i=1,imax do 4614 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),6) dsorcx(i,k)=dsrce(ixo(i,k),6) 4614 continue do 4612 k=1,lp1 do 4612 i=1,imax sorc(i,k,6)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 4612 continue c---source function for combined band 7 do 4714 i=1,imax do 4714 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),7) dsorcx(i,k)=dsrce(ixo(i,k),7) 4714 continue do 4712 k=1,lp1 do 4712 i=1,imax sorc(i,k,7)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 4712 continue c---source function for combined band 8 do 4814 i=1,imax do 4814 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),8) dsorcx(i,k)=dsrce(ixo(i,k),8) 4814 continue do 4812 k=1,lp1 do 4812 i=1,imax sorc(i,k,8)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 4812 continue c---source function for band 9 (560-670 cm-1) do 4914 i=1,imax do 4914 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),9) dsorcx(i,k)=dsrce(ixo(i,k),9) 4914 continue do 4912 k=1,lp1 do 4912 i=1,imax sorc(i,k,9)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 4912 continue c---source function for band 10 (670-800 cm-1) do 5014 i=1,imax do 5014 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),10) dsorcx(i,k)=dsrce(ixo(i,k),10) 5014 continue do 5012 k=1,lp1 do 5012 i=1,imax sorc(i,k,10)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 5012 continue c---source function for band 11 (800-900 cm-1) do 5114 i=1,imax do 5114 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),11) dsorcx(i,k)=dsrce(ixo(i,k),11) 5114 continue do 5112 k=1,lp1 do 5112 i=1,imax sorc(i,k,11)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 5112 continue c---source function for band 12 (900-990 cm-1) do 5214 i=1,imax do 5214 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),12) dsorcx(i,k)=dsrce(ixo(i,k),12) 5214 continue do 5212 k=1,lp1 do 5212 i=1,imax sorc(i,k,12)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 5212 continue c---source function for band 13 (990-1070 cm-1) do 5314 i=1,imax do 5314 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),13) dsorcx(i,k)=dsrce(ixo(i,k),13) 5314 continue do 5312 k=1,lp1 do 5312 i=1,imax sorc(i,k,13)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 5312 continue c---source function for band 14 (1070-1200 cm-1) do 5414 i=1,imax do 5414 k=1,lp1 vtmp3x(i,k)=source(ixo(i,k),14) dsorcx(i,k)=dsrce(ixo(i,k),14) 5414 continue do 5412 k=1,lp1 do 5412 i=1,imax sorc(i,k,14)=vtmp3x(i,k)+dt(i,k)*dsorcx(i,k) 5412 continue c c c the following subroutine obtains nlte source function for co2 c c c---obtain special source functions for the 15 um band (csour),the c 9.6 um band (osour) and the window region (ss1) CSJP do 131 i=1,imax*lp1 CSJP ss1(i,1)=sorc(i,1,11)+sorc(i,1,12)+sorc(i,1,14) do 131 j = 1, lp1 do 131 i = 1, imax ss1(i, j) = sorc(i, j, 11) + sorc(i, j, 12) + sorc(i, j, 14) 131 continue CSJP do 143 i=1,imax*lp1 CSJP csour(i,1)=sorc(i,1,9)+sorc(i,1,10) CSJP osour(i,1)=sorc(i,1,13) do 143 j = 1, lp1 do 143 i = 1, imax csour(i, j) = sorc(i, j, 9) + sorc(i, j, 10) osour(i, j) = sorc(i, j, 13) 143 continue c c c second section produces 4 matrices for later use in program: c 1) avephi(i,j) is the scaled water mass for use in tables (this c is u(p,p') in ref.(4); c 2) over (i,j) is the water transmission function (using c "emissivity" approximation) in the 560-800 cm-1 band; c 3) to3(i,j) is the exact ozone transmission function (using c parameters for the 990-1070 cm-1 band from the 1982 afgl catalog) c 4)emiss2(i,j) is the emissivity h20 transmission function due c to the 10 um continuum,treated as one band. c do 201 i=1,imax totphi(i,1)=zero toto3(i,1)=zero tphio3(i,1)=zero totvo2(i,1)=zero 201 continue do 203 k=2,lp1 do 203 i=1,imax totphi(i,k)=totphi(i,k-1)+var2(i,k-1) toto3(i,k)=toto3(i,k-1)+var3(i,k-1) tphio3(i,k)=tphio3(i,k-1)+var4(i,k-1) totvo2(i,k)=totvo2(i,k-1)+cntval(i,k-1) 203 continue c c the calculational procedure used here is: c 1) from the 1-d vectors (totphi,toto3,etc) construct a 1-d c arrays including all information for the upper (i>j) triangle c of the symmetric matrices(lp1,lp1); c 2)perform computations on these arrays to obtain to31d,over1d, c avephj: upper triangle of the symmetric matrices in 1-d form c 3) fill up the lower triangle c the diagram below illustrates the relationship between the 1-d c array and the 2-d symmetric matrix for a 4x4 matrix. c c i c 1 2 3 4 c -------------------------- c 1 1 2 3 the nos. are the c j 2 4 5 positions in the c 3 6 1-d array c 4 c c c compute "upper triangle" transmission functions for c the 9.6 um band (vtmp1) and the 15 um band (vtmp3). also, c the c stage 1....compute o3 ,over transmission fctns and avephi c---do k=1 calculation (from flux layer kk to the top) separately c as vectorization is improved,and ozone cts transmissivity c may be extracted here. CSJP do 302 i=1,imax*l CSJP avmo3(i,1)=toto3(i,2) CSJP avpho3(i,1)=tphio3(i,2) CSJP avephj(i,1)=totphi(i,2) CSJP avvo2(i,1)=totvo2(i,2) CSJP fac1(i,1)=bo3rnd(2)*avpho3(i,1)/avmo3(i,1) CSJP vtmp2(i,1)=haf*(fac1(i,1)* CSJP & (sqrt(one+(four*ao3rnd(2)*avmo3(i,1))/fac1(i,1))-one)) CSJP vtmp1(i,1)=exp(hm1ez*(vtmp2(i,1)+sko3r*avvo2(i,1))) CSJP vtmp3(i,1)=exp(hm1ez*(sqrt(ab15wd*avephj(i,1))+ CSJP & skc1r*avvo2(i,1))) do 302 j = 1, l do 302 i = 1, imax avmo3(i, j) = toto3(i, j+1) avpho3(i, j) = tphio3(i, j+1) avephj(i, j) = totphi(i, j+1) avvo2(i, j) = totvo2(i, j+1) fac1(i, j) = bo3rnd(2) * avpho3(i, j) / avmo3(i, j) vtmp2(i, j) = haf * (fac1(i, j) * & (sqrt(one+(four*ao3rnd(2)*avmo3(i, j))/fac1(i, j))-one)) vtmp1(i, j) = exp(hm1ez*(vtmp2(i, j)+sko3r*avvo2(i, j))) vtmp3(i, j) = exp(hm1ez*(sqrt(ab15wd*avephj(i, j)) + & skc1r*avvo2(i, j))) c---the next 4 lines fill up the upper triangle of the to3,avephi c and over matrices, and all of the to3spc array. the notation c is obscure; note that an equivalent statement for line 1 would be: c for i=1 to imax and kp=k+1(=2) to lp1: to3(i,kp,1)=vtmp1(i,kp) c c I *think* this above comment is wrong, and that the correct c equivalent statement would be: to3(i, kp, 1) = vtmp1(i, kp-1) c SJP 2003/03/08 CSJP to3(i,2,1)=vtmp1(i,1) CSJP avephi(i,2,1)=avephj(i,1) CSJP over(i,2,1)=vtmp3(i,1) CSJP to3spc(i,1)=vtmp2(i,1) to3(i, j+1, 1) = vtmp1(i, j) avephi(i, j+1, 1) = avephj(i, j) over(i, j+1, 1) = vtmp3(i, j) to3spc(i, j) = vtmp2(i, j) 302 continue c---fill in lower triangle of to3,over arrays. it is unnecessary to c fill in avephi array at this time. do 307 kp=2,lp1 do 307 i=1,imax to3(i,1,kp)=vtmp1(i,kp-1) over(i,1,kp)=vtmp3(i,kp-1) 307 continue c---now repeat for the k=2..l cases. do 321 k=2,l do 322 kk=1,lp1-k do 322 i=1,imax avmo3(i,kk)=toto3(i,kk+k)-toto3(i,k) avpho3(i,kk)=tphio3(i,kk+k)-tphio3(i,k) avephj(i,kk)=totphi(i,kk+k)-totphi(i,k) avvo2(i,kk)=totvo2(i,kk+k)-totvo2(i,k) 322 continue CSJP Former machine dependence at this point do kk=1,lp1-k do i=1,imax fac1(i,kk)=bo3rnd(2)*avpho3(i,kk)/avmo3(i,kk) vtmp2(i,kk)=haf*(fac1(i,kk)* & (sqrt(one+(four*ao3rnd(2)*avmo3(i,kk))/fac1(i,kk))-one)) vtmp1(i,kk)=exp(hm1ez*(vtmp2(i,kk)+sko3r*avvo2(i,kk))) vtmp3(i,kk)=exp(hm1ez*(sqrt(ab15wd*avephj(i,kk))+ & skc1r*avvo2(i,kk))) to3(i,k+kk,k)=vtmp1(i,kk) avephi(i,k+kk,k)=avephj(i,kk) over(i,k+kk,k)=vtmp3(i,kk) end do end do do 327 kp=k+1,lp1 do 327 i=1,imax to3(i,k,kp)=vtmp1(i,kp-k) over(i,k,kp)=vtmp3(i,kp-k) 327 continue 321 continue c---initialize diagonal elements of to3,over (not needed for avephi) do 309 k=1,lp1 do 309 i=1,imax to3(i,k,k)=1. over(i,k,k)=1. 309 continue c*****stage 4....compute one-band continuum transmissivities (emiss2) c CSJP do 481 i=1,imax*lp1 CSJP vtmp3(i,1)=exp(hm1ez*totvo2(i,1)) CSJP totevv(i,1)=one/vtmp3(i,1) do 481 j = 1, lp1 do 481 i = 1, imax vtmp3(i, j) = exp(hm1ez*totvo2(i, j)) if (vtmp3(i, j) .eq. 0.0) then write (*, *) write (*, *) "ABORTING: Fatal error in FST88" write (*, *) write (*, *) "i = ", i write (*, *) "j = ", j write (*, *) write (*, *) "totvo2(i, j) = ", totvo2(i, j) write (*, *) stop end if totevv(i, j) = one / vtmp3(i, j) 481 continue do 501 k=1,lp1 do 501 kp=1,k do 503 i=1,imax emiss2(i,kp,k)=vtmp3(i,k)*totevv(i,kp) 503 continue 501 continue do 505 k=1,l do 505 kp=k+1,lp1 do 507 i=1,imax emiss2(i,kp,k)=vtmp3(i,kp)*totevv(i,k) 507 continue 505 continue c the third section calculates boundary layer and nearby layer c corrections to the transmission functions obtained above. methods c are given in ref. (4). c combine co21,over into co21; before making nbl corrections, c load (1,k) values for use in exact cts calculations in spa88. CSJP do 605 i=1,imax*lp1*lp1 CSJP co21(i,1,1)=co21(i,1,1)*over(i,1,1) do 605 k = 1, lp1 do 605 j = 1, lp1 do 605 i = 1, imax co21(i, j, k) = co21(i, j, k) * over(i, j, k) 605 continue do 607 k=1,lp1 do 607 i=1,imax co2sp(i,k)=co21(i,1,k) 607 continue c the following ratios are used in various nbl calculations: c do 619 k=1,l do 619 i=1,imax temp_sjp = over(i, k, k+1) * co2nbl(i, k) if (temp_sjp .lt. 0.0) then write (*, *) write (*, *) "ABORTING: Fatal error in FST88" write (*, *) write (*, *) "i = ", i write (*, *) "k = ", k write (*, *) write (*, *) "over(i, k, k+1) = ", over(i, k, k+1) write (*, *) "co2nbl(i, k) = ", co2nbl(i, k) write (*, *) stop end if rlog(i,k)=log(temp_sjp) 619 continue CSJP do 601 i=1,imax*lm1 CSJP delpr1(i,2)=delp(i,2)*(press(i,2)-p(i,2)) CSJP alpsq1(i,2)=sqrt(delpr1(i,2)) CSJP alp(i,lp1)=-alpsq1(i,2)*rlog(i,2) do 601 j = 1, lm1 do 601 i = 1, imax delpr1(i, j+1) = delp(i, j+1) * (press(i, j+1) - p(i, j+1)) alpsq1(i, j+1) = sqrt(delpr1(i, j+1)) alp(i, j+lp1-1) = -alpsq1(i,j+1) * rlog(i, j+1) 601 continue CSJP do 603 i=1,imax*l CSJP delpr2(i,2)=delp(i,1)*(p(i,2)-press(i,1)) CSJP alpsq2(i,2)=sqrt(delpr2(i,2)) CSJP alp(i,1)=-alpsq2(i,2)*rlog(i,1) do 603 j = 1, l do 603 i = 1, imax delpr2(i, j+1) = delp(i, j) * (p(i, j+1) - press(i, j)) alpsq2(i, j+1) = sqrt(delpr2(i, j+1)) alp(i, j) = -alpsq2(i, j+1) * rlog(i, j) 603 continue do 625 i=1,imax alp(i,ll)=-rlog(i,l) alp(i,llp1)=-rlog(i,l)*sqrt(delp(i,l)*(p(i,lp1)-press(i,lm1))) 625 continue c the first computation is for the 15 um band,with the c for the combined h2o and co2 transmission function. c c perform nbl computations for the 15 um band c***the statement function sf in prev. versions is now explicitly c evaluated. CSJP do 631 i=1,imax*llp1 c!!!!!!!!!! c c(i,1)=alp(i,1)*(hmp66667+alp(i,1)*(quartr+alp(i,1)*hm6666m2)) CSJP c(i,1)= 2.0/alp(i,1)**2 * CSJP & (1.0-exp(-alp(i,1)) * (alp(i,1)+1.0)) - 1.0 do 631 j = 1, llp1 do 631 i = 1, imax c(i, j) = 2.0 / alp(i, j)**2 * & (1.0-exp(-alp(i, j)) * (alp(i, j)+1.0)) - 1.0 631 continue do 641 i=1,imax co21(i,lp1,lp1)=one+c(i,l) co21(i,lp1,l)=one+(delp2(i,l)*c(i,ll)-(press(i,l)-p(i,l))* & c(i,llm1))/(p(i,lp1)-press(i,l)) co21(i,l,lp1)=one+((p(i,lp1)-press(i,lm1))*c(i,llp1)- & (p(i,lp1)-press(i,l))*c(i,l))/(press(i,l)-press(i,lm1)) 641 continue c*** the k indices in the following loop run from 21 to 40 in the c l40 skyhi code version do 643 k=2,l do 643 i=1,imax co21(i,k,k)=one+haf*(c(i,lm1+k)+c(i,k-1)) 643 continue c c compute nearby-layer transmissivities for the o3 band and for the c one-band continuum band (to3 and emiss2). the sf2 function is c used. the method is the same as described for co2 in ref (4). CSJP do 651 i=1,imax*lm1 CSJP csub(i,2)=cntval(i,2)*delpr1(i,2) do 651 j = 1, lm1 do 651 i = 1, imax csub(i, j+1) = cntval(i, j+1) * delpr1(i, j+1) 651 continue CSJP do i=1,imax*l CSJP csub(i,lp1)=cntval(i,1)*delpr2(i,2) do j = 1, l do i = 1, imax csub(i, j+lp1-1) = cntval(i, j) * delpr2(i, j+1) end do end do c---the sf2 function in prev. versions is now explicitly evaluated c!!!!!!!!!!! c do 655 i=1,imax*llm2 CSJP do 655 i=1,imax*llm1 CSJP csub2(i,2)=sko3r*csub(i,2) CSJP c(i,2)=csub(i,2)*(hmp5+csub(i,2)*(hp166666-csub(i,2)*h41666m2)) CSJP c2(i,2)=csub2(i,2)*(hmp5+csub2(i,2)* CSJP & (hp166666-csub2(i,2)*h41666m2)) do 655 j = 1, llm1 do 655 i = 1, imax csub2(i, j+1) = sko3r * csub(i, j+1) c(i, j+1) = csub(i, j+1) * (hmp5 + csub(i, j+1) * (hp166666 - & csub(i, j+1) * h41666m2)) c2(i, j+1) = csub2(i, j+1) * (hmp5 + csub2(i, j+1) * & (hp166666 - csub2(i, j+1) * h41666m2)) 655 continue do 661 i=1,imax c!!!!!!!!!!! c emiss2(i,lp1,lp1)=one+c(i,llm1) c to3(i,lp1,lp1)=one+c2(i,llm1) emiss2(i,lp1,lp1)=one+c(i,ll) to3(i,lp1,lp1)=one+c2(i,ll) 661 continue do 663 k=2,l do 663 i=1,imax emiss2(i,k,k)=one+haf*(c(i,k)+c(i,lm1+k)) to3(i,k,k)=one+haf*(c2(i,k)+c2(i,lm1+k)) 663 continue c c fourth section obtains water transmission functions c used in q(approx) calculations and also makes nbl corrections: c 1) emiss (i,j) is the transmission function matrix obtained c by calling subroutine e1e288; c 2) "nearby layer" corrections (emiss(i,i)) are obtained c using subroutine e3v88; c 3) special values at the surface (emiss(l,lp1),emiss(lp1,l), c emiss(lp1,lp1)) are calculated. c c emxx,av1,and empl are computed before avephi is modified c c obtain arguments for e1e288 and e3v88: c do 801 i=1,imax emx1(i)=qh2o(i,l)*press(i,l)*(press(i,l)-p(i,l))*gp0inv 801 continue do 803 k=1,lm1 do 803 i=1,imax emxx(i,k)=avephi(i,l,k)+emx1(i) 803 continue CSJP do 811 i=1,imax*l CSJP empl(i,2)=qh2o(i,1)*p(i,2)*(p(i,2)-press(i,1))*gp0inv do 811 j = 1, l do 811 i = 1, imax empl(i, j+1) = qh2o(i, j) * p(i, j+1) * (p(i, j+1) - & press(i, j)) * gp0inv 811 continue CSJP do 812 i=1,imax*lm1 CSJP empl(i,lp2)=qh2o(i,2)*p(i,2)*(press(i,2)-p(i,2))*gp0inv do 812 j = 1, lm1 do 812 i = 1, imax empl(i, j+lp2-1) = qh2o(i, j+1) * p(i, j+1) * & (press(i, j+1) - p(i, j+1)) * gp0inv 812 continue do 821 i=1,imax empl(i,1)=avephi(i,lp1,l) empl(i,llp1)=empl(i,ll) tpl(i,1)=temp(i,l) tpl(i,lp1)=haf*(t(i,lp1)+temp(i,l)) tpl(i,llp1)=haf*(t(i,l)+temp(i,l)) 821 continue do 823 k=2,l do 823 i=1,imax tpl(i,k)=t(i,k) tpl(i,k+l)=t(i,k) 823 continue do 827 k=2,l do 829 i=1,imax avephi(i,k,k)=emxx(i,k-1) 829 continue 827 continue c do 833 i=1,imax avephi(i,lp1,lp1)=avephi(i,lp1,l)+empl(i,l) c---a suitable value is assigned to avephi(i,1,1);it is unused by the c code and needed only to avoid indefinite values avephi(i,1,1)=h101m16 833 continue c compute logs of water mass arguments for e1e288; c the corresponding quantity for e3v88 (empl) must be obtained c within that subroutine, as empl is used after e3v88 is called. c c---only the upper triangle values of avephi exist ;only these values c have their logarithms taken. c from Hal 990916 do 841 k=1,lp1 if(lw.eq.22)then do j=k,lp1 do i=1,imax avephi(i,j,k)=log10(avephi(i,j,k)) end do end do else CSJP do i=1,imax*(lp2-k) CSJP avephi(i,k,k)=log10(avephi(i,k,k)) CSJP enddo do j = 1, lp2-k do i = 1, imax avephi(i, j+k-1, k) = log10(avephi(i, j+k-1, k)) end do end do endif c---now fill in the lower triangle elements of avephi do 843 kp=k+1,lp1 do 843 i=1,imax avephi(i,k,kp)=avephi(i,kp,k) 843 continue 841 continue c c call e1e288 for emissivity transmission fctns for h2o call e1e288(e1cts1,e1cts2,e1flx,e1ctw1,e1ctw2,fxo,dt) c do 853 k=1,lm1 do 855 i=1,imax emiss(i,lp1,k)=haf*(emiss(i,k+1,k+1)+emiss(i,lp1,k)) 855 continue 853 continue c c call e3v88 for nbl h2o transmissivities call e3v88(emd,tpl,empl) c c compute nearby layer and special-case transmissivities for emiss c using methods for h2o given in ref. (4) do 851 k=2,l do 851 i=1,imax emiss(i,k,k)=emd(i,k+l)+emd(i,k) 851 continue c do 861 i=1,imax emiss(i,l,lp1)=(emd(i,1)*empl(i,1)-emd(i,lp1)*empl(i,lp1))/ & emx1(i) + quartr*(emiss(i,l,lp1)+emiss(i,lp1,lp1)) emiss(i,lp1,lp1)=two*emd(i,lp1) emiss(i,lp1,l)=two*(emd(i,1)*empl(i,1)-emd(i,llp1)*empl(i,llp1))/ & (qh2o(i,l)*press(i,l)*(p(i,lp1)-press(i,l))*gp0inv) 861 continue c c subroutine spa88 is called to obtain exact cts for water c co2 and o3, and approximate cts co2 and o3 calculations. c call spa88 c c this section performs the calculation of "emissivity" c fluxes for the 4 components comprising the lw frequency region. c (emiss = the 0-160,1200-2200 cm-1 band; emiss2 the 800-990, c 1070-1200 cm-1 band; to3 the 990-1070 cm-1 band; co21 the 560-800 c cm-1 band). emisst is the combined exchange term and flx the c combined net flux. c CSJP do 901 i=1,imax*lp1 CSJP tc(i,1)=(temp(i,1)*temp(i,1))**2 do 901 j = 1, lp1 do 901 i = 1, imax tc(i, j) = (temp(i, j) * temp(i, j))**2 901 continue CSJP do 903 i=1,imax*l CSJP oss(i,2)=osour(i,2)-osour(i,1) CSJP css(i,2)=csour(i,2)-csour(i,1) CSJP dtc(i,2)=tc(i,2)-tc(i,1) CSJP ss2(i,2)=ss1(i,2)-ss1(i,1) do 903 j = 1, l do 903 i = 1, imax oss(i, j+1) = osour(i, j+1) - osour(i, j) css(i, j+1) = csour(i, j+1) - csour(i, j) dtc(i, j+1) = tc(i, j+1) - tc(i, j) ss2(i, j+1) = ss1(i, j+1) - ss1(i, j) 903 continue c do 905 k=1,lp1 do 905 i=1,imax emisst(i,1,k)=tc(i,1)*e1flx(i,k)+ss1(i,1)* & emiss2(i,k,1)+sorc(i,1,13)*to3(i,k,1)+csour(i,1)*co2sp(i,k) 905 continue do 909 k=1,lp1 CSJP do 909 i=1,imax*l CSJP emisst(i,2,k)=dtc(i,2)*emiss(i,2,k)+ CSJP & ss2(i,2)*emiss2(i,2,k)+ CSJP & oss(i,2)*to3(i,2,k) +css(i,2)*co21(i,2,k) do 909 j = 1, l do 909 i = 1, imax emisst(i, j+1, k) = dtc(i, j+1) * emiss(i, j+1, k) + & ss2(i, j+1) * emiss2(i, j+1, k) + oss(i,j+1) * & to3(i, j+1, k) + css(i, j+1) * co21(i, j+1, k) 909 continue c do 912 k=1,lp1 do 912 i=1,imax flx(i,k)=emisst(i,1,k)*cldfac(i,1,k) flxclr(i,k)=emisst(i,1,k) 912 continue do 911 kp=2,lp1 do 911 k=1,lp1 do 911 i=1,imax flx(i,k)=flx(i,k)+emisst(i,kp,k)*cldfac(i,kp,k) flxclr(i,k)=flxclr(i,k)+emisst(i,kp,k) 911 continue c this section computes the emissivity cts heating rates for 2 c emissivity bands: the 0-160,1200-2200 cm-1 band and the 800- c 990,1070-1200 cm-1 band. the remaining cts comtributions are c contained in ctso3, computed in spa88. c CSJP do 999 i=1,imax*lp1 CSJP vtmp1(i,1)=emiss2(i,1,1)*cldfac(i,1,1) CSJP vtmp1clr(i,1)=emiss2(i,1,1) do 999 j = 1, lp1 do 999 i = 1, imax vtmp1(i, j) = emiss2(i, j, 1) * cldfac(i, j, 1) vtmp1clr(i, j) = emiss2(i, j, 1) 999 continue CSJP do 1001 i=1,imax*l CSJP cts(i,1)=radcon*delp(i,1)*(tc(i,1)* CSJP & (e1ctw2(i,1)*cldfac(i,2,1)-e1ctw1(i,1)*cldfac(i,1,1)) + CSJP & ss1(i,1)*(vtmp1(i,2)-vtmp1(i,1))) CSJP ctsclr(i,1)=radcon*delp(i,1)*(tc(i,1)* CSJP & (e1ctw2(i,1)-e1ctw1(i,1)) + CSJP & ss1(i,1)*(vtmp1clr(i,2)-vtmp1clr(i,1))) do 1001 j = 1, l do 1001 i = 1, imax cts(i, j) = radcon * delp(i, j) * (tc(i, j) * & (e1ctw2(i, j) * cldfac(i, j+1, 1) - e1ctw1(i, j) * & cldfac(i, j, 1)) + ss1(i, j) * (vtmp1(i, j+1) - & vtmp1(i, j))) ctsclr(i, j) = radcon * delp(i, j) * (tc(i, j) * & (e1ctw2(i, j) - e1ctw1(i, j)) + & ss1(i, j) * (vtmp1clr(i, j+1) - vtmp1clr(i, j))) 1001 continue c CSJP do 1011 i=1,imax*l CSJP ceval(i,1)=tc(i,1)*(cldfac(i,1,1)*(e1cts1(i,1)-e1ctw1(i,1)) - CSJP & cldfac(i,2,1)*(e1cts2(i,1)-e1ctw2(i,1))) CSJP cevalclr(i,1)=tc(i,1)*((e1cts1(i,1)-e1ctw1(i,1)) - CSJP & (e1cts2(i,1)-e1ctw2(i,1))) do 1011 j = 1, l do 1011 i = 1, imax ceval(i, j) = tc(i, j) * (cldfac(i, j, 1) * (e1cts1(i, j) - & e1ctw1(i, j)) - cldfac(i, j+1, 1) * & (e1cts2(i, j) - e1ctw2(i, j))) cevalclr(i, j) = tc(i, j) * ((e1cts1(i, j) - e1ctw1(i, j)) - & (e1cts2(i, j) - e1ctw2(i, j))) 1011 continue do 1012 i=1,imax flx1e1(i)=tc(i,lp1)*cldfac(i,lp1,1)* & (e1cts1(i,lp1)-e1ctw1(i,lp1)) flx1e1clr(i)=tc(i,lp1)*(e1cts1(i,lp1)-e1ctw1(i,lp1)) 1012 continue do 1014 k=1,l do 1013 i=1,imax flx1e1(i)=flx1e1(i)+ceval(i,k) flx1e1clr(i)=flx1e1clr(i)+cevalclr(i,k) 1013 continue 1014 continue c c final section obtains emissivity heating rates, c total heating rates and the flux at the ground c c ....calculate the emissivity heating rates CSJP do 1101 i=1,imax*l CSJP heatem(i,1)=radcon*(flx(i,2)-flx(i,1))*delp(i,1) CSJP heatemclr(i,1)=radcon*(flxclr(i,2)-flxclr(i,1))*delp(i,1) do 1101 j = 1, l do 1101 i = 1, imax heatem(i, j) = radcon * (flx(i, j+1) - flx(i, j)) * & delp(i, j) heatemclr(i, j) = radcon * (flxclr(i, j+1) - & flxclr(i, j)) * delp(i, j) 1101 continue c ....calculate the total heating rates CSJP do 1103 i=1,imax*l CSJP heatra(i,1)=heatem(i,1)-cts(i,1)-ctso3(i,1)+excts(i,1) CSJP heatraclr(i,1)=heatemclr(i,1)-ctsclr(i,1)-ctso3clr(i,1) CSJP & +exctsclr(i,1) do 1103 j = 1, l do 1103 i = 1, imax heatra(i, j) = heatem(i, j) - cts(i, j) - ctso3(i, j) + & excts(i, j) heatraclr(i, j) = heatemclr(i, j) - ctsclr(i, j) - & ctso3clr(i, j) + exctsclr(i, j) 1103 continue c ....calculate the flux at each flux level using the flux at the c top (flx1e1+gxcts) and the integral of the heating rates (vsum1) CSJP do 1111 i=1,imax*l CSJP vsum1(i,1)=heatra(i,1)*delp2(i,1)*radcon1 CSJP vsum1clr(i,1)=heatraclr(i,1)*delp2(i,1)*radcon1 do 1111 j = 1, l do 1111 i = 1, imax vsum1(i, j) = heatra(i, j) * delp2(i, j) * radcon1 vsum1clr(i, j) = heatraclr(i, j) * delp2(i, j) * radcon1 1111 continue do 1115 i=1,imax topflx(i)=flx1e1(i)+gxcts(i) flxnet(i,1)=topflx(i) grnflxclr(i)=flx1e1clr(i)+gxctsclr(i) 1115 continue c---only the surface value of flux (grnflx) is needed unless c the thick cloud section is invoked. do 1123 k=2,lp1 do 1123 i=1,imax flxnet(i,k)=flxnet(i,k-1)+vsum1(i,k-1) grnflxclr(i)=grnflxclr(i)+vsum1clr(i,k-1) 1123 continue do 1125 i=1,imax grnflx(i)=flxnet(i,lp1) 1125 continue c c this is the thick cloud section.optionally,if thick cloud c fluxes are to be "convectively adjusted",ie,df/dp is constant, c for cloudy part of grid point, the following code is executed. c***first,count the number of clouds along the lat. row. skip the c entire thick cloud computation if there are no clouds. icnt=0 do 1301 i=1,imax icnt=icnt+nclds(i) 1301 continue if (icnt.eq.0) go to 6999 c---find the maximum number of clouds in the latitude row kclds=nclds(1) do 2106 i=2,imax kclds=max(nclds(i),kclds) 2106 continue c c c***obtain the pressures and fluxes of the top and bottom of c the nc'th cloud (it is assumed that all ktop and kbtm's have c been defined!). do 1361 kk=1,kclds do 1363 i=1,imax j1=ktop(i,kk+1) c if (j1.eq.1) go to 1363 j3=kbtm(i,kk+1)+1 if ((j3-1).gt.j1) then ptop(i)=p(i,j1) pbot(i)=p(i,j3) ftop(i)=flxnet(i,j1) fbot(i)=flxnet(i,j3) c***obtain the "flux derivative" df/dp (delptc) delptc(i)=(ftop(i)-fbot(i))/(ptop(i)-pbot(i)) c***calculate the tot. flux chg. from the top of the cloud, for c all levels. do 1365 k=j1+1,j3-1 z1(i,k)=(p(i,k)-ptop(i))*delptc(i)+ftop(i) flxnet(i,k)=flxnet(i,k)*(one-camt(i,kk+1)*emcld(i,kk+1)) + & z1(i,k)*camt(i,kk+1)*emcld(i,kk+1) 1365 continue endif 1363 continue 1361 continue c***using this flux chg. in the cloudy part of the grid box, obtain c the new fluxes, weighting the clear and cloudy fluxes:again, only c the fluxes in thick-cloud levels will eventually be used. c do 6051 k=1,lp1 c do 6051 i=1,imax c flxnet(i,k)=flxnet(i,k)*(one-camt(i,nc)*emcld(i,nc)) + c 1 z1(i,k)*camt(i,nc)*emcld(i,nc) c051 continue c***merge flxthk into flxnet for appropriate levels. c do 1401 k=1,lp1 c do 1401 i=1,imax c if (k.gt.itop(i) .and. k.le.ibot(i) c 1 .and. (nc-1).le.nclds(i)) then c flxnet(i,k)=flxthk(i,k) c endif c401 continue c c******end of cloud loop***** 6999 continue c***the final step is to recompute the heating rates based on the c revised fluxes: CSJP do 6101 i=1,imax*l c heatra(i,1)=radcon*(flxnet(i,2)-flxnet(i,1))*delp(i,1) CSJP heatra(i,1)=(flxnet(i,2)-flxnet(i,1)) do 6101 j = 1, l do 6101 i = 1, imax heatra(i, j) = (flxnet(i, j+1) - flxnet(i, j)) 6101 continue c the thick cloud section ends here. return end