program RADIA c c This program estimates the longwave and shotwave c radiative fluxes using the emissivity method . c CAUTION : On some machines this program needs c to be compiled with an option that initialize all c unitialized variables to zero.For most compilers c this procedure is the default . With Unix operating c system that option is : c f77 -K filename.f c parameter ( nk =1 2) parameter ( hz = 12.,day = 20,month = 12) parameter (ilevv = 12,ilev2v = ilevv+2) real irtop dimension t(nk),q(nk),sh(nk),s(nk),declsc(2) dimension rtrter(nk),rtrsol(nk) common/wetc/eps,ae(2),be(2),ht(2),tc common/rrr/ fup(ilev2v,8),fdw(ilev2v,8),rw(ilev2v) &,fx(ilev2v),tl(ilev2v),aso(ilev2v,8),rso(ilev2v,8) &,albc(8),rsa,secz,coe(8),dtc(ilev2v,8),fnet(ilev2v,8) &,cld(ilev2v,8),rhx(ilevv),dp(ilevv),qz(ilevv),p(ilev2v) &,tz(ilevv),dpz(ilev2v),qs(ilevv),sso(8),dtw(ilev2v,8) c c Define the logical switches for radiation calls. c logical swc,swr data swr,swc /.true.,.true./ c c Define the grid distance in the x and y directions c and some other required constants. c data dphi,dlambda,stebol /5.5376,5.625,5.6678e-8 / c c Calculation are done on sigma levels . c sigma = p/ps . c data s /.1,.2,.3,.4,.5,.6,.7,.8,.85,.9,.95,.99/ c pi = 3.14159 picon = pi/180. c c Open output file c open (20,file='radia.out',status='unknown') c c Define temperature ,specific humidity for the c vertical column. c data t/220.16,218.06,223.71,235.78,247.31,256.73, & 264.38,269.16,271.98,274.71,277.10,278.41/ c data q/ 4.0412000e-05,1.7426601e-05,2.6705300e-05 & ,9.5105803e-05,2.6963701e-04,6.8341498e-04 & ,1.6110500e-03,2.1570900e-03,2.3708199e-03 & ,2.5780201e-03,2.9228700e-03,3.2207901e-03/ c c Define surface parameters. These variables are explained c in SIMFLX . c data ts,ps /278.96,88230.5/ data rlong,sinlat,coslat,albedo/0.00,0.587,0.810,8.05e-02/ data vent/178.22/ c c Compute declination angle from day and month c call CONRAD (declsc,day,month) c c Compute intermediate sigma levels (sh) c do 9100 k = 1, nk-1 sh(k) = sqrt ( s(k)*s(k+1) ) 9100 continue sh(nk) = sqrt ( s(nk) ) c c Compute Corilis parameter c f = 2.*7.29e-05*sinlat c c Define moisture constants c call WETCNS c c Define the grid point position (not necessary for column model) c ic = 6 jc = 7 c c Call radiation routine c call RAD (rtrsol,rtrter,flxvis,flxiri,ts,t,q,ps,sh,nk, & rlong,sinlat,coslat,s,hz,albedo,vent,declsc, & swr,swc,irtop,ic,jc) c c Write output of shortwave and longwave heating at surface(Watt/m2) c and vertical profiles of the shortwave and longwave fluxes(deg/day) c fact = 86400. 1000 format(2x,2f12.4) print *,flxvis,flxiri do 9102 k = 1, nk write(20,1000)rtrsol(k)*fact,rtrter(k)*fact 9102 continue stop end subroutine WETCNS c c This subroutine computes latent heat c of evaporation with respect to water c and ice . c common /wet/eps,ae(2),be(2),ht(2),tc data cp,elv,eli/1004.5,2.500e+6,2.834e+6/ data eps,e0,t0,rgocp/.622,610.71,273.,.286/ tc = 263. ht(1) = elv/cp ht(2) = eli / cp do 8110 i = 1, 2 be(i) = eps * ht(i) / rgocp 8110 ae(i) = be(i) / t0 + alog(e0) return end subroutine RAD (warm,cool,fsg,flg,tg,t,q,ps,sh,nk, & rlong,sinl,cosl,s,hz,albedo,vent,declsc, & swr,swc,irtop,ic,jc) c c This subroutine computes the shortwave and longwave c radiation . c parameter (ilgv = 17,ilatv = 10,ilevv = 12) parameter (ilev2v = ilevv+2) logical swr,swc real irtop(ilgv,ilatv) real t(nk),q(nk),sh(nk),declsc(2),pz(ilevv),albsc(8) & ,s(nk),warm(nk),cool(nk) common /rrr/fup(ilev2v,8),fdw(ilev2v,8),rw(ilev2v) & ,f(ilev2v),tl(ilev2v),aso(ilev2v,8),rso(ilev2v,8) & ,albc(8),rsa,secz,coe(8),dtc(ilev2v,8),fnet(ilev2v,8) & ,cld(ilev2v,8),rh(ilevv),dp(ilevv),qz(ilevv),p(ilev2v) & ,tz(ilevv),dpz(ilev2v),qs(ilevv),sso(8),dtw(ilev2v,8) data pclb,pclt,pcmb,pcmt,pchb,pcht/900.,2*700., & 2*500.,300.0/ data albc /0.,2*0.5,0.6,0.5,0.46,0.5,0.46 / data albsc /0.,2*0.66,0.76,0.66,0.5,0.66,0.54/ data el,g,sigma,pi,cp /0.75,9.81,5.667e-8,3.1415927,1004.5 / data iclr,iccl,iclm,iclmh,iccm,icch,icmh,iclh /1,2,3,4,5,6,7,8/ data rhcl,rhcm,rhch /0.66,0.5,0.4/ c c rwtop and ttop represent rw and t at the top wind level c in this case, the top level is chosen as the 100 mb level c data rwtop,ttop /4.7e-5,190.0/ el = 0.0 do 9100 k = 1, nk do 9100 jr= 1, 8 cld(k,jr) = 1. 9100 continue do 9102 k = 1, nk warm(k)= 0. cool(k)= 0. 9102 continue if(.not.swr) return c c Define some onstants needed. rs is the solar constant c and is set to 1.997 . psurf is the surface pressure in c mks and is to be multiplied by .01 for mb units . c rs = 1.997*4.186e4/60. romega = 360./24. albc(1) = albedo np = nk+1 n1 = nk ginv = 1./g picon = pi/180. gam = g/cp psurf = ps*0.01 c c Invert the fields so that the surface is level 1 and c top is level is level nk . c do 9104 i = 1, nk k = nk-i+1 tz(k) = t(i) qz(k) = q(i) pz(k) = sh(i)*psurf p(k+1) = s(i)*psurf 9104 continue p(1) = psurf do 9106 i = 2, nk tl(i) = tz(i-1)+(tz(i)-tz(i-1))*(alog(p(i))- & alog(pz(i-1)))/(alog(pz(i))-alog(pz(i-1))) c dpz(i) = pz(i-1)-pz(i) dp(i-1) = p(i-1)-p(i) 9106 continue dpz(1) = p(1) - p(2) dpz(np) = pz(nk) - p(np) dp(nk) = p(nk)-p(np) ts = tz(1)+(tz(1)-tz(2))*(alog(psurf)- & alog(pz(1)))/(alog(pz(1))-alog(pz(2))) tl(1) = tg * (1. - el) + ts * el tl(np) = ttop rw(np) = rwtop f(np) = sigma*tl(np)**4 do 9108 ir= 1, nk k = nk-ir+1 rw(k) = rw(k+1)+qz(k)*dp(k)*10.*ginv*(pz(k)* & .001)**.85*(ts/tz(k))**0.5 f(k) = sigma*tz(k)**4 rqs = humspc(tz(k),tz(k),pz(k)*100.) rh(k) = qz(k)/rqs if (rh(k).lt.1.0) go to 9108 rh(k) = 1.0 qs(k) = rqs 9108 continue icl = 1 iclp = 1 icm = 1 icmp = 1 ich = 1 ichp = 1 cl = 0. cm = 0. ch = 0. if (.not.swc) go to 72 cloud base and top for all cloud types rhl = 0. rhm = 0. rhh = 0. icl = n1 icm = n1 ich = n1 iclp = 0 icmp = 0 ichp = 0 dpl = 0. dpm = 0. dph = 0. do 9110 i = 1, nk if( pz(i).le.pclb.and.pz(i).gt.pclt) goto 22 goto 23 22 if (rh(i).ge.rhcl) goto 25 goto 9110 25 icl = min0(i,icl) iclp = i+1 dpl = dpl + dp(i) rhl = rhl + rh(i) * dp(i) go to 9110 23 if (pz(i).le.pcmb .and. pz(i).gt. pcmt) goto 26 goto 27 26 if (rh(i).ge.rhcm) goto 28 goto 9110 28 icm = min0(i,icm) icmp = i+1 dpm = dpm + dp(i) rhm = rhm + rh(i) * dp(i) go to 9110 27 if (pz(i).le.pchb .and.pz(i).gt.pcht) goto 29 goto 9110 29 if (rh(i).ge.rhch) goto 31 goto 9110 31 ich = min0(i,ich) ichp = i+1 dph = dph + dp(i) rhh = rhh + rh(i) * dp(i) 9110 continue c print 515,icl,iclp,icm,icmp,ich,ichp 515 format(1x,'icl,iclp,icm,icmp,ich,ichp',6i5) c c Calculate cloud amount for 8 configurations c do 9112 jr= 1, 8 do 9112 k = 1, np 9112 cld(k,jr) = 1. if (icl.gt.iclp) go to 50 cl = (rhl/dpl-rhcl)/(1.-rhcl) cl = cl*cl if (cl.ge.1.) cl = 1. iclp1 = iclp-1 do 9114 k = icl,iclp1 cld(k,8) = 0. do 9114 l = 2,4 cld(k,l) = 0. 9114 continue 50 if (icm.gt.icmp) go to 60 cm = (rhm/dpm-rhcm)/(1.-rhcm) cm = cm*cm if (cm.ge.1.) cm = 1. icmp1 = icmp-1 do 9116 k = icm, icmp1 do 9116 l = 3, 7 if (l.ne.6)cld(k,l) = 0. 9116 continue 60 if (ich.gt.ichp) go to 70 ch = (rhh/dph-rhch)/(1.-rhch) ch = ch*ch ch = ch*0.4 if (ch.ge.1.) ch = 1. ichp1 = ichp-1 do 9118 k = ich,ichp1 do 9118 l = 4,8,2 cld(k,l) = 0. if(l.eq.6)cld(k,l+1) =0. 9118 continue 70 if (iclp .eq.0)icl = 1 if (iclp .eq.0)iclp = 1 if (icmp.eq.0) icm = 1 if (icmp.eq.0) icmp = 1 if (ichp.eq.0) ich = 1 if (ichp.eq.0) ichp = 1 c c Print low medium and high cloud fractions c c print 71,cl,cm,ch 71 format(1x,'cl,cm,ch',3f10.2) 72 cl1 = 1.0-cl cm1 = 1.0-cm ch1 = 1.0-ch clr = cl1*cm1*ch1 coe(1) = clr ccl = cl*cm1*ch1 coe(2) = ccl clm = cl*cm*ch1 coe(3) = clm clmh = cl*cm*ch coe(4) = clmh ccm = cm*cl1*ch1 coe(5) = ccm cch = ch*cm1*cl1 coe(6) = cch cmh = cm*ch*cl1 coe(7) = cmh clh = cl*ch*cm1 coe(8) = clh do 9120 k = 1, np do 9120 jr= 1, 8 fup(k,jr) = 0. fdw(k,jr) = 0. if (k.ne.np) dtc(k,jr) = 0.0 9120 continue c c Compute the downward flux of long wave radiation . c call RLW (iclr,nk,icl,iclp,icm,icmp,ich,ichp,0) if (clr.ne.0. ) call RLW (iclr,nk,1,1,1,1,1,1,1) if (ccl.ne.0. ) call RLW (iccl,nk,icl,iclp,1,1,1,1,1) if (clm.ne.0. ) call RLW (iclm,nk,icl,iclp,icm,icmp,1,1,1) if (ccm.ne.0. ) call RLW (iccm,nk,1,1,icm,icmp,1,1,1) if (clmh.ne.0.) call RLW (iclmh,nk,icl,iclp,icm,icmp,ich,ichp,1) if (cch.ne.0. ) call RLW (icch,nk,1,1,1,1,ich,ichp,1) if (cmh.ne.0. ) call RLW (icmh,nk,1,1,icm,icmp,ich,ichp,1) if (clh.ne.0. ) call RLW (iclh,nk,icl,iclp,1,1,ich,ichp,1) c c Compute the downward long wave flux at surface (flg). c flg = 0. do 9122 jr= 1, 8 fdg = fdw(1,jr) if (fdg.eq.0.) fdg = sigma*ts**4 flg = flg + fdg*coe(jr) 9122 continue c c Compute the solar heating cosz is the cosine of the zenith angle c hz = hour,rlong = long in radians sinlat,coslat,sind,cosd from main prog sind = declsc(1) cosd = declsc(2) ang = hz*romega*picon+rlong-pi cosz = sind*sinl+cosd*cosl*cos(ang) if (cosz.le..01) cosz = 0.01 secz = 1./cosz rss = 0.651*rs*cosz rsa = 0.349*rs*cosz do 9124 k = 1, np do 9124 jr= 1, 8 if (k.ne.np) dtw(k,jr) = 0.0 aso(k,jr) = 0. rso(k,jr) = 0. 9124 continue if (cosz.le.0.01)go to 120 if (clr.ne.0. )call SLR (iclr,nk,1,1,1,1,1,1,pz) if (ccl.ne.0. )call SLR (iccl,nk,icl,iclp,1,1,1,1,pz) if (cch.ne.0. )call SLR (icch,nk,1,1,1,1,ich,ichp,pz) if (clm.ne.0. )call SLR (iclm,nk,icl,iclp,icm,icmp,1,1,pz) if (clmh.ne.0.)call SLR (iclmh,nk,icl,iclp,icm,icmp,ich,ichp,pz) if (ccm.ne.0. )call SLR (iccm,nk,1,1,icm,icmp,1,1,pz) if (cmh.ne.0. )call SLR (icmh,nk,1,1,icm,icmp,ich,ichp,pz) if (clh.ne.0. )call SLR (iclh,nk,icl,iclp,1,1,ich,ichp,pz) do 9126 jr= 1, 8 do 9126 k = 2, np dtw(k-1,jr)= gam*1.e-2/dp(k-1)*((aso(k,jr) & -aso(k-1,jr))+(rso(k-1,jr)-rso(k,jr))) 9126 continue 120 do 9130 k = 1, n1 ii = nk-k+1 warm(ii) = 0. do 9128 jr= 1, 8 warm(ii) = warm(ii)+dtw(k,jr)*coe(jr) 9128 continue 9130 continue fsga = 0. do 9132 jr= 1, 8 fsga = fsga +aso(1,jr)*coe(jr) 9132 continue resg = humspc(ts,ts,ps) c c Compute specific humidity at surface qsurf c assuming neutral mixed layer. c qss = qz(1) gw = 2.94*qss/resg-1.94 if (gw.lt.0 ) gw = 0. if (gw.gt.1.) gw = 1. if (cosz.le..01) go to 135 albair = 0.085-0.245*alog10(cosz) if (albair.gt.1.) albair = 1. fsgs = 0. do 9134 jr= 1, 8 r1 = albsc(jr)*(1.-cld(icl ,jr)) r2 = albsc(jr)*(1.-cld(icm ,jr)) r3 = albsc(jr)*(1.-cld(ich ,jr)) r123 = 1.-(1.-r1)*(1.-r2)*(1.-r3)/(1. & -r1*r2-r2*r3-r3*r1+2.*r1*r2*r3) x = 1.-(1.-r123)*(1.-albair) sso(jr)= rss*(1.-x)/(1.-x*albedo) fsgs = fsgs+sso(jr)*coe(jr) 9134 continue fsg = fsgs+fsga 135 if (cosz.le..01) fsg = 0. c c Computation of upward flux of longwave radiation c if (clr.ne.0. ) call RLW (iclr,nk,1,1,1,1,1,1,2) if (ccl.ne.0. ) call RLW (iccl,nk,icl,iclp,1,1,1,1,2) if (clm.ne.0. ) call RLW (iclm,nk,icl,iclp,icm,icmp,1,1,2) if (clmh.ne.0.) call RLW (iclmh,nk,icl,iclp,icm,icmp,ich,ichp,2) if (ccm.ne.0. ) call RLW (iccm,nk,1,1,icm,icmp,1,1,2) if (cmh.ne.0. ) call RLW (icmh,nk,1,1,icm,icmp,ich,ichp,2) if (cch.ne.0. ) call RLW (icch,nk,1,1,1,1,ich,ichp,2) if (clh.ne.0. ) call RLW (iclh,nk,icl,iclp,1,1,ich,ichp,2) do 9138 k = 1, np do 9136 jr= 1, 8 fnet(k,jr)= fup(k,jr)+fdw(k,jr) 9136 continue 9138 continue fuptop = 0. do 9140 jr= 1, 8 9140 fuptop = fuptop-fup(np,jr)*coe(jr) irtop(ic,jc) = fuptop do 9142 k = 2, np k1 = k-1 do 9142 jr= 1, 8 dtc(k1,jr)= -gam*1.e-2/dp(k1)*(fnet(k1,jr) & -fnet(k,jr))*cld(k1,jr) 9142 continue do 9146 k = 1, nk ii = nk-k+1 cool(ii) = 0. do 9144 jr= 1, 8 cool(ii) = cool(ii)+dtc(k,jr)*coe(jr) 9144 continue 9146 continue return end subroutine CONRAD (declsc,day,month) c c This subroutine computes the declination angle . c The Julian month is assumed to have length of c 365/12=30.4 c real declsc(2) data pi /3.14159/ ajour = 30.4 * (month-1) + day rdecl = .412 *cos((ajour+10.)*2.*pi /365.25 - pi) declsc(1) = sin (rdecl) declsc(2) = cos (rdecl) return end subroutine SLR (icat,nk,iclb,iclt,icmb,icmt,ichb,icht,pz) c c This subroutine computes the shortwave radiation. c parameter (ilevv = 12,ilev2v = ilevv+2) common /rrr/ fup(ilev2v,8),fdw(ilev2v,8),rw(ilev2v), & f(ilev2v),tl(ilev2v),aso(ilev2v,8),rso(ilev2v,8), & albc(8),rsa,secz,coe(8),dtc(ilev2v,8),fnet(ilev2v,8), & cld(ilev2v,8),rh(ilevv),dp(ilevv),qz(ilevv),p(ilev2v), & tz(ilevv),dpz(ilev2v),qs(ilevv),sso(8),dtw(ilev2v,8) real pz(1) c c Definition of variables . c c index 1 with high cloud c index 2 no high cloud c index 3 no high or middle clouds c aso direct solar insolation c rso reflected solar insolation c 1.4 equivalent optical depth of c a cloud of 200 mb thickness c n = nk+1 icb1 = n itop = n do 9120 i = 1, n rso(i,icat)= 0. 9120 continue aso(itop,icat) = rsa c c picks up the cloud top . c ict1 = max0(iclt,icmt,icht) index = 1 if (ict1.eq.icmt) index = 2 if (ict1.eq.iclt) index = 3 c c Computation for levels outside the cloud : c first above ,second below and third in the cloud . c 100 do 9121 i = ict1, icb1 if (i .eq. icb1 .and. icb1 .ne. n) goto 9121 c c clear case is eliminated. c dw = rw(i)*secz if (itop.ne.n)dw = 1.66*(rw(i)-rw(itop))+1.4*(p(icb1)-p(itop))/200. aso(i,icat) = (aso(itop,icat)-rso(itop,icat))*(1.-fcn(dw)) 9121 continue reaso = albc(icat)*aso(ict1,icat) c c This is alpha*s(absorbed)(1-a(wct sec(zeta)) c if (ict1.eq.1) reaso = albc(ict1)*aso(ict1,icat) do 9122 i = ict1, icb1 dw = 1.66*(rw(ict1)-rw(i)) rso(i,icat) = reaso*(1.-fcn(dw)) c c alpha c * sa *(1-a(wctsec(zeta))((1.66(wct-wi)) c 9122 continue if (ict1 .eq.1) return if (index-2) 10, 11, 12 10 icb1 = max0(iclb,icmb,ichb) if (ichb.eq.icmt) icb1 = icmb if (ichb.eq.icmt) index = index + 1 if (ichb.eq.icmt.and.icmb.eq.iclt) icb1 = iclb if (ichb.eq.icmt.and.icmb.eq.iclt) index = index + 1 go to 3 11 icb1 = max0(iclb,icmb) if (iclt.eq.icmb) icb1 = iclb if (iclt.eq.icmb) index = index + 1 go to 3 12 icb1 = iclb c c Perform levels inside the clouds . c 3 do 9123 i = icb1, ict1 if (i.eq.ict1) go to 9123 eqcw = 1.4*(p(i)-p(ict1))/200. dw = eqcw+1.66*(rw(i)-rw(ict1)) aso(i,icat) = (1.-albc(icat))*aso(ict1,icat)*(1.-fcn(dw)) 9123 continue c c For computations below clouds use the following. c if (icb1.eq.ichb) ict1 = max0(1,iclt,icmt) if (icb1.eq.ichb) itop = icht if (icb1.eq.icmb) ict1 = max0(1,iclt) if (icb1.eq.icmb) itop = icmt if (icb1.eq.iclb) ict1 = 1 if (icb1.eq.iclb) itop = iclt index = index+1 go to 100 end subroutine RLW (icat,nk,iclb,iclt,icmb,icmt,ichb,icht,ntyp) c c This subroutine computes the longwave radiative fluxes . c c ntyp = 0 initialize call c 1 downward flux c 2 upward flux c parameter (ilevv =12,ilev2v=ilevv+2,ilev2sq =ilev2v*ilev2v) common/rrr/fup(ilev2v,8),fdw(ilev2v,8),rw(ilev2v), & f(ilev2v),tl(ilev2v),aso(ilev2v,8),rso(ilev2v,8), & albc(8),rsa,secz,coe(8),dtc(ilev2v,8),fnet(ilev2v,8), & cld(ilev2v,8),rh(ilevv),dp(ilevv),qz(ilevv),p(ilev2v), & tz(ilevv),dpz(ilev2v),qs(ilevv),sso(8),dtw(ilev2v,8) real ff(ilev2v), em(ilev2v,ilev2v), tem(ilev2v) real dw(ilev2sq),emw(ilev2sq) data sigma / 5.667e-8 / a(x) = sigma * x ** 4 if (ntyp-1) 100, 200, 300 100 n = nk + 1 do 9130 k = 1, n 9130 em(k,k)= 0. j = 0 do 9131 k1= 1, nk k11 = k1 + 1 do 9131 k2= k11, n j = j+1 9131 dw(j) = rw(k1)-rw(k2) c c Define the emissivity function c call EMTAB (dw,emw,j) j = 0 do 9132 k1= 1, n k11 = k1+1 do 9132 k2= k11,n j = j+1 9132 em(k1,k2) = emw(j) call EMTAB (rw,tem,n) do 9133 k = 2, nk 9133 ff(k) = 0. ff(1) = a(tl(1)) ff(n) = f(n) if (iclt.le.iclb) go to 150 ff(iclb) = a(tl(iclb)) ff(iclt) = a(tl(iclt)) 150 if (icmt.le.icmb) go to 160 ff(icmb) = a(tl(icmb)) ff(icmt) = a(tl(icmt)) 160 if (icht.le.ichb) return ff(ichb) = a(tl(ichb)) ff(icht) = a(tl(icht)) return 200 index = 1 kkk = 1 icb1 = max0(iclt,icmt,icht) if (icb1.eq.icmt) index = 2 if (icb1.eq.iclt) index = 3 ict1 = n fdw(ict1,icat) = ff(ict1) * tem(ict1) do 9135 k = icb1, nk fd = ff(ict1) * (tem(k) - em(k,ict1)) do 9134 i = k, nk 9134 fd = fd + f(i) * (em(k,i+1) - em(k,i)) 9135 fdw(k,icat)= fd go to 260 230 fdw (ict1,icat) = ff(ict1) ict1m = ict1 - 1 do 9137 k = icb1, ict1m fd = ff(ict1) * (1. - em(k,ict1)) do 9136 i = k, ict1m 9136 fd = fd + f(i) * (em(k,i+1) - em(k,i)) 9137 fdw(k,icat)= fd 260 if (icb1.eq.1) return 270 if (index - 2) 275, 280, 285 275 index = 2 if (ichb .le. icmt) goto 280 ict1 = max0 (iclb,icmb,ichb) icb1 = max0 (1,iclt,icmt) if (kkk.eq.1) go to 230 go to 310 280 index = 3 if (icmb.eq.iclt) go to 285 ict1 = max0(iclb,icmb) icb1 = max0(1,iclt) if (kkk.eq.1) goto 230 goto 310 285 ict1 = iclb icb1 = 1 if (kkk.eq.1) goto 230 goto 310 300 index = 1 kkk = 2 icb1 = max0(iclt,icmt,icht) if (icb1.eq.icmt) index = 2 if (icb1.eq.iclt) index = 3 if (icb1.eq.1 ) index = 4 ict1 = n 310 do 9139 k = icb1, ict1 fu = ff(icb1) * (1. - em(icb1,k)) k1 = k - 1 do 9138 i = icb1, k1 9138 fu = fu + f(i) * (em(i,k) - em(i+1,k)) 9139 fup(k,icat)= -fu go to 260 end function HUMSPC (td,t,p) c c This subroutine computes the specific c humidity as a function of dewpoint td c and pressure p depending on phase. c common/wet/eps,ae(2),be(2),ht(2),tc iph = 1 if (t.le.tc) iph = 2 e = exp(ae(iph) - be(iph)/td ) humspc = eps * e / p return end subroutine EMTAB (dw,emit,nlev) c c This subroutines returns the emissivity values c as function of the optical depth in logarithmic scale c parameter (ilevv=12,ilev2v=ilevv+2,ilev2sq=ilev2v*ilev2v) real eps(18),wlog(19),dw(nlev),emit(nlev) real dwx(ilev2sq),wd(ilev2sq) data eps/0.0186,0.0258,0.0411,0.0572,0.0781,0.114,0.146,0.183, & 0.236,0.277,0.319,0.374,0.417,0.462,0.529,0.59,0.63,0.65/ data wlog/-5.,-4.7,-4.3,-4.0,-3.7,-3.3,-3.0,-2.7,-2.3,-2.0, & -1.7,-1.3,-1.0,-0.7,-0.3,-0.0,0.3,0.7,1.0/ data n,n1/18,17/ do 9140 j = 1, nlev emit(j)= 0. if (dw(j).ge.0.) go to 15 print 1000,dw dw(j) = 0. 15 dwx(j) = dw(j) if (dwx(j).eq.0.) dwx(j) = 1. 9140 continue do 9141 j = 1, nlev 9141 wd(j) = alog10(dwx(j)) c do 9142 j = 1, nlev if (dw(j).eq. 0.) goto 9142 if (wd(j) .gt.wlog(1)) goto 2 1 emit(j)= eps(1) goto 9142 2 do 9143 i = 1, n1 k = i if (wd(j).gt.wlog(i).and.wd(j).le.wlog(i+1)) go to 4 9143 continue 4 demow = (eps(k+1)-eps(k))/(wlog(k+1)-wlog(k)) emit(j) = eps(k)+(wd(j)-wlog(k))*demow 9142 continue return 1000 format(1x,'negative dw',e16.8) end function FCN (x) fcn = 0.271*x**0.303 return end