Back to home page

MITgcm
 
 

    

File indexing completed on 2023-03-03 06:09:59 UTC

view on githubraw file Latest commit 06d4643e on 2023-01-18 18:18:37 UTC
a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
180c87c6d9 Andr*0002       subroutine swrio (nymd,nhms,bi,bj,ndswr,myid,istrip,npcs,
082e38725b Andr*0003      .        low_level,mid_level,im,jm,lm,
d2aaec7e02 Andr*0004      .        pz,plz,plze,dpres,pkht,pkz,tz,qz,oz,co2,
7d59649d1b Andr*0005      .        albvisdr,albvisdf,albnirdr,albnirdf,
2a3ae9099b Andr*0006      .        dtradsw,dtswclr,radswg,swgclr,
7d59649d1b Andr*0007      .        fdifpar,fdirpar,osr,osrclr,
082e38725b Andr*0008      .        ptop,nswcld,cldsw,cswmo,nswlz,swlz,
7d59649d1b Andr*0009      .        lpnt,imstturb,qliqave,fccave,landtype,xlats,xlons)
1662f365b2 Andr*0010 
                0011       implicit none
                0012 
                0013 c Input Variables
                0014 c ---------------
d2aaec7e02 Andr*0015       integer nymd,nhms,bi,bj,ndswr,myid,istrip,npcs
                0016       integer mid_level,low_level
62e9d9f553 Jean*0017       integer im,jm,lm
a456aa407c Andr*0018       _RL  ptop
                0019       _RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1),dpres(im,jm,lm)
                0020       _RL pkht(im,jm,lm+1),pkz(im,jm,lm)
                0021       _RL tz(im,jm,lm),qz(im,jm,lm)
                0022       _RL oz(im,jm,lm)
                0023       _RL co2
                0024       _RL albvisdr(im,jm),albvisdf(im,jm),albnirdr(im,jm)
                0025       _RL albnirdf(im,jm)
                0026       _RL radswg(im,jm),swgclr(im,jm),fdifpar(im,jm),fdirpar(im,jm)
daaf4a8d7f Andr*0027       _RL osr(im,jm),osrclr(im,jm),dtradsw(im,jm,lm)
                0028       _RL dtswclr(im,jm,lm)
62e9d9f553 Jean*0029       integer nswcld,nswlz
                0030       _RL cldsw(im,jm,lm),cswmo(im,jm,lm),swlz(im,jm,lm)
                0031       logical lpnt
                0032       integer imstturb
                0033       _RL qliqave(im,jm,lm),fccave(im,jm,lm)
                0034       integer landtype(im,jm)
a456aa407c Andr*0035       _RL xlats(im,jm),xlons(im,jm)
1662f365b2 Andr*0036 
                0037 c Local Variables
                0038 c ---------------
d2aaec7e02 Andr*0039       integer   i,j,L,nn,nsecf
4a402f223d Andr*0040       integer   ntmstp,nymd2,nhms2
a456aa407c Andr*0041       _RL      getcon,grav,cp,undef
                0042       _RL      ra,alf,reffw,reffi,tminv
1662f365b2 Andr*0043 
62e9d9f553 Jean*0044       parameter ( reffw = 10.0 )
                0045       parameter ( reffi = 65.0 )
1662f365b2 Andr*0046 
a456aa407c Andr*0047       _RL tdry(im,jm,lm)
                0048       _RL TEMP1(im,jm)
                0049       _RL TEMP2(im,jm)
                0050       _RL zenith (im,jm)
                0051       _RL cldtot (im,jm,lm)
                0052       _RL cldmxo (im,jm,lm)
                0053       _RL totcld (im,jm)
                0054       _RL cldlow (im,jm)
                0055       _RL cldmid (im,jm)
                0056       _RL cldhi  (im,jm)
                0057       _RL taulow (im,jm)
                0058       _RL taumid (im,jm)
                0059       _RL tauhi  (im,jm)
                0060       _RL tautype(im,jm,lm,3)
                0061       _RL tau(im,jm,lm)
62e9d9f553 Jean*0062       _RL albedo(im,jm)
a456aa407c Andr*0063 
                0064       _RL PK(ISTRIP,lm)
                0065       _RL qzl(ISTRIP,lm),CLRO(ISTRIP,lm)
                0066       _RL TZL(ISTRIP,lm)
                0067       _RL OZL(ISTRIP,lm)
                0068       _RL PLE(ISTRIP,lm+1)
                0069       _RL COSZ(ISTRIP)
                0070       _RL dpstrip(ISTRIP,lm)
                0071 
                0072       _RL albuvdr(ISTRIP),albuvdf(ISTRIP)
                0073       _RL albirdr(ISTRIP),albirdf(ISTRIP)
                0074       _RL difpar (ISTRIP),dirpar (ISTRIP)
                0075 
                0076       _RL fdirir(istrip),fdifir(istrip)
                0077       _RL fdiruv(istrip),fdifuv(istrip)
                0078 
                0079       _RL flux(istrip,lm+1)
                0080       _RL fluxclr(istrip,lm+1)
                0081       _RL dtsw(istrip,lm)
                0082       _RL dtswc(istrip,lm)
                0083 
                0084       _RL taul   (istrip,lm)
                0085       _RL reff   (istrip,lm,2)
                0086       _RL tauc   (istrip,lm,2)
                0087       _RL taua   (istrip,lm)
                0088       _RL tstrip (istrip)
8b150b4af9 Andr*0089 #ifdef ALLOW_DIAGNOSTICS
                0090       logical  diagnostics_is_on
                0091       external diagnostics_is_on
                0092       _RL tmpdiag(im,jm),tmpdiag2(im,jm)
                0093 #endif
d2aaec7e02 Andr*0094 
                0095       logical first
                0096       data first /.true./
1662f365b2 Andr*0097 
                0098 C **********************************************************************
                0099 C ****                       INITIALIZATION                         ****
                0100 C **********************************************************************
                0101 
                0102       grav  = getcon('GRAVITY')
                0103       cp    = getcon('CP')
                0104       undef = getcon('UNDEF')
                0105 
                0106       NTMSTP = nsecf(NDSWR)
                0107       TMINV  = 1./float(ntmstp)
                0108 
                0109 C Compute Temperature from Theta
                0110 C ------------------------------
                0111       do L=1,lm
                0112       do j=1,jm
                0113       do i=1,im
                0114       tdry(i,j,L) = tz(i,j,L)*pkz(i,j,L)
                0115       enddo
                0116       enddo
                0117       enddo
                0118 
9524fe64b5 Andr*0119       if (first .and. myid.eq.1 ) then
1662f365b2 Andr*0120       print *
                0121       print *,'Low-Level Clouds are Grouped between levels: ',
                0122      .         lm,' and ',low_level
                0123       print *,'Mid-Level Clouds are Grouped between levels: ',
                0124      .         low_level-1,' and ',mid_level
                0125       print *
                0126       first = .false.
                0127       endif
                0128 
                0129 C **********************************************************************
                0130 C ****             CALCULATE COSINE OF THE ZENITH ANGLE             ****
                0131 C **********************************************************************
                0132 
62e9d9f553 Jean*0133 #ifdef FIZHI_USE_FIXED_DAY
                0134       CALL ASTRO ( 20040321,  NHMS,  XLATS,XLONS, im*jm, TEMP1,RA )
                0135                    NYMD2 = NYMD
                0136                    NHMS2 = NHMS
                0137       CALL TICK  ( NYMD2,  NHMS2, NTMSTP )
                0138       CALL ASTRO ( 20040321,  NHMS2, XLATS,XLONS, im*jm, TEMP2,RA )
                0139 #else
d2aaec7e02 Andr*0140       CALL ASTRO ( NYMD,   NHMS,  XLATS,XLONS, im*jm, TEMP1,RA )
1662f365b2 Andr*0141                    NYMD2 = NYMD
                0142                    NHMS2 = NHMS
                0143       CALL TICK  ( NYMD2,  NHMS2, NTMSTP )
d2aaec7e02 Andr*0144       CALL ASTRO ( NYMD2,  NHMS2, XLATS,XLONS, im*jm, TEMP2,RA )
62e9d9f553 Jean*0145 #endif
1662f365b2 Andr*0146 
                0147       do j = 1,jm
                0148       do i = 1,im
                0149           zenith(I,j) = TEMP1(I,j) + TEMP2(I,j)
                0150       IF( zenith(I,j) .GT. 0.0 ) THEN
                0151           zenith(I,j) = (2./3.)*( zenith(I,j)-TEMP2(I,j)*TEMP1(I,j)
                0152      .                          / zenith(I,j) )
                0153       ENDIF
                0154       ENDDO
                0155       ENDDO
                0156 
                0157 C **********************************************************************
                0158 c ****        Compute Two-Dimension Total Cloud Fraction (0-1)      ****
                0159 C **********************************************************************
                0160 
                0161 c Initialize Clear Sky Probability for Low, Mid, and High Clouds
                0162 c --------------------------------------------------------------
                0163       do j =1,jm
                0164       do i =1,im
                0165       cldlow(i,j) = 0.0
                0166       cldmid(i,j) = 0.0
                0167       cldhi (i,j) = 0.0
                0168       enddo
                0169       enddo
                0170 
                0171 c Adjust cloud fractions and cloud liquid water due to moist turbulence
                0172 c ---------------------------------------------------------------------
                0173       if(imstturb.ne.0) then
                0174         do L =1,lm
                0175         do j =1,jm
                0176         do i =1,im
e7070f537c Cons*0177            cldtot(i,j,L)=min(1.0 _d 0,max(cldsw(i,j,L),fccave(i,j,L)/
                0178      $          imstturb))
                0179            cldmxo(i,j,L)=min(1.0 _d 0,cswmo(i,j,L))
7d59649d1b Andr*0180            swlz(i,j,L)=swlz(i,j,L)+qliqave(i,j,L)/imstturb
1662f365b2 Andr*0181         enddo
                0182         enddo
                0183         enddo
                0184       else
                0185         do L =1,lm
                0186         do j =1,jm
                0187         do i =1,im
32362b8fa8 Cons*0188          cldtot(i,j,L) =  min( 1.0 _d 0,cldsw(i,j,L) )
                0189          cldmxo(i,j,L) =  min( 1.0 _d 0,cswmo(i,j,L) )
1662f365b2 Andr*0190         enddo
                0191         enddo
                0192         enddo
                0193       endif
                0194 
                0195 c Compute 3-D Cloud Fractions
                0196 c ---------------------------
                0197       if( nswcld.ne.0 ) then
                0198       do L = 1,lm
                0199       do j = 1,jm
                0200       do i = 1,im
                0201 c Compute Low-Mid-High Maximum Overlap Cloud Fractions
                0202 c ----------------------------------------------------
                0203       if(     L.lt.mid_level ) then
                0204              cldhi (i,j) = max( cldtot(i,j,L),cldhi (i,j) )
62e9d9f553 Jean*0205       elseif( L.ge.mid_level .and.
1662f365b2 Andr*0206      .        L.lt.low_level ) then
                0207              cldmid(i,j) = max( cldtot(i,j,L),cldmid(i,j) )
                0208       elseif( L.ge.low_level ) then
                0209              cldlow(i,j) = max( cldtot(i,j,L),cldlow(i,j) )
                0210       endif
                0211 
                0212       enddo
                0213       enddo
                0214       enddo
                0215       endif
                0216 
                0217 c Totcld => Product of Clear Sky Probabilities for Low, Mid, and High Clouds
                0218 c --------------------------------------------------------------------------
                0219       do j = 1,jm
                0220       do i = 1,im
                0221       totcld(i,j) = 1.0 - (1.-cldhi (i,j))
                0222      .                  * (1.-cldmid(i,j))
                0223      .                  * (1.-cldlow(i,j))
                0224       enddo
                0225       enddo
                0226 
8b150b4af9 Andr*0227 #ifdef ALLOW_DIAGNOSTICS
                0228 
1662f365b2 Andr*0229 c Compute Cloud Diagnostics
                0230 c -------------------------
8b150b4af9 Andr*0231       if(diagnostics_is_on('CLDFRC  ',myid) ) then
                0232        call diagnostics_fill(totcld,'CLDFRC  ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0233       endif
                0234 
8b150b4af9 Andr*0235       if(diagnostics_is_on('CLDRAS  ',myid) ) then
                0236        do L=1,lm
                0237        do j=1,jm
                0238        do i=1,im
                0239         tmpdiag(i,j) = cswmo(i,j,L)
                0240        enddo
                0241        enddo
                0242        call diagnostics_fill(tmpdiag,'CLDRAS  ',L,1,3,bi,bj,myid)
1662f365b2 Andr*0243       enddo
                0244       endif
                0245 
8b150b4af9 Andr*0246       if(diagnostics_is_on('CLDTOT  ',myid) ) then
                0247        do L=1,lm
                0248        do j=1,jm
                0249        do i=1,im
                0250         tmpdiag(i,j) = cldtot(i,j,L)
                0251        enddo
                0252        enddo
                0253        call diagnostics_fill(tmpdiag,'CLDTOT  ',L,1,3,bi,bj,myid)
1662f365b2 Andr*0254       enddo
                0255       endif
                0256 
8b150b4af9 Andr*0257       if(diagnostics_is_on('CLDLOW  ',myid) ) then
                0258        call diagnostics_fill(cldlow,'CLDLOW  ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0259       endif
                0260 
8b150b4af9 Andr*0261       if(diagnostics_is_on('CLDMID  ',myid) ) then
                0262        call diagnostics_fill(cldmid,'CLDMID  ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0263       endif
                0264 
8b150b4af9 Andr*0265       if(diagnostics_is_on('CLDHI   ',myid) ) then
                0266        call diagnostics_fill(cldhi,'CLDHI   ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0267       endif
                0268 
8b150b4af9 Andr*0269       if(diagnostics_is_on('LZRAD   ',myid) ) then
                0270        do L=1,lm
                0271        do j=1,jm
                0272        do i=1,im
                0273         tmpdiag(i,j) = swlz(i,j,L) * 1.0e6
                0274        enddo
                0275        enddo
                0276        call diagnostics_fill(tmpdiag,'LZRAD   ',L,1,3,bi,bj,myid)
                0277        enddo
1662f365b2 Andr*0278       endif
                0279 
                0280 c Albedo Diagnostics
                0281 c ------------------
8b150b4af9 Andr*0282       if(diagnostics_is_on('ALBVISDR',myid) ) then
                0283        call diagnostics_fill(albvisdr,'ALBVISDR',0,1,3,bi,bj,myid)
1662f365b2 Andr*0284       endif
                0285 
8b150b4af9 Andr*0286       if(diagnostics_is_on('ALBVISDF',myid) ) then
                0287        call diagnostics_fill(albvisdf,'ALBVISDF',0,1,3,bi,bj,myid)
1662f365b2 Andr*0288       endif
                0289 
8b150b4af9 Andr*0290       if(diagnostics_is_on('ALBNIRDR',myid) ) then
                0291        call diagnostics_fill(albnirdr,'ALBNIRDR',0,1,3,bi,bj,myid)
1662f365b2 Andr*0292       endif
                0293 
8b150b4af9 Andr*0294       if(diagnostics_is_on('ALBNIRDF',myid) ) then
                0295        call diagnostics_fill(albnirdf,'ALBNIRDF',0,1,3,bi,bj,myid)
1662f365b2 Andr*0296       endif
                0297 
8b150b4af9 Andr*0298 #endif
                0299 
1662f365b2 Andr*0300 C Compute Optical Thicknesses and Diagnostics
                0301 C -------------------------------------------
7d59649d1b Andr*0302       call opthk(tdry,plz,plze,swlz,cldtot,cldmxo,landtype,im,jm,lm,
                0303      .                                                          tautype)
1662f365b2 Andr*0304 
                0305       do L = 1,lm
                0306       do j = 1,jm
                0307       do i = 1,im
7d59649d1b Andr*0308       tau(i,j,L) = tautype(i,j,L,1)+tautype(i,j,L,2)+tautype(i,j,L,3)
1662f365b2 Andr*0309       enddo
                0310       enddo
                0311       enddo
                0312 
8b150b4af9 Andr*0313 #ifdef ALLOW_DIAGNOSTICS
                0314       if(diagnostics_is_on('TAUAVE  ',myid) ) then
                0315        do L=1,lm
                0316        do j=1,jm
                0317        do i=1,im
                0318         tmpdiag(i,j) = tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L))
                0319        enddo
                0320        enddo
                0321        call diagnostics_fill(tmpdiag,'TAUAVE  ',L,1,3,bi,bj,myid)
1662f365b2 Andr*0322       enddo
                0323       endif
                0324 
8b150b4af9 Andr*0325       if(diagnostics_is_on('TAUCLD  ',myid) .and.
                0326      .                 diagnostics_is_on('TAUCLDC ',myid) ) then
                0327        do L=1,lm
                0328        do j=1,jm
                0329        do i=1,im
                0330         if( cldtot(i,j,L).ne.0.0 ) then
                0331         tmpdiag(i,j) = tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L))
                0332         tmpdiag2(i,j) = 1.
                0333         else
                0334         tmpdiag(i,j) = 0.
                0335         tmpdiag2(i,j) = 0.
                0336         endif
                0337        enddo
                0338        enddo
                0339        call diagnostics_fill(tmpdiag,'TAUCLD  ',L,1,3,bi,bj,myid)
                0340        call diagnostics_fill(tmpdiag,'TAUCLDC ',L,1,3,bi,bj,myid)
1662f365b2 Andr*0341       enddo
                0342       endif
                0343 
                0344 c Compute Low, Mid, and High Cloud Optical Depth Diagnostics
                0345 c ----------------------------------------------------------
8b150b4af9 Andr*0346       if(diagnostics_is_on('TAULOW  ',myid) .and.
                0347      .                 diagnostics_is_on('TAULOWC ',myid) ) then
180c87c6d9 Andr*0348        do j = 1,jm
                0349        do i = 1,im
8b150b4af9 Andr*0350         taulow(i,j) =  0.0
180c87c6d9 Andr*0351         if( cldlow(i,j).ne.0.0 ) then
                0352          do L = low_level,lm
                0353           taulow(i,j) = taulow(i,j) + tau(i,j,L)
                0354          enddo
8b150b4af9 Andr*0355          tmpdiag2(i,j) = 1.
                0356         else
                0357          tmpdiag(i,j) = 0.
180c87c6d9 Andr*0358         endif
                0359        enddo
                0360        enddo
8b150b4af9 Andr*0361        call diagnostics_fill(taulow,'TAULOW  ',0,1,3,bi,bj,myid)
                0362        call diagnostics_fill(tmpdiag2,'TAULOWC ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0363       endif
                0364 
8b150b4af9 Andr*0365       if(diagnostics_is_on('TAUMID  ',myid) .and.
                0366      .                 diagnostics_is_on('TAUMIDC ',myid) ) then
180c87c6d9 Andr*0367        do j = 1,jm
                0368        do i = 1,im
8b150b4af9 Andr*0369         taumid(i,j) =  0.0
180c87c6d9 Andr*0370         if( cldmid(i,j).ne.0.0 ) then
                0371          do L = mid_level,low_level+1
                0372           taumid(i,j) = taumid(i,j) + tau(i,j,L)
                0373          enddo
8b150b4af9 Andr*0374          tmpdiag2(i,j) = 1.
                0375         else
                0376          tmpdiag(i,j) = 0.
180c87c6d9 Andr*0377         endif
                0378        enddo
                0379        enddo
8b150b4af9 Andr*0380        call diagnostics_fill(taumid,'TAUMID  ',0,1,3,bi,bj,myid)
                0381        call diagnostics_fill(tmpdiag2,'TAUMIDC ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0382       endif
                0383 
8b150b4af9 Andr*0384       if(diagnostics_is_on('TAUHI   ',myid) .and.
                0385      .                 diagnostics_is_on('TAUHIC  ',myid) ) then
180c87c6d9 Andr*0386        do j = 1,jm
                0387        do i = 1,im
8b150b4af9 Andr*0388         tauhi(i,j) =  0.0
180c87c6d9 Andr*0389         if( cldhi(i,j).ne.0.0 ) then
                0390          do L = 1,mid_level+1
                0391           tauhi(i,j) = tauhi(i,j) + tau(i,j,L)
                0392          enddo
8b150b4af9 Andr*0393          tmpdiag2(i,j) = 1.
                0394         else
                0395          tmpdiag(i,j) = 0.
180c87c6d9 Andr*0396         endif
                0397        enddo
                0398        enddo
8b150b4af9 Andr*0399        call diagnostics_fill(tauhi,'TAUHI   ',0,1,3,bi,bj,myid)
                0400        call diagnostics_fill(tmpdiag2,'TAUHIC  ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0401       endif
8b150b4af9 Andr*0402 #endif
1662f365b2 Andr*0403 
                0404 C***********************************************************************
                0405 C ****                     LOOP OVER GLOBAL REGIONS                 ****
                0406 C **********************************************************************
                0407 
                0408       do 1000 nn = 1,npcs
                0409 
                0410 C **********************************************************************
                0411 C ****                   VARIABLE INITIALIZATION                    ****
                0412 C **********************************************************************
                0413 
                0414       CALL STRIP ( zenith,COSZ,im*jm,ISTRIP,1,NN )
                0415 
d2aaec7e02 Andr*0416       CALL STRIP ( plze,  ple   ,im*jm,ISTRIP,lm+1,NN)
                0417       CALL STRIP ( pkz ,  pk    ,im*jm,ISTRIP,lm  ,NN)
                0418       CALL STRIP ( dpres,dpstrip,im*jm,ISTRIP,lm  ,NN)
                0419       CALL STRIP ( tdry,  tzl   ,im*jm,ISTRIP,lm  ,NN)
                0420       CALL STRIP ( qz  ,  qzl   ,im*jm,ISTRIP,lm  ,NN)
                0421       CALL STRIP ( oz  ,  ozl   ,im*jm,ISTRIP,lm  ,NN)
                0422       CALL STRIP ( tau ,  taul  ,im*jm,ISTRIP,lm  ,NN)
1662f365b2 Andr*0423 
                0424       CALL STRIP ( albvisdr,albuvdr,im*jm,ISTRIP,1,NN )
                0425       CALL STRIP ( albvisdf,albuvdf,im*jm,ISTRIP,1,NN )
                0426       CALL STRIP ( albnirdr,albirdr,im*jm,ISTRIP,1,NN )
                0427       CALL STRIP ( albnirdf,albirdf,im*jm,ISTRIP,1,NN )
                0428 
                0429       call strip ( cldtot,clro,im*jm,istrip,lm,nn )
                0430 
                0431 c Partition Tau between Water and Ice Particles
                0432 c ---------------------------------------------
                0433       do L= 1,lm
                0434       do i= 1,istrip
32362b8fa8 Cons*0435               alf = min( max((tzl(i,l)-253.15)/20.,0. _d 0) ,1. _d 0)
1662f365b2 Andr*0436       taua(i,L)   = 0.
                0437 
                0438       if( alf.ne.0.0 .and. alf.ne.1.0 ) then
                0439       tauc(i,L,1) = taul(i,L)/(1.+alf/(1-alf) * (reffi/reffw*0.80) )
                0440       tauc(i,L,2) = taul(i,L)/(1.+(1-alf)/alf * (reffw/reffi*1.25) )
                0441       endif
                0442 
                0443       if( alf.eq.0.0 ) then
                0444       tauc(i,L,1) = taul(i,L)
                0445       tauc(i,L,2) = 0.0
                0446       endif
                0447 
                0448       if( alf.eq.1.0 ) then
                0449       tauc(i,L,1) = 0.0
                0450       tauc(i,L,2) = taul(i,L)
                0451       endif
                0452 
                0453       reff(i,L,1) = reffi
                0454       reff(i,L,2) = reffw
                0455       enddo
                0456       enddo
                0457 
                0458       call sorad ( istrip,1,1,lm,ple,tzl,qzl,ozl,co2,
                0459      .             tauc,reff,clro,mid_level,low_level,taua,
                0460      .             albirdr,albirdf,albuvdr,albuvdf,cosz,
                0461      .             flux,fluxclr,fdirir,fdifir,dirpar,difpar,
                0462      .             fdiruv,fdifuv )
                0463 
                0464 C **********************************************************************
                0465 C ****     Compute Mass-Weighted Theta Tendencies from Fluxes       ****
                0466 C **********************************************************************
                0467 
                0468       do l=1,lm
                0469       do i=1,istrip
032fb71841 Andr*0470       alf = grav*(ple(i,Lm+1)-ptop)/(cp*dpstrip(i,L)*100)
1662f365b2 Andr*0471       dtsw (i,L) = alf*( flux   (i,L)-flux   (i,L+1) )/pk(i,L)
                0472       dtswc(i,L) = alf*( fluxclr(i,L)-fluxclr(i,L+1) )/pk(i,L)
                0473       enddo
                0474       enddo
62e9d9f553 Jean*0475 
1662f365b2 Andr*0476       call paste ( dtsw , dtradsw ,istrip,im*jm,lm,nn )
                0477       call paste ( dtswc, dtswclr ,istrip,im*jm,lm,nn )
                0478 
                0479       call paste ( flux   (1,1),osr   ,istrip,im*jm,1,nn )
                0480       call paste ( fluxclr(1,1),osrclr,istrip,im*jm,1,nn )
                0481 
                0482       call paste ( flux   (1,lm+1),radswg,istrip,im*jm,1,nn )
                0483       call paste ( fluxclr(1,lm+1),swgclr,istrip,im*jm,1,nn )
                0484 
                0485       call paste ( dirpar,fdirpar,istrip,im*jm,1,nn )
                0486       call paste ( difpar,fdifpar,istrip,im*jm,1,nn )
                0487 
                0488 c Calculate Mean Albedo
                0489 c ---------------------
                0490       do i=1,istrip
7d59649d1b Andr*0491        if( cosz(i).gt.0.0 ) then
                0492         tstrip(i) = 1.0 - flux(i,lm+1)/
                0493      . ( fdirir(i)+fdifir(i)+dirpar(i)+difpar(i) + fdiruv(i)+fdifuv(i) )
1662f365b2 Andr*0494         if( tstrip(i).lt.0.0 ) tstrip(i) = undef
7d59649d1b Andr*0495        else
                0496         tstrip(i) = undef
                0497        endif
1662f365b2 Andr*0498       enddo
                0499       call paste ( tstrip,albedo,istrip,im*jm,1,nn )
                0500 
                0501  1000 CONTINUE
                0502 
8b150b4af9 Andr*0503 #ifdef ALLOW_DIAGNOSTICS
                0504 
1662f365b2 Andr*0505 c Mean Albedo Diagnostic
                0506 c ----------------------
8b150b4af9 Andr*0507       if(diagnostics_is_on('ALBEDO  ',myid) .and.
                0508      .                 diagnostics_is_on('ALBEDOC ',myid) ) then
                0509        do j=1,jm
                0510        do i=1,im
                0511         if( albedo(i,j).ne.undef ) then
                0512          tmpdiag(i,j) = albedo(i,j)
                0513          tmpdiag2(i,j) = 1.
                0514         else
                0515          tmpdiag(i,j) = 0.
                0516          tmpdiag2(i,j) = 0.
                0517         endif
                0518        enddo
                0519        enddo
                0520        call diagnostics_fill(tmpdiag,'ALBEDO  ',0,1,3,bi,bj,myid)
                0521        call diagnostics_fill(tmpdiag2,'ALBEDOC ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0522       endif
8b150b4af9 Andr*0523 #endif
1662f365b2 Andr*0524 
                0525 C **********************************************************************
                0526 C ****                 ZERO-OUT OR BUMP COUNTERS                    ****
                0527 C **********************************************************************
                0528 
                0529       nswlz    = 0
                0530       nswcld   = 0
                0531       imstturb = 0
                0532 
                0533       do L = 1,lm
                0534       do j = 1,jm
                0535       do i = 1,im
                0536          fccave(i,j,L) = 0.
                0537         qliqave(i,j,L) = 0.
                0538       enddo
                0539       enddo
                0540       enddo
                0541 
                0542       return
                0543       end
                0544       subroutine opthk ( tl,pl,ple,lz,cf,cfm,lwi,im,jm,lm,tau )
                0545 C***********************************************************************
                0546 C
                0547 C  PURPOSE:
                0548 C  ========
                0549 C    Compute cloud optical thickness using temperature and
                0550 C    cloud fractions.
                0551 C
                0552 C  INPUT:
                0553 C  ======
                0554 C    tl ......... Temperature at Model Layers (K)
                0555 C    pl ......... Model Layer Pressures (mb)
                0556 C    ple ........ Model Edge  Pressures (mb)
                0557 C    lz ......... Diagnosed Convective and Large-Scale Cloud Water (g/g)
                0558 C    cf ......... Total Cloud Fraction (Random + Maximum Overlap)
                0559 C    cfm ........ Maximum Overlap Cloud Fraction
                0560 C    lwi ........ Surface Land Type
                0561 C    im ......... Longitudinal dimension
                0562 C    jm ......... Latitudinal  dimension
                0563 C    lm ......... Vertical     dimension
                0564 C
                0565 C  OUTPUT:
                0566 C  =======
                0567 C    tau ........ Cloud Optical Thickness (non-dimensional)
                0568 C                  tau(im,jm,lm,1):  Suspended Ice
                0569 C                  tau(im,jm,lm,2):  Suspended Water
                0570 C                  tau(im,jm,lm,3):  Raindrops
                0571 C
                0572 C***********************************************************************
                0573 
                0574       implicit none
                0575 
                0576       integer  im,jm,lm,i,j,L
                0577 
a456aa407c Andr*0578       _RL  tl(im,jm,lm)
                0579       _RL  pl(im,jm,lm)
                0580       _RL ple(im,jm,lm+1)
                0581       _RL  lz(im,jm,lm)
                0582       _RL  cf(im,jm,lm)
                0583       _RL cfm(im,jm,lm)
                0584       _RL tau(im,jm,lm,3)
1662f365b2 Andr*0585       integer lwi(im,jm)
                0586 
a456aa407c Andr*0587       _RL dp, alf, fracls, fraccu
                0588       _RL tauice, tauh2o, tauras
1662f365b2 Andr*0589 
                0590 c Compute Cloud Optical Depths
                0591 c ----------------------------
                0592       do L=1,lm
                0593       do j=1,jm
                0594       do i=1,im
e7070f537c Cons*0595          alf   =  min( max(( tl(i,j,L)-233.15)/20.,0. _d 0), 1. _d 0)
                0596          dp   =  ple(i,j,L+1)-ple(i,j,L)
1662f365b2 Andr*0597          tau(i,j,L,1)  = 0.0
                0598          tau(i,j,L,2)  = 0.0
                0599          tau(i,j,L,3)  = 0.0
                0600       if( cf(i,j,L).ne.0.0 ) then
                0601 
                0602 c Determine fraction of large-scale and cumulus clouds
                0603 c ----------------------------------------------------
                0604          fracls = ( cf(i,j,L)-cfm(i,j,L) )/cf(i,j,L)
                0605          fraccu = 1.0-fracls
                0606 
                0607 c Define tau for large-scale ice and water clouds
62e9d9f553 Jean*0608 c Note:  tauice is scaled between (0.02 & .2) for:  0 < lz <  2 mg/kg over Land
                0609 c Note:  tauh2o is scaled between (0.20 & 20) for:  0 < lz <  5 mg/kg over Land
1662f365b2 Andr*0610 c Note:  tauh2o is scaled between (0.20 & 12) for:  0 < lz < 50 mg/kg over Ocean
                0611 c -------------------------------------------------------------------------------
                0612 
                0613 c Large-Scale Ice
                0614 c ---------------
e7070f537c Cons*0615          tauice = max( 0.0002 _d 0, 0.002*min(500*lz(i,j,L)*1000,
                0616      $        1.0 _d 0) )
1662f365b2 Andr*0617          tau(i,j,L,1) = fracls*(1-alf)*tauice*dp
                0618 
                0619 c Large-Scale Water
                0620 c -----------------
7d59649d1b Andr*0621 C Over Land
1662f365b2 Andr*0622          if( lwi(i,j).le.10 ) then
e7070f537c Cons*0623             tauh2o = max( 0.0020 _d 0, 0.200*min(200*lz(i,j,L)*1000,
62e9d9f553 Jean*0624      $           1.0 _d 0) )
7d59649d1b Andr*0625           tau(i,j,L,3) = fracls*alf*tauh2o*dp
1662f365b2 Andr*0626          else
7d59649d1b Andr*0627 C Non-Precipitation Clouds Over Ocean
                0628           if( lz(i,j,L).eq.0.0 ) then
62e9d9f553 Jean*0629            tauh2o = .12
7d59649d1b Andr*0630            tau(i,j,L,2) = fracls*alf*tauh2o*dp
                0631           else
                0632 C Over Ocean
e7070f537c Cons*0633              tauh2o = max( 0.0020 _d 0, 0.120*min( 20*lz(i,j,L)*1000,
                0634      $            1.0 _d 0) )
7d59649d1b Andr*0635            tau(i,j,L,3) = fracls*alf*tauh2o*dp
                0636           endif
1662f365b2 Andr*0637          endif
                0638 
                0639 c Sub-Grid Convective
                0640 c -------------------
                0641          if( tl(i,j,L).gt.210.0 ) then
                0642          tauras = .16
                0643          tau(i,j,L,3) = tau(i,j,L,3) + fraccu*tauras*dp
                0644          else
                0645          tauras = tauice
                0646          tau(i,j,L,1) = tau(i,j,L,1) + fraccu*tauras*dp
                0647          endif
                0648 
                0649 c Define tau for large-scale and cumulus clouds
                0650 c ---------------------------------------------
                0651 ccc      tau(i,j,L) = ( fracls*( alf*tauh2o + (1-alf)*tauice )
                0652 ccc  .                + fraccu*tauras )*dp
                0653 
                0654       endif
                0655 
                0656       enddo
                0657       enddo
                0658       enddo
                0659 
                0660       return
                0661       end
                0662       subroutine sorad(m,n,ndim,np,pl,ta,wa,oa,co2,
                0663      *                 taucld,reff,fcld,ict,icb,
                0664      *                 taual,rsirbm,rsirdf,rsuvbm,rsuvdf,cosz,
                0665      *                 flx,flc,fdirir,fdifir,fdirpar,fdifpar,
                0666      *                 fdiruv,fdifuv)
                0667 c********************************************************************
62e9d9f553 Jean*0668 c
1662f365b2 Andr*0669 c This routine computes solar fluxes due to the absoption by water
                0670 c  vapor, ozone, co2, o2, clouds, and aerosols and due to the
                0671 c  scattering by clouds, aerosols, and gases.
                0672 c
                0673 c This is a vectorized code. It computes the fluxes simultaneous for
                0674 c  (m x n) soundings, which is a subset of the (m x ndim) soundings.
                0675 c  In a global climate model, m and ndim correspond to the numbers of
                0676 c  grid boxes in the zonal and meridional directions, respectively.
                0677 c
                0678 c Ice and liquid cloud particles are allowed to co-exist in any of the
                0679 c  np layers. Two sets of cloud parameters are required as inputs, one
                0680 c  for ice paticles and the other for liquid particles.  These parameters
                0681 c  are optical thickness (taucld) and effective particle size (reff).
                0682 c
                0683 c If no information is available for reff, a default value of
                0684 c  10 micron for liquid water and 75 micron for ice can be used.
                0685 c
                0686 c Clouds are grouped into high, middle, and low clouds separated by the
                0687 c  level indices ict and icb.  For detail, see the subroutine cldscale.
                0688 c
62e9d9f553 Jean*0689 c----- Input parameters:
1662f365b2 Andr*0690 c                                                   units      size
                0691 c   number of soundings in zonal direction (m)       n/d        1
                0692 c   number of soundings in meridional direction (n)  n/d        1
                0693 c   maximum number of soundings in                   n/d        1
                0694 c           meridional direction (ndim)
                0695 c   number of atmospheric layers (np)                n/d        1
                0696 c   level pressure (pl)                              mb       m*ndim*(np+1)
                0697 c   layer temperature (ta)                           k        m*ndim*np
                0698 c   layer specific humidity (wa)                     gm/gm    m*ndim*np
                0699 c   layer ozone concentration (oa)                   gm/gm    m*ndim*np
                0700 c   co2 mixing ratio by volumn (co2)               parts/part   1
                0701 c   cloud optical thickness (taucld)                 n/d      m*ndim*np*2
                0702 c                index 1 for ice particles
                0703 c                index 2 for liquid drops
                0704 c   effective cloud-particle size (reff)           micrometer m*ndim*np*2
                0705 c                index 1 for ice particles
                0706 c                index 2 for liquid drops
                0707 c   cloud amount (fcld)                            fraction   m*ndim*np
                0708 c   level index separating high and middle           n/d        1
                0709 c                clouds (ict)
                0710 c   level index separating middle and low clouds     n/d        1
                0711 c                clouds (icb)
62e9d9f553 Jean*0712 c   aerosol optical thickness (taual)                n/d      m*ndim*np
1662f365b2 Andr*0713 c   solar ir surface albedo for beam                fraction   m*ndim
62e9d9f553 Jean*0714 c                radiation (rsirbm)
1662f365b2 Andr*0715 c   solar ir surface albedo for diffuse             fraction   m*ndim
62e9d9f553 Jean*0716 c                radiation (rsirdf)
1662f365b2 Andr*0717 c   uv + par surface albedo for beam                     fraction   m*ndim
62e9d9f553 Jean*0718 c                radiation (rsuvbm)
1662f365b2 Andr*0719 c   uv + par surface albedo for diffuse                  fraction   m*ndim
62e9d9f553 Jean*0720 c                radiation (rsuvdf)
1662f365b2 Andr*0721 c   cosine of solar zenith angle (cosz)            n/d        m*ndim
                0722 c
                0723 c----- Output parameters
                0724 c
                0725 c   all-sky flux (downward minus upward) (flx)     fraction   m*ndim*(np+1)
                0726 c   clear-sky flux (downward minus upward) (flc)   fraction   m*ndim*(np+1)
                0727 c   all-sky direct downward ir (0.7-10 micron)
                0728 c                flux at the surface (fdirir)      fraction   m*ndim
                0729 c   all-sky diffuse downward ir flux at
                0730 c                the surface (fdifir)              fraction   m*ndim
                0731 c   all-sky direct downward par (0.4-0.7 micron)
                0732 c                flux at the surface (fdirpar)     fraction   m*ndim
                0733 c   all-sky diffuse downward par flux at
                0734 c                the surface (fdifpar)             fraction   m*ndim
                0735 *
                0736 c----- Notes:
                0737 c
                0738 c    (1) The unit of flux is fraction of the incoming solar radiation
                0739 c        at the top of the atmosphere.  Therefore, fluxes should
                0740 c        be equal to flux multiplied by the extra-terrestrial solar
                0741 c        flux and the cosine of solar zenith angle.
                0742 c    (2) Clouds and aerosols can be included in any layers by specifying
62e9d9f553 Jean*0743 c        fcld(i,j,k), taucld(i,j,k,*) and taual(i,j,k), k=1,np.
1662f365b2 Andr*0744 c        For an atmosphere without clouds and aerosols,
                0745 c        set fcld(i,j,k)=taucld(i,j,k,*)=taual(i,j,k)=0.0.
                0746 c    (3) Aerosol single scattering albedos and asymmetry
                0747 c        factors are specified in the subroutines solir and soluv.
                0748 c    (4) pl(i,j,1) is the pressure at the top of the model, and
                0749 c        pl(i,j,np+1) is the surface pressure.
                0750 c    (5) the pressure levels ict and icb correspond approximately
                0751 c        to 400 and 700 mb.
62e9d9f553 Jean*0752 c
1662f365b2 Andr*0753 c**************************************************************************
62e9d9f553 Jean*0754 
1662f365b2 Andr*0755       implicit none
                0756 
                0757 c-----Explicit Inline Directives
                0758 
06d4643e1f Jean*0759 #ifdef FIZHI_CRAY
4e1c053948 Jean*0760 #ifdef FIZHI_F77_COMPIL
1662f365b2 Andr*0761 cfpp$ expand (expmn)
                0762 #endif
                0763 #endif
a456aa407c Andr*0764       _RL expmn
1662f365b2 Andr*0765 
                0766 c-----input parameters
                0767 
                0768       integer m,n,ndim,np,ict,icb
a456aa407c Andr*0769       _RL pl(m,ndim,np+1),ta(m,ndim,np),wa(m,ndim,np),oa(m,ndim,np)
                0770       _RL  taucld(m,ndim,np,2),reff(m,ndim,np,2)
                0771       _RL  fcld(m,ndim,np),taual(m,ndim,np)
                0772       _RL  rsirbm(m,ndim),rsirdf(m,ndim),
1662f365b2 Andr*0773      *     rsuvbm(m,ndim),rsuvdf(m,ndim),cosz(m,ndim),co2
                0774 
                0775 c-----output parameters
                0776 
a456aa407c Andr*0777       _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
                0778       _RL  fdirir(m,ndim),fdifir(m,ndim)
                0779       _RL  fdirpar(m,ndim),fdifpar(m,ndim)
                0780       _RL  fdiruv(m,ndim),fdifuv(m,ndim)
1662f365b2 Andr*0781 
                0782 c-----temporary array
62e9d9f553 Jean*0783 
175684e43e Andr*0784       integer i,j,k
a456aa407c Andr*0785       _RL  cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np)
                0786       _RL  dp(m,n,np),wh(m,n,np),oh(m,n,np),scal(m,n,np)
                0787       _RL  swh(m,n,np+1),so2(m,n,np+1),df(m,n,np+1)
                0788       _RL  sdf(m,n),sclr(m,n),csm(m,n),x
62e9d9f553 Jean*0789 
1662f365b2 Andr*0790 c-----------------------------------------------------------------
                0791 
62e9d9f553 Jean*0792       do j= 1, n
                0793        do i= 1, m
1662f365b2 Andr*0794 
62e9d9f553 Jean*0795          swh(i,j,1)=0.
                0796          so2(i,j,1)=0.
1662f365b2 Andr*0797 
                0798 c-----csm is the effective secant of the solar zenith angle
62e9d9f553 Jean*0799 c     see equation (12) of Lacis and Hansen (1974, JAS)
                0800 
1662f365b2 Andr*0801          csm(i,j)=35./sqrt(1224.*cosz(i,j)*cosz(i,j)+1.)
                0802 
62e9d9f553 Jean*0803        enddo
1662f365b2 Andr*0804       enddo
                0805 
                0806       do k= 1, np
                0807        do j= 1, n
                0808          do i= 1, m
                0809 
62e9d9f553 Jean*0810 c-----compute layer thickness and pressure-scaling function.
1662f365b2 Andr*0811 c     indices for the surface level and surface layer
                0812 c     are np+1 and np, respectively.
62e9d9f553 Jean*0813 
1662f365b2 Andr*0814           dp(i,j,k)=pl(i,j,k+1)-pl(i,j,k)
                0815           scal(i,j,k)=dp(i,j,k)*(.5*(pl(i,j,k)+pl(i,j,k+1))/300.)**.8
62e9d9f553 Jean*0816 
1662f365b2 Andr*0817 c-----compute scaled water vapor amount, unit is g/cm**2
                0818 
                0819           wh(i,j,k)=1.02*wa(i,j,k)*scal(i,j,k)*
                0820      *              (1.+0.00135*(ta(i,j,k)-240.))
                0821           swh(i,j,k+1)=swh(i,j,k)+wh(i,j,k)
                0822 
                0823 c-----compute ozone amount, unit is (cm-atm)stp.
62e9d9f553 Jean*0824 
1662f365b2 Andr*0825           oh(i,j,k)=1.02*oa(i,j,k)*dp(i,j,k)*466.7
                0826 
                0827         enddo
                0828        enddo
                0829       enddo
                0830 
                0831 c-----scale cloud optical thickness in each layer from taucld (with
                0832 c     cloud amount fcld) to tauclb and tauclf (with cloud amount cc).
                0833 c     tauclb is the scaled optical thickness for beam radiation and
                0834 c     tauclf is for diffuse radiation.
                0835 
                0836       call cldscale(m,n,ndim,np,cosz,fcld,taucld,ict,icb,
                0837      *              cc,tauclb,tauclf)
                0838 
                0839 c-----initialize fluxes for all-sky (flx), clear-sky (flc), and
                0840 c     flux reduction (df)
                0841 
                0842       do k=1, np+1
                0843        do j=1, n
                0844         do i=1, m
                0845           flx(i,j,k)=0.
                0846           flc(i,j,k)=0.
                0847           df(i,j,k)=0.
                0848         enddo
                0849        enddo
                0850       enddo
                0851 
                0852 c-----compute solar ir fluxes
                0853 
                0854       call solir (m,n,ndim,np,wh,taucld,tauclb,tauclf,reff,ict,icb
                0855      *     ,fcld,cc,taual,csm,rsirbm,rsirdf,flx,flc,fdirir,fdifir)
                0856 
                0857 c-----compute and update uv and par fluxes
                0858 
                0859       call soluv (m,n,ndim,np,oh,dp,taucld,tauclb,tauclf,reff,ict,icb
                0860      *     ,fcld,cc,taual,csm,rsuvbm,rsuvdf,flx,flc
                0861      *     ,fdirpar,fdifpar,fdiruv,fdifuv)
                0862 
                0863 c-----compute scaled amount of o2 (so2), unit is (cm-atm)stp.
                0864 
                0865       do k= 1, np
                0866        do j= 1, n
                0867         do i= 1, m
                0868           so2(i,j,k+1)=so2(i,j,k)+165.22*scal(i,j,k)
                0869         enddo
                0870        enddo
                0871       enddo
                0872 
                0873 c-----compute flux reduction due to oxygen following
                0874 c      chou (J. climate, 1990). The fraction 0.0287 is the
                0875 c      extraterrestrial solar flux in the o2 bands.
                0876 
                0877        do k= 2, np+1
                0878         do j= 1, n
                0879          do i= 1, m
                0880            x=so2(i,j,k)*csm(i,j)
                0881            df(i,j,k)=df(i,j,k)+0.0287*(1.-expmn(-0.00027*sqrt(x)))
                0882          enddo
                0883         enddo
62e9d9f553 Jean*0884        enddo
1662f365b2 Andr*0885 
                0886 c-----compute scaled amounts for co2 (so2). unit is (cm-atm)stp.
                0887 
                0888       do k= 1, np
                0889        do j= 1, n
                0890         do i= 1, m
                0891          so2(i,j,k+1)=so2(i,j,k)+co2*789.*scal(i,j,k)
                0892         enddo
                0893        enddo
                0894       enddo
                0895 
                0896 c-----compute and update flux reduction due to co2 following
                0897 c     chou (J. Climate, 1990)
                0898 
                0899       call flxco2(m,n,np,so2,swh,csm,df)
                0900 
                0901 c-----adjust for the effect of o2 cnd co2 on clear-sky fluxes.
                0902 
                0903       do k= 2, np+1
                0904        do j= 1, n
                0905         do i= 1, m
                0906           flc(i,j,k)=flc(i,j,k)-df(i,j,k)
                0907         enddo
                0908        enddo
                0909       enddo
                0910 
                0911 c-----adjust for the all-sky fluxes due to o2 and co2.  It is
                0912 c     assumed that o2 and co2 have no effects on solar radiation
                0913 c     below clouds. clouds are assumed randomly overlapped.
                0914 
                0915       do j=1,n
                0916        do i=1,m
                0917         sdf(i,j)=0.0
                0918         sclr(i,j)=1.0
                0919        enddo
                0920       enddo
                0921 
                0922       do k=1,np
                0923        do j=1,n
                0924         do i=1,m
                0925 
                0926 c-----sclr is the fraction of clear sky.
                0927 c     sdf is the flux reduction below clouds.
                0928 
                0929          if(fcld(i,j,k).gt.0.01) then
                0930           sdf(i,j)=sdf(i,j)+df(i,j,k)*sclr(i,j)*fcld(i,j,k)
                0931           sclr(i,j)=sclr(i,j)*(1.-fcld(i,j,k))
                0932          endif
                0933           flx(i,j,k+1)=flx(i,j,k+1)-sdf(i,j)-df(i,j,k+1)*sclr(i,j)
                0934 
                0935         enddo
                0936        enddo
                0937       enddo
                0938 
                0939 c-----adjust for the direct downward ir flux.
                0940       do j= 1, n
                0941        do i= 1, m
                0942          x = (1.-rsirdf(i,j))*fdifir(i,j) + (1.-rsirbm(i,j))*fdirir(i,j)
                0943          x = (-sdf(i,j)-df(i,j,np+1)*sclr(i,j))/x
                0944          fdirir(i,j)=fdirir(i,j)*(1.+x)
                0945          fdifir(i,j)=fdifir(i,j)*(1.+x)
                0946        enddo
                0947       enddo
                0948 
                0949       return
62e9d9f553 Jean*0950       end
1662f365b2 Andr*0951 
                0952 c********************************************************************
                0953 
                0954       subroutine cldscale (m,n,ndim,np,cosz,fcld,taucld,ict,icb,
                0955      *                     cc,tauclb,tauclf)
                0956 
                0957 c********************************************************************
                0958 c
                0959 c   This subroutine computes the covers of high, middle, and
                0960 c    low clouds and scales the cloud optical thickness.
                0961 c
                0962 c   To simplify calculations in a cloudy atmosphere, clouds are
                0963 c    grouped into high, middle and low clouds separated by the levels
                0964 c    ict and icb (level 1 is the top of the atmosphere).
                0965 c
                0966 c   Within each of the three groups, clouds are assumed maximally
                0967 c    overlapped, and the cloud cover (cc) of a group is the maximum
                0968 c    cloud cover of all the layers in the group.  The optical thickness
                0969 c    (taucld) of a given layer is then scaled to new values (tauclb and
                0970 c    tauclf) so that the layer reflectance corresponding to the cloud
                0971 c    cover cc is the same as the original reflectance with optical
                0972 c    thickness taucld and cloud cover fcld.
                0973 c
                0974 c---input parameters
                0975 c
                0976 c    number of grid intervals in zonal direction (m)
                0977 c    number of grid intervals in meridional direction (n)
                0978 c    maximum number of grid intervals in meridional direction (ndim)
                0979 c    number of atmospheric layers (np)
                0980 c    cosine of the solar zenith angle (cosz)
                0981 c    fractional cloud cover (fcld)
                0982 c    cloud optical thickness (taucld)
                0983 c    index separating high and middle clouds (ict)
                0984 c    index separating middle and low clouds (icb)
                0985 c
                0986 c---output parameters
                0987 c
                0988 c    fractional cover of high, middle, and low clouds (cc)
                0989 c    scaled cloud optical thickness for beam radiation (tauclb)
                0990 c    scaled cloud optical thickness for diffuse radiation (tauclf)
                0991 c
                0992 c********************************************************************
                0993 
                0994       implicit none
                0995 
                0996 c-----input parameters
                0997 
                0998       integer m,n,ndim,np,ict,icb
a456aa407c Andr*0999       _RL  cosz(m,ndim),fcld(m,ndim,np),taucld(m,ndim,np,2)
1662f365b2 Andr*1000 
                1001 c-----output parameters
                1002 
a456aa407c Andr*1003       _RL  cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np)
1662f365b2 Andr*1004 
                1005 c-----temporary variables
                1006 
                1007       integer i,j,k,im,it,ia,kk
a456aa407c Andr*1008       _RL   fm,ft,fa,xai,taux
1662f365b2 Andr*1009 
                1010 c-----pre-computed table
                1011 
                1012       integer   nm,nt,na
62e9d9f553 Jean*1013       parameter (nm=11,nt=9,na=11)
a456aa407c Andr*1014       _RL   dm,dt,da,t1,caib(nm,nt,na),caif(nt,na)
1662f365b2 Andr*1015       parameter (dm=0.1,dt=0.30103,da=0.1,t1=-0.9031)
                1016 
                1017 c-----include the pre-computed table for cai
                1018 
4a402f223d Andr*1019 #include "cai-dat.h"
3f93544ef8 Andr*1020       save caib,caif
1662f365b2 Andr*1021 
                1022 c-----clouds within each of the high, middle, and low clouds are
                1023 c     assumed maximally overlapped, and the cloud cover (cc)
                1024 c     for a group is the maximum cloud cover of all the layers
                1025 c     in the group
                1026 
                1027       do j=1,n
                1028        do i=1,m
                1029           cc(i,j,1)=0.0
                1030           cc(i,j,2)=0.0
                1031           cc(i,j,3)=0.0
                1032        enddo
                1033       enddo
                1034 
                1035       do k=1,ict-1
                1036        do j=1,n
                1037         do i=1,m
                1038            cc(i,j,1)=max(cc(i,j,1),fcld(i,j,k))
                1039         enddo
                1040        enddo
                1041       enddo
                1042 
                1043       do k=ict,icb-1
                1044        do j=1,n
                1045         do i=1,m
                1046            cc(i,j,2)=max(cc(i,j,2),fcld(i,j,k))
                1047         enddo
                1048        enddo
                1049       enddo
                1050 
                1051       do k=icb,np
                1052        do j=1,n
                1053         do i=1,m
                1054            cc(i,j,3)=max(cc(i,j,3),fcld(i,j,k))
                1055         enddo
                1056        enddo
                1057       enddo
                1058 
                1059 c-----scale the cloud optical thickness.
                1060 c     taucld(i,j,k,1) is the optical thickness for ice particles, and
                1061 c     taucld(i,j,k,2) is the optical thickness for liquid particles.
62e9d9f553 Jean*1062 
1662f365b2 Andr*1063       do k=1,np
                1064 
                1065          if(k.lt.ict) then
                1066             kk=1
                1067          elseif(k.ge.ict .and. k.lt.icb) then
                1068             kk=2
                1069          else
                1070             kk=3
                1071          endif
                1072 
                1073        do j=1,n
                1074         do i=1,m
                1075 
                1076          tauclb(i,j,k) = 0.0
                1077          tauclf(i,j,k) = 0.0
                1078 
                1079          taux=taucld(i,j,k,1)+taucld(i,j,k,2)
                1080          if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then
                1081 
                1082 c-----normalize cloud cover
                1083 
                1084            fa=fcld(i,j,k)/cc(i,j,kk)
                1085 
                1086 c-----table look-up
                1087 
32362b8fa8 Cons*1088            taux=min(taux,32. _d 0)
1662f365b2 Andr*1089 
                1090            fm=cosz(i,j)/dm
                1091            ft=(log10(taux)-t1)/dt
                1092            fa=fa/da
62e9d9f553 Jean*1093 
1662f365b2 Andr*1094            im=int(fm+1.5)
                1095            it=int(ft+1.5)
                1096            ia=int(fa+1.5)
62e9d9f553 Jean*1097 
1662f365b2 Andr*1098            im=max(im,2)
                1099            it=max(it,2)
                1100            ia=max(ia,2)
62e9d9f553 Jean*1101 
1662f365b2 Andr*1102            im=min(im,nm-1)
                1103            it=min(it,nt-1)
                1104            ia=min(ia,na-1)
                1105 
                1106            fm=fm-float(im-1)
                1107            ft=ft-float(it-1)
                1108            fa=fa-float(ia-1)
                1109 
                1110 c-----scale cloud optical thickness for beam radiation.
                1111 c     the scaling factor, xai, is a function of the solar zenith
                1112 c     angle, optical thickness, and cloud cover.
62e9d9f553 Jean*1113 
1662f365b2 Andr*1114            xai=    (-caib(im-1,it,ia)*(1.-fm)+
                1115      *      caib(im+1,it,ia)*(1.+fm))*fm*.5+caib(im,it,ia)*(1.-fm*fm)
62e9d9f553 Jean*1116 
1662f365b2 Andr*1117            xai=xai+(-caib(im,it-1,ia)*(1.-ft)+
                1118      *      caib(im,it+1,ia)*(1.+ft))*ft*.5+caib(im,it,ia)*(1.-ft*ft)
                1119 
                1120            xai=xai+(-caib(im,it,ia-1)*(1.-fa)+
                1121      *     caib(im,it,ia+1)*(1.+fa))*fa*.5+caib(im,it,ia)*(1.-fa*fa)
                1122 
                1123            xai= xai-2.*caib(im,it,ia)
32362b8fa8 Cons*1124            xai=max(xai,0.0 _d 0)
62e9d9f553 Jean*1125 
1662f365b2 Andr*1126            tauclb(i,j,k) = taux*xai
                1127 
                1128 c-----scale cloud optical thickness for diffuse radiation.
                1129 c     the scaling factor, xai, is a function of the cloud optical
                1130 c     thickness and cover but not the solar zenith angle.
                1131 
                1132            xai=    (-caif(it-1,ia)*(1.-ft)+
                1133      *      caif(it+1,ia)*(1.+ft))*ft*.5+caif(it,ia)*(1.-ft*ft)
                1134 
                1135            xai=xai+(-caif(it,ia-1)*(1.-fa)+
                1136      *      caif(it,ia+1)*(1.+fa))*fa*.5+caif(it,ia)*(1.-fa*fa)
                1137 
                1138            xai= xai-caif(it,ia)
32362b8fa8 Cons*1139            xai=max(xai,0.0 _d 0)
62e9d9f553 Jean*1140 
1662f365b2 Andr*1141            tauclf(i,j,k) = taux*xai
                1142 
                1143          endif
                1144 
                1145         enddo
                1146        enddo
                1147       enddo
                1148 
                1149       return
                1150       end
                1151 c***********************************************************************
                1152 
                1153       subroutine solir (m,n,ndim,np,wh,taucld,tauclb,tauclf,reff,
                1154      *  ict,icb,fcld,cc,taual,csm,rsirbm,rsirdf,flx,flc,fdirir,fdifir)
62e9d9f553 Jean*1155 
1662f365b2 Andr*1156 c************************************************************************
                1157 c  compute solar flux in the infrared region. The spectrum is divided
                1158 c   into three bands:
                1159 c
                1160 c          band   wavenumber(/cm)  wavelength (micron)
                1161 c           1       1000-4400         2.27-10.0
                1162 c           2       4400-8200         1.22-2.27
                1163 c           3       8200-14300        0.70-1.22
                1164 c
                1165 c----- Input parameters:                            units      size
                1166 c
                1167 c   number of soundings in zonal direction (m)       n/d        1
                1168 c   number of soundings in meridional direction (n)  n/d        1
                1169 c   maximum number of soundings in                   n/d        1
                1170 c          meridional direction (ndim)
                1171 c   number of atmospheric layers (np)                n/d        1
                1172 c   layer water vapor content (wh)                 gm/cm^2    m*n*np
                1173 c   cloud optical thickness (taucld)                 n/d      m*ndim*np*2
                1174 c          index 1 for ice paticles
                1175 c          index 2 for liquid particles
                1176 c   scaled cloud optical thickness                   n/d      m*n*np
                1177 c          for beam radiation (tauclb)
                1178 c   scaled cloud optical thickness                   n/d      m*n*np
                1179 c          for diffuse radiation  (tauclf)
                1180 c   effective cloud-particle size (reff)           micrometer m*ndim*np*2
                1181 c          index 1 for ice paticles
                1182 c          index 2 for liquid particles
                1183 c   level index separating high and                  n/d      m*n
                1184 c          middle clouds (ict)
                1185 c   level index separating middle and                n/d      m*n
                1186 c          low clouds (icb)
                1187 c   cloud amount (fcld)                            fraction   m*ndim*np
                1188 c   cloud amount of high, middle, and                n/d      m*n*3
                1189 c          low clouds (cc)
                1190 c   aerosol optical thickness (taual)                n/d      m*ndim*np
                1191 c   cosecant of the solar zenith angle (csm)         n/d      m*n
                1192 c   near ir surface albedo for beam                fraction   m*ndim
                1193 c                radiation (rsirbm)
                1194 c   near ir surface albedo for diffuse             fraction   m*ndim
                1195 c                radiation (rsirdf)
                1196 c
                1197 c----- output (updated) parameters:
                1198 c
                1199 c   all-sky flux (downward-upward) (flx)           fraction   m*ndim*(np+1)
                1200 c   clear-sky flux (downward-upward) (flc)         fraction   m*ndim*(np+1)
                1201 c   all-sky direct downward ir flux at
                1202 c          the surface (fdirir)                    fraction   m*ndim
                1203 c   all-sky diffuse downward ir flux at
                1204 c          the surface (fdifir)                    fraction   m*ndim
                1205 c
                1206 c----- note: the following parameters must be specified by users:
                1207 c   aerosol single scattering albedo (ssaal)         n/d      nband
                1208 c   aerosol asymmetry factor (asyal)                 n/d      nband
                1209 c
                1210 c*************************************************************************
                1211 
                1212       implicit none
                1213 
62e9d9f553 Jean*1214 c-----Explicit Inline Directives
                1215 
06d4643e1f Jean*1216 #ifdef FIZHI_CRAY
4e1c053948 Jean*1217 #ifdef FIZHI_F77_COMPIL
1662f365b2 Andr*1218 cfpp$ expand (deledd)
                1219 cfpp$ expand (sagpol)
62e9d9f553 Jean*1220 cfpp$ expand (expmn)
                1221 #endif
                1222 #endif
a456aa407c Andr*1223       _RL expmn
1662f365b2 Andr*1224 
                1225 c-----input parameters
                1226 
                1227       integer m,n,ndim,np,ict,icb
a456aa407c Andr*1228       _RL  taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np)
                1229       _RL  tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3)
                1230       _RL  rsirbm(m,ndim),rsirdf(m,ndim)
                1231       _RL  wh(m,n,np),taual(m,ndim,np),csm(m,n)
1662f365b2 Andr*1232 
                1233 c-----output (updated) parameters
                1234 
a456aa407c Andr*1235       _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
                1236       _RL  fdirir(m,ndim),fdifir(m,ndim)
1662f365b2 Andr*1237 
                1238 c-----static parameters
                1239 
                1240       integer nk,nband
                1241       parameter (nk=10,nband=3)
a456aa407c Andr*1242       _RL  xk(nk),hk(nband,nk),ssaal(nband),asyal(nband)
                1243       _RL  aia(nband,3),awa(nband,3),aig(nband,3),awg(nband,3)
1662f365b2 Andr*1244 
                1245 c-----temporary array
                1246 
                1247       integer ib,ik,i,j,k
a456aa407c Andr*1248       _RL  ssacl(m,n,np),asycl(m,n,np)
                1249       _RL  rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2),
1662f365b2 Andr*1250      *       rs(m,n,np+1,2),ts(m,n,np+1,2)
a456aa407c Andr*1251       _RL  fall(m,n,np+1),fclr(m,n,np+1)
                1252       _RL  fsdir(m,n),fsdif(m,n)
1662f365b2 Andr*1253 
a456aa407c Andr*1254       _RL  tauwv,tausto,ssatau,asysto,tauto,ssato,asyto
                1255       _RL  taux,reff1,reff2,w1,w2,g1,g2
                1256       _RL  ssaclt(m,n),asyclt(m,n)
                1257       _RL  rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)
                1258       _RL  rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)
1662f365b2 Andr*1259 
                1260 c-----water vapor absorption coefficient for 10 k-intervals.
                1261 c     unit: cm^2/gm
                1262 
62e9d9f553 Jean*1263       data xk/
                1264      1  0.0010, 0.0133, 0.0422, 0.1334, 0.4217,
                1265      2  1.334,  5.623,  31.62,  177.8,  1000.0/
1662f365b2 Andr*1266 
                1267 c-----water vapor k-distribution function,
                1268 c     the sum of hk is 0.52926. unit: fraction
                1269 
                1270       data hk/
                1271      1 .01074,.08236,.20673,  .00360,.01157,.03497,
                1272      2 .00411,.01133,.03011,  .00421,.01143,.02260,
                1273      3 .00389,.01240,.01336,  .00326,.01258,.00696,
                1274      4 .00499,.01381,.00441,  .00465,.00650,.00115,
                1275      5 .00245,.00244,.00026,  .00145,.00094,.00000/
                1276 
                1277 c-----aerosol single-scattering albedo and asymmetry factor
                1278 
                1279       data ssaal,asyal/nband*0.999,nband*0.850/
62e9d9f553 Jean*1280 
1662f365b2 Andr*1281 c-----coefficients for computing the single scattering albedo of
                1282 c     ice clouds from ssa=1-(aia(*,1)+aia(*,2)*reff+aia(*,3)*reff**2)
                1283 
                1284       data aia/
                1285      1  .08938331, .00215346,-.00000260,
                1286      2  .00299387, .00073709, .00000746,
                1287      3 -.00001038,-.00000134, .00000000/
                1288 
                1289 c-----coefficients for computing the single scattering albedo of
                1290 c     liquid clouds from ssa=1-(awa(*,1)+awa(*,2)*reff+awa(*,3)*reff**2)
                1291 
                1292       data awa/
                1293      1  .01209318,-.00019934, .00000007,
                1294      2  .01784739, .00088757, .00000845,
                1295      3 -.00036910,-.00000650,-.00000004/
                1296 
                1297 c-----coefficients for computing the asymmetry factor of ice clouds
                1298 c     from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2
                1299 
                1300       data aig/
                1301      1  .84090400, .76098937, .74935228,
                1302      2  .00126222, .00141864, .00119715,
                1303      3 -.00000385,-.00000396,-.00000367/
                1304 
                1305 c-----coefficients for computing the asymmetry factor of liquid clouds
                1306 c     from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2
                1307 
                1308       data awg/
                1309      1  .83530748, .74513197, .79375035,
                1310      2  .00257181, .01370071, .00832441,
                1311      3  .00005519,-.00038203,-.00023263/
                1312 
                1313 c-----initialize surface fluxes, reflectances, and transmittances
                1314 
                1315       do j= 1, n
                1316        do i= 1, m
                1317          fdirir(i,j)=0.0
                1318          fdifir(i,j)=0.0
                1319          rr(i,j,np+1,1)=rsirbm(i,j)
                1320          rr(i,j,np+1,2)=rsirbm(i,j)
                1321          rs(i,j,np+1,1)=rsirdf(i,j)
                1322          rs(i,j,np+1,2)=rsirdf(i,j)
                1323          td(i,j,np+1,1)=0.0
                1324          td(i,j,np+1,2)=0.0
                1325          tt(i,j,np+1,1)=0.0
                1326          tt(i,j,np+1,2)=0.0
                1327          ts(i,j,np+1,1)=0.0
                1328          ts(i,j,np+1,2)=0.0
                1329        enddo
                1330       enddo
                1331 
                1332 c-----integration over spectral bands
                1333 
                1334       do 100 ib=1,nband
                1335 
                1336 c-----compute cloud single scattering albedo and asymmetry factor
                1337 c     for a mixture of ice and liquid particles.
                1338 
                1339        do k= 1, np
                1340 
                1341         do j= 1, n
                1342          do i= 1, m
                1343 
                1344            ssaclt(i,j)=1.0
                1345            asyclt(i,j)=1.0
                1346 
                1347            taux=taucld(i,j,k,1)+taucld(i,j,k,2)
                1348           if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then
                1349 
32362b8fa8 Cons*1350            reff1=min(reff(i,j,k,1),130. _d 0)
                1351            reff2=min(reff(i,j,k,2),20.0 _d 0)
1662f365b2 Andr*1352 
                1353            w1=(1.-(aia(ib,1)+(aia(ib,2)+
                1354      *         aia(ib,3)*reff1)*reff1))*taucld(i,j,k,1)
                1355            w2=(1.-(awa(ib,1)+(awa(ib,2)+
                1356      *         awa(ib,3)*reff2)*reff2))*taucld(i,j,k,2)
                1357            ssaclt(i,j)=(w1+w2)/taux
                1358 
                1359            g1=(aig(ib,1)+(aig(ib,2)+aig(ib,3)*reff1)*reff1)*w1
                1360            g2=(awg(ib,1)+(awg(ib,2)+awg(ib,3)*reff2)*reff2)*w2
                1361            asyclt(i,j)=(g1+g2)/(w1+w2)
                1362 
                1363           endif
                1364 
                1365          enddo
                1366         enddo
                1367 
                1368         do j=1,n
                1369          do i=1,m
                1370             ssacl(i,j,k)=ssaclt(i,j)
                1371          enddo
                1372         enddo
                1373         do j=1,n
                1374          do i=1,m
                1375             asycl(i,j,k)=asyclt(i,j)
                1376          enddo
                1377         enddo
                1378 
                1379        enddo
                1380 
                1381 c-----integration over the k-distribution function
                1382 
                1383          do 200 ik=1,nk
                1384 
                1385           do 300 k= 1, np
                1386 
                1387            do j= 1, n
                1388             do i= 1, m
                1389 
                1390              tauwv=xk(ik)*wh(i,j,k)
9524fe64b5 Andr*1391 
1662f365b2 Andr*1392 c-----compute total optical thickness, single scattering albedo,
                1393 c     and asymmetry factor.
62e9d9f553 Jean*1394 
1662f365b2 Andr*1395              tausto=tauwv+taual(i,j,k)+1.0e-8
                1396              ssatau=ssaal(ib)*taual(i,j,k)
                1397              asysto=asyal(ib)*ssaal(ib)*taual(i,j,k)
62e9d9f553 Jean*1398 
1662f365b2 Andr*1399 c-----compute reflectance and transmittance for cloudless layers
                1400 
                1401              tauto=tausto
                1402              ssato=ssatau/tauto+1.0e-8
                1403 
                1404 c            if (ssato .gt. 0.001) then
                1405 
32362b8fa8 Cons*1406 c             ssato=min(ssato,0.999999 _d 0)
1662f365b2 Andr*1407 c             asyto=asysto/(ssato*tauto)
                1408 
62e9d9f553 Jean*1409 c             call deledd(tauto,ssato,asyto,csm(i,j),
1662f365b2 Andr*1410 c     *                   rr1t(i,j),tt1t(i,j),td1t(i,j))
                1411 
                1412 c             call sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j))
                1413 
                1414 c            else
                1415 
                1416              td1t(i,j)=expmn(-tauto*csm(i,j))
                1417              ts1t(i,j)=expmn(-1.66*tauto)
                1418              tt1t(i,j)=0.0
                1419              rr1t(i,j)=0.0
                1420              rs1t(i,j)=0.0
                1421 
                1422 c            endif
                1423 
                1424 c-----compute reflectance and transmittance for cloud layers
                1425 
                1426              if (tauclb(i,j,k) .lt. 0.01) then
                1427 
                1428               rr2t(i,j)=rr1t(i,j)
                1429               tt2t(i,j)=tt1t(i,j)
                1430               td2t(i,j)=td1t(i,j)
                1431               rs2t(i,j)=rs1t(i,j)
                1432               ts2t(i,j)=ts1t(i,j)
                1433 
                1434              else
                1435 
                1436               tauto=tausto+tauclb(i,j,k)
                1437               ssato=(ssatau+ssacl(i,j,k)*tauclb(i,j,k))/tauto+1.0e-8
32362b8fa8 Cons*1438               ssato=min(ssato,0.999999 _d 0)
1662f365b2 Andr*1439               asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclb(i,j,k))/
                1440      *              (ssato*tauto)
                1441 
62e9d9f553 Jean*1442               call deledd(tauto,ssato,asyto,csm(i,j),
1662f365b2 Andr*1443      *                    rr2t(i,j),tt2t(i,j),td2t(i,j))
                1444 
                1445               tauto=tausto+tauclf(i,j,k)
                1446               ssato=(ssatau+ssacl(i,j,k)*tauclf(i,j,k))/tauto+1.0e-8
32362b8fa8 Cons*1447               ssato=min(ssato,0.999999 _d 0)
1662f365b2 Andr*1448               asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclf(i,j,k))/
                1449      *              (ssato*tauto)
                1450 
                1451               call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j))
                1452 
                1453              endif
                1454 
                1455             enddo
                1456            enddo
                1457 
                1458 c-----the following code appears to be awkward, but it is efficient
                1459 c     in certain parallel processors.
                1460 
                1461            do j=1,n
                1462             do i=1,m
                1463                rr(i,j,k,1)=rr1t(i,j)
                1464             enddo
                1465            enddo
                1466            do j=1,n
                1467             do i=1,m
                1468                tt(i,j,k,1)=tt1t(i,j)
                1469             enddo
                1470            enddo
                1471            do j=1,n
                1472             do i=1,m
                1473                td(i,j,k,1)=td1t(i,j)
                1474             enddo
                1475            enddo
                1476            do j=1,n
                1477             do i=1,m
                1478                rs(i,j,k,1)=rs1t(i,j)
                1479             enddo
                1480            enddo
                1481            do j=1,n
                1482             do i=1,m
                1483                ts(i,j,k,1)=ts1t(i,j)
                1484             enddo
                1485            enddo
62e9d9f553 Jean*1486 
1662f365b2 Andr*1487            do j=1,n
                1488             do i=1,m
                1489                rr(i,j,k,2)=rr2t(i,j)
                1490             enddo
                1491            enddo
                1492            do j=1,n
                1493             do i=1,m
                1494                tt(i,j,k,2)=tt2t(i,j)
                1495             enddo
                1496            enddo
                1497            do j=1,n
                1498             do i=1,m
                1499                td(i,j,k,2)=td2t(i,j)
                1500             enddo
                1501            enddo
                1502            do j=1,n
                1503             do i=1,m
                1504                rs(i,j,k,2)=rs2t(i,j)
                1505             enddo
                1506            enddo
                1507            do j=1,n
                1508             do i=1,m
                1509                ts(i,j,k,2)=ts2t(i,j)
                1510             enddo
                1511            enddo
                1512 
                1513  300  continue
                1514 
                1515 c-----flux calculations
62e9d9f553 Jean*1516 
9524fe64b5 Andr*1517        do k= 1, np+1
                1518         do j= 1, n
                1519          do i= 1, m
                1520           fclr(i,j,k) = 0.
                1521           fall(i,j,k) = 0.
                1522          enddo
                1523         enddo
                1524        enddo
                1525        do j= 1, n
                1526         do i= 1, m
                1527          fsdir(i,j) = 0.
                1528          fsdif(i,j) = 0.
                1529         enddo
                1530        enddo
                1531 
082e38725b Andr*1532         call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,
                1533      *               fclr,fall,fsdir,fsdif)
1662f365b2 Andr*1534 
                1535        do k= 1, np+1
                1536         do j= 1, n
                1537          do i= 1, m
                1538           flx(i,j,k) = flx(i,j,k)+fall(i,j,k)*hk(ib,ik)
                1539          enddo
                1540         enddo
                1541         do j= 1, n
                1542          do i= 1, m
                1543           flc(i,j,k) = flc(i,j,k)+fclr(i,j,k)*hk(ib,ik)
                1544          enddo
                1545         enddo
                1546        enddo
                1547 
                1548        do j= 1, n
                1549         do i= 1, m
                1550           fdirir(i,j) = fdirir(i,j)+fsdir(i,j)*hk(ib,ik)
                1551           fdifir(i,j) = fdifir(i,j)+fsdif(i,j)*hk(ib,ik)
                1552         enddo
                1553        enddo
                1554 
                1555   200 continue
                1556   100 continue
62e9d9f553 Jean*1557 
1662f365b2 Andr*1558       return
                1559       end
                1560 
                1561 c************************************************************************
                1562 
                1563       subroutine soluv (m,n,ndim,np,oh,dp,taucld,tauclb,tauclf,reff,
                1564      *  ict,icb,fcld,cc,taual,csm,rsuvbm,rsuvdf,flx,flc
                1565      * ,fdirpar,fdifpar,fdiruv,fdifuv)
                1566 
                1567 c************************************************************************
                1568 c  compute solar fluxes in the uv+visible region. the spectrum is
                1569 c  grouped into 8 bands:
62e9d9f553 Jean*1570 c
1662f365b2 Andr*1571 c              Band     Micrometer
                1572 c
                1573 c       UV-C    1.     .175 - .225
                1574 c               2.     .225 - .245
                1575 c                      .260 - .280
                1576 c               3.     .245 - .260
                1577 c
                1578 c       UV-B    4.     .280 - .295
                1579 c               5.     .295 - .310
                1580 c               6.     .310 - .320
62e9d9f553 Jean*1581 c
1662f365b2 Andr*1582 c       UV-A    7.     .320 - .400
62e9d9f553 Jean*1583 c
1662f365b2 Andr*1584 c       PAR     8.     .400 - .700
                1585 c
                1586 c----- Input parameters:                            units      size
                1587 c
                1588 c   number of soundings in zonal direction (m)       n/d        1
                1589 c   number of soundings in meridional direction (n)  n/d        1
                1590 c   maximum number of soundings in                   n/d        1
                1591 c           meridional direction (ndim)
                1592 c   number of atmospheric layers (np)                n/d        1
                1593 c   layer ozone content (oh)                      (cm-atm)stp m*n*np
                1594 c   layer pressure thickness (dp)                    mb       m*n*np
                1595 c   cloud optical thickness (taucld)                 n/d      m*ndim*np*2
                1596 c          index 1 for ice paticles
                1597 c          index 2 for liquid particles
                1598 c   scaled cloud optical thickness                   n/d      m*n*np
                1599 c          for beam radiation (tauclb)
                1600 c   scaled cloud optical thickness                   n/d      m*n*np
                1601 c          for diffuse radiation  (tauclf)
                1602 c   effective cloud-particle size (reff)           micrometer m*ndim*np*2
                1603 c          index 1 for ice paticles
                1604 c          index 2 for liquid particles
                1605 c   level indiex separating high and                 n/d      m*n
                1606 c          middle clouds (ict)
                1607 c   level indiex separating middle and               n/d      m*n
                1608 c          low clouds (icb)
                1609 c   cloud amount (fcld)                            fraction   m*ndim*np
                1610 c   cloud amounts of high, middle, and               n/d      m*n*3
                1611 c          low clouds (cc)
                1612 c   aerosol optical thickness (taual)                n/d      m*ndim*np
                1613 c   cosecant of the solar zenith angle (csm)         n/d      m*n
                1614 c   uv+par surface albedo for beam                 fraction   m*ndim
                1615 c           radiation (rsuvbm)
                1616 c   uv+par surface albedo for diffuse              fraction   m*ndim
                1617 c           radiation (rsuvdf)
                1618 c
                1619 c----- output (updated) parameters:
                1620 c
                1621 c   all-sky net downward flux (flx)                fraction   m*ndim*(np+1)
                1622 c   clear-sky net downward flux (flc)              fraction   m*ndim*(np+1)
                1623 c   all-sky direct downward par flux at
                1624 c          the surface (fdirpar)                   fraction   m*ndim
                1625 c   all-sky diffuse downward par flux at
                1626 c          the surface (fdifpar)                   fraction   m*ndim
                1627 c
                1628 c----- note: the following parameters must be specified by users:
                1629 c
                1630 c   aerosol single scattering albedo (ssaal)         n/d        1
                1631 c   aerosol asymmetry factor (asyal)                 n/d        1
                1632 c
                1633 *
                1634 c***********************************************************************
                1635 
                1636       implicit none
                1637 
62e9d9f553 Jean*1638 c-----Explicit Inline Directives
                1639 
06d4643e1f Jean*1640 #ifdef FIZHI_CRAY
4e1c053948 Jean*1641 #ifdef FIZHI_F77_COMPIL
62e9d9f553 Jean*1642 cfpp$ expand (deledd)
                1643 cfpp$ expand (sagpol)
                1644 #endif
1662f365b2 Andr*1645 #endif
                1646 
                1647 c-----input parameters
                1648 
                1649       integer m,n,ndim,np,ict,icb
a456aa407c Andr*1650       _RL  taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np)
                1651       _RL  tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3)
                1652       _RL  oh(m,n,np),dp(m,n,np),taual(m,ndim,np)
                1653       _RL  rsuvbm(m,ndim),rsuvdf(m,ndim),csm(m,n)
1662f365b2 Andr*1654 
                1655 c-----output (updated) parameter
                1656 
a456aa407c Andr*1657       _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
                1658       _RL  fdirpar(m,ndim),fdifpar(m,ndim)
                1659       _RL  fdiruv(m,ndim),fdifuv(m,ndim)
1662f365b2 Andr*1660 
                1661 c-----static parameters
                1662 
                1663       integer nband
                1664       parameter (nband=8)
a456aa407c Andr*1665       _RL  hk(nband),xk(nband),ry(nband)
                1666       _RL  asyal(nband),ssaal(nband),aig(3),awg(3)
1662f365b2 Andr*1667 
                1668 c-----temporary array
                1669 
                1670       integer i,j,k,ib
a456aa407c Andr*1671       _RL  taurs,tauoz,tausto,ssatau,asysto,tauto,ssato,asyto
                1672       _RL  taux,reff1,reff2,g1,g2,asycl(m,n,np)
                1673       _RL  td(m,n,np+1,2),rr(m,n,np+1,2),tt(m,n,np+1,2),
1662f365b2 Andr*1674      *       rs(m,n,np+1,2),ts(m,n,np+1,2)
a456aa407c Andr*1675       _RL  fall(m,n,np+1),fclr(m,n,np+1),fsdir(m,n),fsdif(m,n)
                1676       _RL  asyclt(m,n)
                1677       _RL  rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)
                1678       _RL  rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)
1662f365b2 Andr*1679 
                1680 c-----hk is the fractional extra-terrestrial solar flux.
                1681 c     the sum of hk is 0.47074.
                1682 
                1683       data hk/.00057, .00367, .00083, .00417,
                1684      *        .00600, .00556, .05913, .39081/
                1685 
                1686 c-----xk is the ozone absorption coefficient. unit: /(cm-atm)stp
                1687 
                1688       data xk /30.47, 187.2,  301.9,   42.83,
                1689      *         7.09,  1.25,   0.0345,  0.0539/
                1690 
                1691 c-----ry is the extinction coefficient for Rayleigh scattering.
                1692 c     unit: /mb.
                1693 
                1694       data ry /.00604, .00170, .00222, .00132,
                1695      *         .00107, .00091, .00055, .00012/
                1696 
                1697 c-----aerosol single-scattering albedo and asymmetry factor
                1698 
                1699       data ssaal,asyal/nband*0.999,nband*0.850/
                1700 
                1701 c-----coefficients for computing the asymmetry factor of ice clouds
                1702 c     from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2
                1703 
                1704       data aig/.74625000,.00105410,-.00000264/
                1705 
                1706 c-----coefficients for computing the asymmetry factor of liquid
                1707 c     clouds from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2
                1708 
                1709       data awg/.82562000,.00529000,-.00014866/
                1710 
                1711 c-----initialize surface reflectances and transmittances
62e9d9f553 Jean*1712 
1662f365b2 Andr*1713       do j= 1, n
62e9d9f553 Jean*1714        do i= 1, m
1662f365b2 Andr*1715          rr(i,j,np+1,1)=rsuvbm(i,j)
                1716          rr(i,j,np+1,2)=rsuvbm(i,j)
                1717          rs(i,j,np+1,1)=rsuvdf(i,j)
                1718          rs(i,j,np+1,2)=rsuvdf(i,j)
                1719          td(i,j,np+1,1)=0.0
                1720          td(i,j,np+1,2)=0.0
                1721          tt(i,j,np+1,1)=0.0
                1722          tt(i,j,np+1,2)=0.0
                1723          ts(i,j,np+1,1)=0.0
                1724          ts(i,j,np+1,2)=0.0
                1725        enddo
                1726       enddo
                1727 
                1728 c-----compute cloud asymmetry factor for a mixture of
                1729 c     liquid and ice particles.  unit of reff is micrometers.
                1730 
                1731       do k= 1, np
                1732 
                1733        do j= 1, n
                1734         do i= 1, m
                1735 
                1736            asyclt(i,j)=1.0
                1737 
                1738            taux=taucld(i,j,k,1)+taucld(i,j,k,2)
                1739           if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then
                1740 
32362b8fa8 Cons*1741            reff1=min(reff(i,j,k,1),130. _d 0)
                1742            reff2=min(reff(i,j,k,2),20.0 _d 0)
1662f365b2 Andr*1743 
                1744            g1=(aig(1)+(aig(2)+aig(3)*reff1)*reff1)*taucld(i,j,k,1)
                1745            g2=(awg(1)+(awg(2)+awg(3)*reff2)*reff2)*taucld(i,j,k,2)
                1746            asyclt(i,j)=(g1+g2)/taux
                1747 
                1748           endif
                1749 
                1750         enddo
                1751        enddo
                1752 
                1753        do j=1,n
                1754         do i=1,m
                1755            asycl(i,j,k)=asyclt(i,j)
                1756         enddo
                1757        enddo
                1758 
                1759       enddo
62e9d9f553 Jean*1760 
1662f365b2 Andr*1761       do j=1,n
                1762        do i=1,m
                1763         fdiruv(i,j)=0.0
                1764         fdifuv(i,j)=0.0
                1765        enddo
                1766       enddo
62e9d9f553 Jean*1767 
1662f365b2 Andr*1768 c-----integration over spectral bands
                1769 
                1770       do 100 ib=1,nband
                1771 
                1772        do 300 k= 1, np
                1773 
                1774         do j= 1, n
                1775          do i= 1, m
                1776 
                1777 c-----compute ozone and rayleigh optical thicknesses
                1778 
                1779           taurs=ry(ib)*dp(i,j,k)
                1780           tauoz=xk(ib)*oh(i,j,k)
62e9d9f553 Jean*1781 
1662f365b2 Andr*1782 c-----compute total optical thickness, single scattering albedo,
                1783 c     and asymmetry factor
                1784 
                1785           tausto=taurs+tauoz+taual(i,j,k)+1.0e-8
                1786           ssatau=ssaal(ib)*taual(i,j,k)+taurs
                1787           asysto=asyal(ib)*ssaal(ib)*taual(i,j,k)
                1788 
                1789 c-----compute reflectance and transmittance for cloudless layers
                1790 
                1791           tauto=tausto
                1792           ssato=ssatau/tauto+1.0e-8
32362b8fa8 Cons*1793           ssato=min(ssato,0.999999 _d 0)
1662f365b2 Andr*1794           asyto=asysto/(ssato*tauto)
                1795 
62e9d9f553 Jean*1796           call deledd(tauto,ssato,asyto,csm(i,j),
1662f365b2 Andr*1797      *                rr1t(i,j),tt1t(i,j),td1t(i,j))
                1798 
                1799           call sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j))
                1800 
                1801 c-----compute reflectance and transmittance for cloud layers
                1802 
                1803           if (tauclb(i,j,k) .lt. 0.01) then
                1804 
                1805            rr2t(i,j)=rr1t(i,j)
                1806            tt2t(i,j)=tt1t(i,j)
                1807            td2t(i,j)=td1t(i,j)
                1808            rs2t(i,j)=rs1t(i,j)
                1809            ts2t(i,j)=ts1t(i,j)
                1810 
                1811           else
                1812 
                1813            tauto=tausto+tauclb(i,j,k)
                1814            ssato=(ssatau+tauclb(i,j,k))/tauto+1.0e-8
32362b8fa8 Cons*1815            ssato=min(ssato,0.999999 _d 0)
1662f365b2 Andr*1816            asyto=(asysto+asycl(i,j,k)*tauclb(i,j,k))/(ssato*tauto)
                1817 
62e9d9f553 Jean*1818            call deledd(tauto,ssato,asyto,csm(i,j),
1662f365b2 Andr*1819      *                 rr2t(i,j),tt2t(i,j),td2t(i,j))
                1820 
                1821            tauto=tausto+tauclf(i,j,k)
                1822            ssato=(ssatau+tauclf(i,j,k))/tauto+1.0e-8
32362b8fa8 Cons*1823            ssato=min(ssato,0.999999 _d 0)
1662f365b2 Andr*1824            asyto=(asysto+asycl(i,j,k)*tauclf(i,j,k))/(ssato*tauto)
                1825 
                1826            call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j))
                1827 
                1828           endif
                1829 
                1830          enddo
                1831         enddo
                1832 
                1833         do j=1,n
                1834          do i=1,m
                1835             rr(i,j,k,1)=rr1t(i,j)
                1836          enddo
                1837         enddo
                1838         do j=1,n
                1839          do i=1,m
                1840             tt(i,j,k,1)=tt1t(i,j)
                1841          enddo
                1842         enddo
                1843         do j=1,n
                1844          do i=1,m
                1845             td(i,j,k,1)=td1t(i,j)
                1846          enddo
                1847         enddo
                1848         do j=1,n
                1849          do i=1,m
                1850             rs(i,j,k,1)=rs1t(i,j)
                1851          enddo
                1852         enddo
                1853         do j=1,n
                1854          do i=1,m
                1855             ts(i,j,k,1)=ts1t(i,j)
                1856          enddo
                1857         enddo
                1858 
                1859         do j=1,n
                1860          do i=1,m
                1861             rr(i,j,k,2)=rr2t(i,j)
                1862          enddo
                1863         enddo
                1864         do j=1,n
                1865          do i=1,m
                1866             tt(i,j,k,2)=tt2t(i,j)
                1867          enddo
                1868         enddo
                1869         do j=1,n
                1870          do i=1,m
                1871             td(i,j,k,2)=td2t(i,j)
                1872          enddo
                1873         enddo
                1874         do j=1,n
                1875          do i=1,m
                1876             rs(i,j,k,2)=rs2t(i,j)
                1877          enddo
                1878         enddo
                1879         do j=1,n
                1880          do i=1,m
                1881             ts(i,j,k,2)=ts2t(i,j)
                1882          enddo
                1883         enddo
                1884 
                1885  300  continue
                1886 
                1887 c-----flux calculations
62e9d9f553 Jean*1888 
9524fe64b5 Andr*1889        do k= 1, np+1
                1890         do j= 1, n
                1891          do i= 1, m
                1892           fclr(i,j,k) = 0.
                1893           fall(i,j,k) = 0.
                1894          enddo
                1895         enddo
                1896        enddo
                1897        do j= 1, n
                1898         do i= 1, m
                1899          fsdir(i,j) = 0.
                1900          fsdif(i,j) = 0.
                1901         enddo
                1902        enddo
082e38725b Andr*1903         call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,
                1904      *               fclr,fall,fsdir,fsdif)
1662f365b2 Andr*1905 
                1906        do k= 1, np+1
                1907         do j= 1, n
                1908          do i= 1, m
                1909           flx(i,j,k)=flx(i,j,k)+fall(i,j,k)*hk(ib)
                1910          enddo
                1911         enddo
                1912         do j= 1, n
                1913          do i= 1, m
                1914           flc(i,j,k)=flc(i,j,k)+fclr(i,j,k)*hk(ib)
                1915          enddo
                1916         enddo
                1917        enddo
                1918 
                1919        if(ib.eq.nband) then
                1920          do j=1,n
                1921           do i=1,m
                1922            fdirpar(i,j)=fsdir(i,j)*hk(ib)
                1923            fdifpar(i,j)=fsdif(i,j)*hk(ib)
                1924          enddo
                1925         enddo
                1926        else
                1927          do j=1,n
                1928           do i=1,m
                1929            fdiruv(i,j)=fdiruv(i,j)+fsdir(i,j)*hk(ib)
                1930            fdifuv(i,j)=fdifuv(i,j)+fsdif(i,j)*hk(ib)
                1931          enddo
                1932         enddo
                1933        endif
                1934 
                1935  100  continue
                1936 
                1937       return
                1938       end
                1939 
                1940 c*********************************************************************
                1941 
                1942       subroutine deledd(tau,ssc,g0,csm,rr,tt,td)
                1943 
                1944 c*********************************************************************
                1945 c
                1946 c-----uses the delta-eddington approximation to compute the
                1947 c     bulk scattering properties of a single layer
                1948 c     coded following King and Harshvardhan (JAS, 1986)
                1949 c
                1950 c  inputs:
                1951 c
                1952 c     tau: the effective optical thickness
                1953 c     ssc: the effective single scattering albedo
                1954 c     g0:  the effective asymmetry factor
                1955 c     csm: the effective secant of the zenith angle
                1956 c
                1957 c  outputs:
                1958 c
                1959 c     rr: the layer reflection of the direct beam
                1960 c     tt: the layer diffuse transmission of the direct beam
                1961 c     td: the layer direct transmission of the direct beam
                1962 
                1963 c*********************************************************************
                1964 
                1965       implicit none
                1966 
62e9d9f553 Jean*1967 c-----Explicit Inline Directives
                1968 
06d4643e1f Jean*1969 #ifdef FIZHI_CRAY
4e1c053948 Jean*1970 #ifdef FIZHI_F77_COMPIL
1662f365b2 Andr*1971 cfpp$ expand (expmn)
62e9d9f553 Jean*1972 #endif
1662f365b2 Andr*1973 #endif
a456aa407c Andr*1974       _RL expmn
1662f365b2 Andr*1975 
a456aa407c Andr*1976       _RL  zero,one,two,three,four,fourth,seven,tumin
1662f365b2 Andr*1977       parameter (one=1., three=3.)
                1978       parameter (seven=7., two=2.)
                1979       parameter (four=4., fourth=.25)
                1980       parameter (zero=0., tumin=1.e-20)
                1981 
                1982 c-----input parameters
a456aa407c Andr*1983       _RL  tau,ssc,g0,csm
1662f365b2 Andr*1984 
                1985 c-----output parameters
a456aa407c Andr*1986       _RL  rr,tt,td
1662f365b2 Andr*1987 
                1988 c-----temporary parameters
                1989 
a456aa407c Andr*1990       _RL  zth,ff,xx,taup,sscp,gp,gm1,gm2,gm3,akk,alf1,alf2,
af35a8ab47 Andr*1991      *     all1,bll,st7,st8,cll,dll,fll,ell,st1,st2,st3,st4
1662f365b2 Andr*1992 c
                1993                 zth = one / csm
62e9d9f553 Jean*1994 
1662f365b2 Andr*1995 c  delta-eddington scaling of single scattering albedo,
                1996 c  optical thickness, and asymmetry factor,
                1997 c  K & H eqs(27-29)
                1998 
                1999                 ff  = g0*g0
                2000                 xx  = one-ff*ssc
                2001                 taup= tau*xx
                2002                 sscp= ssc*(one-ff)/xx
                2003                 gp  = g0/(one+g0)
62e9d9f553 Jean*2004 
1662f365b2 Andr*2005 c  extinction of the direct beam transmission
62e9d9f553 Jean*2006 
1662f365b2 Andr*2007                 td  = expmn(-taup*csm)
62e9d9f553 Jean*2008 
1662f365b2 Andr*2009 c  gamma1, gamma2, gamma3 and gamma4. see table 2 and eq(26) K & H
                2010 c  ssc and gp are the d-s single scattering
                2011 c  albedo and asymmetry factor.
                2012 
62e9d9f553 Jean*2013                 xx  =  three*gp
1662f365b2 Andr*2014                 gm1 =  (seven - sscp*(four+xx))*fourth
                2015                 gm2 = -(one   - sscp*(four-xx))*fourth
                2016                 gm3 =  (two   -        zth*xx )*fourth
62e9d9f553 Jean*2017 
1662f365b2 Andr*2018 c  akk is k as defined in eq(25) of K & H
62e9d9f553 Jean*2019 
1662f365b2 Andr*2020                 xx  = gm1-gm2
                2021                 akk = sqrt((gm1+gm2)*xx)
62e9d9f553 Jean*2022 
1662f365b2 Andr*2023 c   alf1 and alf2 are alpha1 and alpha2 from eqs (23) & (24) of K & H
62e9d9f553 Jean*2024 
1662f365b2 Andr*2025                 alf1 = gm1 - gm3 * xx
                2026                 alf2 = gm2 + gm3 * xx
62e9d9f553 Jean*2027 
af35a8ab47 Andr*2028 c  all1 is last term in eq(21) of K & H
1662f365b2 Andr*2029 c  bll is last term in eq(22) of K & H
62e9d9f553 Jean*2030 
1662f365b2 Andr*2031                 xx  = akk * two
af35a8ab47 Andr*2032                 all1 = (gm3 - alf2 * zth    )*xx*td
1662f365b2 Andr*2033                 bll = (one - gm3 + alf1*zth)*xx
62e9d9f553 Jean*2034 
1662f365b2 Andr*2035                 xx  = akk * zth
                2036                 st7 = one - xx
                2037                 st8 = one + xx
62e9d9f553 Jean*2038 
1662f365b2 Andr*2039                 xx  = akk * gm3
                2040                 cll = (alf2 + xx) * st7
                2041                 dll = (alf2 - xx) * st8
62e9d9f553 Jean*2042 
1662f365b2 Andr*2043                 xx  = akk * (one-gm3)
                2044                 fll = (alf1 + xx) * st8
                2045                 ell = (alf1 - xx) * st7
62e9d9f553 Jean*2046 
1662f365b2 Andr*2047                 st3 = max(abs(st7*st8),tumin)
                2048                 st3 = sign (st3,st7)
62e9d9f553 Jean*2049 
1662f365b2 Andr*2050                 st2 = expmn(-akk*taup)
                2051                 st4 = st2 * st2
62e9d9f553 Jean*2052 
1662f365b2 Andr*2053                 st1 =  sscp / ((akk+gm1 + (akk-gm1)*st4) * st3)
62e9d9f553 Jean*2054 
1662f365b2 Andr*2055 c   rr is r-hat of eq(21) of K & H
                2056 c   tt is diffuse part of t-hat of eq(22) of K & H
62e9d9f553 Jean*2057 
af35a8ab47 Andr*2058                 rr =   ( cll-dll*st4    -all1*st2)*st1
1662f365b2 Andr*2059                 tt = - ((fll-ell*st4)*td-bll*st2)*st1
62e9d9f553 Jean*2060 
1662f365b2 Andr*2061                 rr = max(rr,zero)
                2062                 tt = max(tt,zero)
                2063 
                2064       return
                2065       end
                2066 
                2067 c*********************************************************************
                2068 
                2069       subroutine sagpol(tau,ssc,g0,rll,tll)
                2070 
                2071 c*********************************************************************
                2072 c-----transmittance (tll) and reflectance (rll) of diffuse radiation
                2073 c     follows Sagan and Pollock (JGR, 1967).
                2074 c     also, eq.(31) of Lacis and Hansen (JAS, 1974).
                2075 c
                2076 c-----input parameters:
                2077 c
                2078 c      tau: the effective optical thickness
                2079 c      ssc: the effective single scattering albedo
                2080 c      g0:  the effective asymmetry factor
                2081 c
                2082 c-----output parameters:
                2083 c
                2084 c      rll: the layer reflection of diffuse radiation
                2085 c      tll: the layer transmission of diffuse radiation
                2086 c
                2087 c*********************************************************************
                2088 
                2089       implicit none
                2090 
62e9d9f553 Jean*2091 c-----Explicit Inline Directives
                2092 
06d4643e1f Jean*2093 #ifdef FIZHI_CRAY
4e1c053948 Jean*2094 #ifdef FIZHI_F77_COMPIL
1662f365b2 Andr*2095 cfpp$ expand (expmn)
62e9d9f553 Jean*2096 #endif
1662f365b2 Andr*2097 #endif
a456aa407c Andr*2098       _RL expmn
1662f365b2 Andr*2099 
a456aa407c Andr*2100       _RL  one,three,four
1662f365b2 Andr*2101       parameter (one=1., three=3., four=4.)
                2102 
                2103 c-----output parameters:
                2104 
a456aa407c Andr*2105       _RL  tau,ssc,g0
1662f365b2 Andr*2106 
                2107 c-----output parameters:
                2108 
a456aa407c Andr*2109       _RL  rll,tll
1662f365b2 Andr*2110 
                2111 c-----temporary arrays
                2112 
a456aa407c Andr*2113       _RL  xx,uuu,ttt,emt,up1,um1,st1
1662f365b2 Andr*2114 
                2115              xx  = one-ssc*g0
                2116              uuu = sqrt( xx/(one-ssc))
                2117              ttt = sqrt( xx*(one-ssc)*three )*tau
                2118              emt = expmn(-ttt)
                2119              up1 = uuu + one
                2120              um1 = uuu - one
                2121              xx  = um1*emt
                2122              st1 = one / ((up1+xx) * (up1-xx))
                2123              rll = up1*um1*(one-emt*emt)*st1
                2124              tll = uuu*four*emt         *st1
                2125 
                2126       return
                2127       end
                2128 
                2129 c*******************************************************************
                2130 
                2131       function expmn(fin)
                2132 
                2133 c*******************************************************************
                2134 c compute exponential for arguments in the range 0> fin > -10.
175684e43e Andr*2135 c*******************************************************************
                2136       implicit none
a456aa407c Andr*2137       _RL  fin,expmn
1662f365b2 Andr*2138 
a456aa407c Andr*2139       _RL one,expmin,e1,e2,e3,e4
1662f365b2 Andr*2140       parameter (one=1.0, expmin=-10.0)
                2141       parameter (e1=1.0,        e2=-2.507213e-1)
                2142       parameter (e3=2.92732e-2, e4=-3.827800e-3)
                2143 
                2144       if (fin .lt. expmin) fin = expmin
                2145       expmn = ((e4*fin + e3)*fin+e2)*fin+e1
                2146       expmn = expmn * expmn
                2147       expmn = one / (expmn * expmn)
                2148 
                2149       return
                2150       end
                2151 
                2152 c*******************************************************************
                2153 
                2154       subroutine cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,
                2155      *           fclr,fall,fsdir,fsdif)
                2156 
                2157 c*******************************************************************
                2158 c  compute upward and downward fluxes using a two-stream adding method
                2159 c  following equations (3)-(5) of Chou (1992, JAS).
                2160 c
62e9d9f553 Jean*2161 c  clouds are grouped into high, middle, and low clouds which are
1662f365b2 Andr*2162 c  assumed randomly overlapped. It involves eight sets of calculations.
                2163 c  In each set of calculations, each atmospheric layer is homogeneous,
                2164 c  either with clouds or without clouds.
                2165 
                2166 c  input parameters:
                2167 c     m:   number of soundings in zonal direction
                2168 c     n:   number of soundings in meridional direction
                2169 c     np:  number of atmospheric layers
                2170 c     ict: the level separating high and middle clouds
                2171 c     icb: the level separating middle and low clouds
                2172 c     cc:  effective cloud covers for high, middle and low clouds
                2173 c     tt:  diffuse transmission of a layer illuminated by beam radiation
                2174 c     td:  direct beam tranmssion
                2175 c     ts:  transmission of a layer illuminated by diffuse radiation
                2176 c     rr:  reflection of a layer illuminated by beam radiation
                2177 c     rs:  reflection of a layer illuminated by diffuse radiation
                2178 c
                2179 c  output parameters:
                2180 c     fclr:  clear-sky flux (downward minus upward)
                2181 c     fall:  all-sky flux (downward minus upward)
                2182 c     fdndir: surface direct downward flux
                2183 c     fdndif: surface diffuse downward flux
                2184 c*********************************************************************c
                2185 
                2186       implicit none
                2187 
                2188 c-----input parameters
                2189 
                2190       integer m,n,np,ict,icb
                2191 
a456aa407c Andr*2192       _RL  rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2)
                2193       _RL  rs(m,n,np+1,2),ts(m,n,np+1,2)
                2194       _RL  cc(m,n,3)
1662f365b2 Andr*2195 
                2196 c-----temporary array
                2197 
                2198       integer i,j,k,ih,im,is
a456aa407c Andr*2199       _RL  rra(m,n,np+1,2,2),tta(m,n,np+1,2,2),tda(m,n,np+1,2,2)
                2200       _RL  rsa(m,n,np+1,2,2),rxa(m,n,np+1,2,2)
                2201       _RL  ch(m,n),cm(m,n),ct(m,n),flxdn(m,n,np+1)
                2202       _RL  fdndir(m,n),fdndif(m,n),fupdif
                2203       _RL  denm,xx
1662f365b2 Andr*2204 
                2205 c-----output parameters
                2206 
a456aa407c Andr*2207       _RL  fclr(m,n,np+1),fall(m,n,np+1)
                2208       _RL  fsdir(m,n),fsdif(m,n)
1662f365b2 Andr*2209 
                2210 c-----initialize all-sky flux (fall) and surface downward fluxes
                2211 
                2212       do k=1,np+1
                2213        do j=1,n
                2214         do i=1,m
                2215            fall(i,j,k)=0.0
                2216         enddo
                2217        enddo
                2218       enddo
                2219 
                2220        do j=1,n
                2221         do i=1,m
                2222            fsdir(i,j)=0.0
                2223            fsdif(i,j)=0.0
                2224         enddo
                2225        enddo
                2226 
                2227 c-----compute transmittances and reflectances for a composite of
                2228 c     layers. layers are added one at a time, going down from the top.
                2229 c     tda is the composite transmittance illuminated by beam radiation
                2230 c     tta is the composite diffuse transmittance illuminated by
                2231 c         beam radiation
                2232 c     rsa is the composite reflectance illuminated from below
                2233 c         by diffuse radiation
                2234 c     tta and rsa are computed from eqs. (4b) and (3b) of Chou
                2235 
                2236 c-----for high clouds. indices 1 and 2 denote clear and cloudy
                2237 c     situations, respectively.
                2238 
                2239       do 10 ih=1,2
                2240 
                2241        do j= 1, n
                2242         do i= 1, m
                2243           tda(i,j,1,ih,1)=td(i,j,1,ih)
                2244           tta(i,j,1,ih,1)=tt(i,j,1,ih)
                2245           rsa(i,j,1,ih,1)=rs(i,j,1,ih)
                2246           tda(i,j,1,ih,2)=td(i,j,1,ih)
                2247           tta(i,j,1,ih,2)=tt(i,j,1,ih)
                2248           rsa(i,j,1,ih,2)=rs(i,j,1,ih)
                2249         enddo
                2250        enddo
                2251 
                2252        do k= 2, ict-1
                2253         do j= 1, n
                2254          do i= 1, m
                2255           denm = ts(i,j,k,ih)/( 1.-rsa(i,j,k-1,ih,1)*rs(i,j,k,ih))
                2256           tda(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*td(i,j,k,ih)
                2257           tta(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*tt(i,j,k,ih)
                2258      *                  +(tda(i,j,k-1,ih,1)*rr(i,j,k,ih)
                2259      *                  *rsa(i,j,k-1,ih,1)+tta(i,j,k-1,ih,1))*denm
                2260           rsa(i,j,k,ih,1)= rs(i,j,k,ih)+ts(i,j,k,ih)
                2261      *                  *rsa(i,j,k-1,ih,1)*denm
                2262           tda(i,j,k,ih,2)= tda(i,j,k,ih,1)
                2263           tta(i,j,k,ih,2)= tta(i,j,k,ih,1)
                2264           rsa(i,j,k,ih,2)= rsa(i,j,k,ih,1)
                2265          enddo
                2266         enddo
                2267        enddo
                2268 
                2269 c-----for middle clouds
                2270 
                2271       do 10 im=1,2
                2272 
                2273        do k= ict, icb-1
                2274         do j= 1, n
                2275          do i= 1, m
                2276           denm = ts(i,j,k,im)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,im))
                2277           tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,im)
                2278           tta(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*tt(i,j,k,im)
                2279      *                  +(tda(i,j,k-1,ih,im)*rr(i,j,k,im)
                2280      *                  *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
                2281           rsa(i,j,k,ih,im)= rs(i,j,k,im)+ts(i,j,k,im)
                2282      *                  *rsa(i,j,k-1,ih,im)*denm
                2283          enddo
                2284         enddo
                2285        enddo
                2286 
                2287   10  continue
                2288 
                2289 c-----layers are added one at a time, going up from the surface.
                2290 c     rra is the composite reflectance illuminated by beam radiation
                2291 c     rxa is the composite reflectance illuminated from above
                2292 c         by diffuse radiation
                2293 c     rra and rxa are computed from eqs. (4a) and (3a) of Chou
62e9d9f553 Jean*2294 
1662f365b2 Andr*2295 c-----for the low clouds
                2296 
                2297       do 20 is=1,2
                2298 
                2299        do j= 1, n
                2300         do i= 1, m
                2301          rra(i,j,np+1,1,is)=rr(i,j,np+1,is)
                2302          rxa(i,j,np+1,1,is)=rs(i,j,np+1,is)
                2303          rra(i,j,np+1,2,is)=rr(i,j,np+1,is)
                2304          rxa(i,j,np+1,2,is)=rs(i,j,np+1,is)
                2305         enddo
                2306        enddo
                2307 
                2308        do k=np,icb,-1
                2309         do j= 1, n
                2310          do i= 1, m
                2311           denm=ts(i,j,k,is)/( 1.-rs(i,j,k,is)*rxa(i,j,k+1,1,is) )
                2312           rra(i,j,k,1,is)=rr(i,j,k,is)+(td(i,j,k,is)
                2313      *        *rra(i,j,k+1,1,is)+tt(i,j,k,is)*rxa(i,j,k+1,1,is))*denm
                2314           rxa(i,j,k,1,is)= rs(i,j,k,is)+ts(i,j,k,is)
                2315      *        *rxa(i,j,k+1,1,is)*denm
                2316           rra(i,j,k,2,is)=rra(i,j,k,1,is)
                2317           rxa(i,j,k,2,is)=rxa(i,j,k,1,is)
                2318          enddo
                2319         enddo
                2320        enddo
                2321 
                2322 c-----for middle clouds
                2323 
                2324       do 20 im=1,2
                2325 
                2326        do k= icb-1,ict,-1
                2327         do j= 1, n
                2328          do i= 1, m
                2329           denm=ts(i,j,k,im)/( 1.-rs(i,j,k,im)*rxa(i,j,k+1,im,is) )
                2330           rra(i,j,k,im,is)= rr(i,j,k,im)+(td(i,j,k,im)
                2331      *        *rra(i,j,k+1,im,is)+tt(i,j,k,im)*rxa(i,j,k+1,im,is))*denm
                2332           rxa(i,j,k,im,is)= rs(i,j,k,im)+ts(i,j,k,im)
                2333      *        *rxa(i,j,k+1,im,is)*denm
                2334          enddo
                2335         enddo
                2336        enddo
                2337 
                2338   20  continue
                2339 
                2340 c-----integration over eight sky situations.
                2341 c     ih, im, is denotes high, middle and low cloud groups.
                2342 
                2343       do 100 ih=1,2
                2344 
62e9d9f553 Jean*2345 c-----clear portion
1662f365b2 Andr*2346 
                2347          if(ih.eq.1) then
                2348            do j=1,n
                2349             do i=1,m
                2350              ch(i,j)=1.0-cc(i,j,1)
                2351             enddo
                2352            enddo
                2353 
                2354           else
                2355 
                2356 c-----cloudy portion
                2357 
                2358            do j=1,n
                2359             do i=1,m
                2360              ch(i,j)=cc(i,j,1)
                2361             enddo
                2362            enddo
                2363 
                2364           endif
                2365 
                2366       do 100 im=1,2
                2367 
                2368 c-----clear portion
                2369 
                2370          if(im.eq.1) then
                2371 
                2372            do j=1,n
                2373             do i=1,m
                2374               cm(i,j)=ch(i,j)*(1.0-cc(i,j,2))
                2375             enddo
                2376            enddo
                2377 
                2378          else
                2379 
                2380 c-----cloudy portion
                2381 
                2382            do j=1,n
                2383             do i=1,m
62e9d9f553 Jean*2384               cm(i,j)=ch(i,j)*cc(i,j,2)
1662f365b2 Andr*2385             enddo
                2386            enddo
                2387 
                2388          endif
                2389 
                2390       do 100 is=1,2
                2391 
                2392 c-----clear portion
                2393 
                2394          if(is.eq.1) then
                2395 
                2396            do j=1,n
                2397             do i=1,m
62e9d9f553 Jean*2398              ct(i,j)=cm(i,j)*(1.0-cc(i,j,3))
1662f365b2 Andr*2399             enddo
                2400            enddo
                2401 
                2402          else
                2403 
                2404 c-----cloudy portion
                2405 
                2406            do j=1,n
                2407             do i=1,m
                2408              ct(i,j)=cm(i,j)*cc(i,j,3)
                2409             enddo
                2410            enddo
                2411 
                2412          endif
                2413 
                2414 c-----add one layer at a time, going down.
                2415 
                2416        do k= icb, np
                2417         do j= 1, n
                2418          do i= 1, m
                2419           denm = ts(i,j,k,is)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,is) )
                2420           tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,is)
                2421           tta(i,j,k,ih,im)=  tda(i,j,k-1,ih,im)*tt(i,j,k,is)
                2422      *         +(tda(i,j,k-1,ih,im)*rr(i,j,k,is)
                2423      *         *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
                2424           rsa(i,j,k,ih,im)= rs(i,j,k,is)+ts(i,j,k,is)
                2425      *         *rsa(i,j,k-1,ih,im)*denm
                2426          enddo
                2427         enddo
                2428        enddo
                2429 
                2430 c-----add one layer at a time, going up.
                2431 
                2432        do k= ict-1,1,-1
                2433         do j= 1, n
                2434          do i= 1, m
                2435           denm =ts(i,j,k,ih)/(1.-rs(i,j,k,ih)*rxa(i,j,k+1,im,is))
                2436           rra(i,j,k,im,is)= rr(i,j,k,ih)+(td(i,j,k,ih)
                2437      *        *rra(i,j,k+1,im,is)+tt(i,j,k,ih)*rxa(i,j,k+1,im,is))*denm
                2438           rxa(i,j,k,im,is)= rs(i,j,k,ih)+ts(i,j,k,ih)
                2439      *        *rxa(i,j,k+1,im,is)*denm
                2440          enddo
                2441         enddo
                2442        enddo
                2443 
                2444 c-----compute fluxes following eq (5) of Chou (1992)
62e9d9f553 Jean*2445 
1662f365b2 Andr*2446 c     fdndir is the direct  downward flux
                2447 c     fdndif is the diffuse downward flux
                2448 c     fupdif is the diffuse upward flux
                2449 
                2450       do k=2,np+1
                2451        do j=1, n
                2452         do i=1, m
                2453          denm= 1./(1.- rxa(i,j,k,im,is)*rsa(i,j,k-1,ih,im))
                2454          fdndir(i,j)= tda(i,j,k-1,ih,im)
                2455          xx = tda(i,j,k-1,ih,im)*rra(i,j,k,im,is)
                2456          fdndif(i,j)= (xx*rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
                2457          fupdif= (xx+tta(i,j,k-1,ih,im)*rxa(i,j,k,im,is))*denm
                2458          flxdn(i,j,k)=fdndir(i,j)+fdndif(i,j)-fupdif
                2459         enddo
                2460        enddo
                2461       enddo
                2462 
                2463        do j=1, n
                2464         do i=1, m
                2465          flxdn(i,j,1)=1.0-rra(i,j,1,im,is)
                2466         enddo
                2467        enddo
                2468 
                2469 c-----summation of fluxes over all (eight) sky situations.
                2470 
                2471        do k=1,np+1
                2472         do j=1,n
                2473          do i=1,m
                2474            if(ih.eq.1 .and. im.eq.1 .and. is.eq.1) then
                2475              fclr(i,j,k)=flxdn(i,j,k)
                2476            endif
                2477              fall(i,j,k)=fall(i,j,k)+flxdn(i,j,k)*ct(i,j)
                2478          enddo
                2479         enddo
                2480        enddo
                2481 
                2482         do j=1,n
                2483          do i=1,m
                2484             fsdir(i,j)=fsdir(i,j)+fdndir(i,j)*ct(i,j)
                2485             fsdif(i,j)=fsdif(i,j)+fdndif(i,j)*ct(i,j)
                2486          enddo
                2487         enddo
                2488 
                2489   100 continue
                2490 
                2491       return
                2492       end
                2493 
                2494 c*****************************************************************
                2495 
                2496       subroutine flxco2(m,n,np,swc,swh,csm,df)
                2497 
                2498 c*****************************************************************
                2499 
                2500 c-----compute the reduction of clear-sky downward solar flux
                2501 c     due to co2 absorption.
                2502 
                2503       implicit none
                2504 
                2505 c-----input parameters
                2506 
                2507       integer m,n,np
a456aa407c Andr*2508       _RL  csm(m,n),swc(m,n,np+1),swh(m,n,np+1),cah(22,19)
1662f365b2 Andr*2509 
                2510 c-----output (undated) parameter
                2511 
a456aa407c Andr*2512       _RL  df(m,n,np+1)
1662f365b2 Andr*2513 
                2514 c-----temporary array
                2515 
62e9d9f553 Jean*2516       integer i,j,k,ic,iw
af35a8ab47 Andr*2517       _RL  xx,clog1,wlog,dc,dw,x1,x2,y2
1662f365b2 Andr*2518 
                2519 c********************************************************************
                2520 c-----include co2 look-up table
                2521 
4a402f223d Andr*2522 #include "cah-dat.h"
7f6041e4b1 Andr*2523 c     save cah
62e9d9f553 Jean*2524 
1662f365b2 Andr*2525 c********************************************************************
                2526 c-----table look-up for the reduction of clear-sky solar
                2527 c     radiation due to co2. The fraction 0.0343 is the
                2528 c     extraterrestrial solar flux in the co2 bands.
                2529 
                2530       do k= 2, np+1
                2531        do j= 1, n
                2532         do i= 1, m
                2533           xx=1./.3
af35a8ab47 Andr*2534           clog1=log10(swc(i,j,k)*csm(i,j))
1662f365b2 Andr*2535           wlog=log10(swh(i,j,k)*csm(i,j))
af35a8ab47 Andr*2536           ic=int( (clog1+3.15)*xx+1.)
1662f365b2 Andr*2537           iw=int( (wlog+4.15)*xx+1.)
                2538           if(ic.lt.2)ic=2
                2539           if(iw.lt.2)iw=2
                2540           if(ic.gt.22)ic=22
62e9d9f553 Jean*2541           if(iw.gt.19)iw=19
af35a8ab47 Andr*2542           dc=clog1-float(ic-2)*.3+3.
62e9d9f553 Jean*2543           dw=wlog-float(iw-2)*.3+4.
1662f365b2 Andr*2544           x1=cah(1,iw-1)+(cah(1,iw)-cah(1,iw-1))*xx*dw
                2545           x2=cah(ic-1,iw-1)+(cah(ic-1,iw)-cah(ic-1,iw-1))*xx*dw
                2546           y2=x2+(cah(ic,iw-1)-cah(ic-1,iw-1))*xx*dc
                2547           if (x1.lt.y2) x1=y2
                2548           df(i,j,k)=df(i,j,k)+0.0343*(x1-y2)
62e9d9f553 Jean*2549         enddo
                2550        enddo
                2551       enddo
1662f365b2 Andr*2552 
                2553       return
                2554       end
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!