/[MITgcm]/MITgcm/pkg/fizhi/fizhi_swrad.F
ViewVC logotype

Diff of /MITgcm/pkg/fizhi/fizhi_swrad.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.4 by molod, Thu Jun 24 19:57:02 2004 UTC revision 1.22 by molod, Tue Dec 14 19:56:45 2004 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "FIZHI_OPTIONS.h"
5        subroutine swrio (nymd,nhms,bi,bj,ndswr,myid,istrip,npcs,        subroutine swrio (nymd,nhms,bi,bj,ndswr,myid,istrip,npcs,
6       .        pz,tz,qz,pkht,oz,co2,       .        low_level,mid_level,im,jm,lm,
7         .        pz,plz,plze,dpres,pkht,pkz,tz,qz,oz,co2,
8       .        albvisdr,albvisdf,albnirdr,albnirdf,       .        albvisdr,albvisdf,albnirdr,albnirdf,
9       .        dtradsw,dtswclr,radswg,swgclr,       .        dtradsw,dtswclr,radswg,swgclr,
10       .        fdifpar,fdirpar,osr,osrclr,       .        fdifpar,fdirpar,osr,osrclr,
11       .        im,jm,lm,sige,sig,dsig,ptop,       .        ptop,nswcld,cldsw,cswmo,nswlz,swlz,
      .        nswcld,cldsw,cswmo,nswlz,swlz,  
12       .        lpnt,imstturb,qliqave,fccave,landtype,xlats,xlons)       .        lpnt,imstturb,qliqave,fccave,landtype,xlats,xlons)
13    
14        implicit none        implicit none
15  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
16  #include "diagnostics.h"  #include "SIZE.h"
17    #include "DIAGNOSTICS_SIZE.h"
18    #include "DIAGNOSTICS.h"
19  #endif  #endif
20    
21  c Input Variables  c Input Variables
22  c ---------------  c ---------------
23        integer nymd,nhms,ndswr,istrip,npcs,bi,bj        integer nymd,nhms,bi,bj,ndswr,myid,istrip,npcs
24          integer mid_level,low_level
25        integer im,jm,lm                integer im,jm,lm        
26        real  ptop                      _RL  ptop
27        real  sige(lm+1)                _RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1),dpres(im,jm,lm)
28        real   sig(lm)                  _RL pkht(im,jm,lm+1),pkz(im,jm,lm)
29        real  dsig(lm)                  _RL tz(im,jm,lm),qz(im,jm,lm)
30          _RL oz(im,jm,lm)
31        real    pz(im,jm)              _RL co2
32        real    tz(im,jm,lm)            _RL albvisdr(im,jm),albvisdf(im,jm),albnirdr(im,jm)
33        real  pkht(im,jm,lm)            _RL albnirdf(im,jm)
34          _RL radswg(im,jm),swgclr(im,jm),fdifpar(im,jm),fdirpar(im,jm)
35        real    co2                    _RL osr(im,jm),osrclr(im,jm),dtradsw(im,jm,lm)
36        real    oz(im,jm,lm)            _RL dtswclr(im,jm,lm)
       real    qz(im,jm,lm)      
   
       real albvisdr(im,jm)      
       real albvisdf(im,jm)      
       real albnirdr(im,jm)      
       real albnirdf(im,jm)      
   
       real   radswg(im,jm)      
       real   swgclr(im,jm)      
       real  fdifpar(im,jm)      
       real  fdirpar(im,jm)      
       real      osr(im,jm)      
       real   osrclr(im,jm)      
       real  dtradsw(im,jm,lm)  
       real  dtswclr(im,jm,lm)  
   
37        integer nswcld,nswlz            integer nswcld,nswlz    
38        real  cldsw(im,jm,lm)          _RL cldsw(im,jm,lm),cswmo(im,jm,lm),swlz(im,jm,lm)  
       real  cswmo(im,jm,lm)    
       real   swlz(im,jm,lm)    
   
39        logical lpnt                    logical lpnt            
40        integer imstturb                integer imstturb        
41        real qliqave(im,jm,lm)          _RL qliqave(im,jm,lm),fccave(im,jm,lm)  
       real  fccave(im,jm,lm)    
   
42        integer landtype(im,jm)        integer landtype(im,jm)
43          _RL xlats(im,jm),xlons(im,jm)
44    
45  c Local Variables  c Local Variables
46  c ---------------  c ---------------
47        integer   i,j,L,nn,nsecf,mid_level, low_level        integer   i,j,L,nn,nsecf
48        integer   nb2,ntmstp,nymd2,nhms2        integer   ntmstp,nymd2,nhms2
49        real      getcon,grav,cp,undef,pcheck        _RL      getcon,grav,cp,undef
50        real      ra,alf,reffw,reffi,tminv        _RL      ra,alf,reffw,reffi,tminv
51    
52        parameter ( reffw = 10.0 )          parameter ( reffw = 10.0 )  
53        parameter ( reffi = 65.0 )          parameter ( reffi = 65.0 )  
54    
55        real      alat(im,jm)        _RL tdry(im,jm,lm)
56        real      alon(im,jm)        _RL TEMP1(im,jm)
57          _RL TEMP2(im,jm)
58          _RL zenith (im,jm)
59          _RL cldtot (im,jm,lm)
60          _RL cldmxo (im,jm,lm)
61          _RL totcld (im,jm)
62          _RL cldlow (im,jm)
63          _RL cldmid (im,jm)
64          _RL cldhi  (im,jm)
65          _RL taulow (im,jm)
66          _RL taumid (im,jm)
67          _RL tauhi  (im,jm)
68          _RL tautype(im,jm,lm,3)
69          _RL tau(im,jm,lm)
70          _RL albedo(im,jm)    
71    
72          _RL PK(ISTRIP,lm)
73          _RL qzl(ISTRIP,lm),CLRO(ISTRIP,lm)
74          _RL TZL(ISTRIP,lm)
75          _RL OZL(ISTRIP,lm)
76          _RL PLE(ISTRIP,lm+1)
77          _RL COSZ(ISTRIP)
78          _RL dpstrip(ISTRIP,lm)
79    
80          _RL albuvdr(ISTRIP),albuvdf(ISTRIP)
81          _RL albirdr(ISTRIP),albirdf(ISTRIP)
82          _RL difpar (ISTRIP),dirpar (ISTRIP)
83    
84          _RL fdirir(istrip),fdifir(istrip)
85          _RL fdiruv(istrip),fdifuv(istrip)
86    
87          _RL flux(istrip,lm+1)
88          _RL fluxclr(istrip,lm+1)
89          _RL dtsw(istrip,lm)
90          _RL dtswc(istrip,lm)
91    
92          _RL taul   (istrip,lm)
93          _RL reff   (istrip,lm,2)
94          _RL tauc   (istrip,lm,2)
95          _RL taua   (istrip,lm)
96          _RL tstrip (istrip)
97    
98        real          PKZ(im,jm,lm)        logical first
99        real          PLZ(im,jm,lm)        data first /.true./
       real         tdry(im,jm,lm)  
       real         PLZE(im,jm,lm+1)  
       real        TEMP1(im,jm)  
       real        TEMP2(im,jm)  
       real      zenith (im,jm)  
       real      cldtot (im,jm,lm)  
       real      cldmxo (im,jm,lm)  
       real      totcld (im,jm)  
       real      cldlow (im,jm)  
       real      cldmid (im,jm)  
       real      cldhi  (im,jm)  
       real      taulow (im,jm)  
       real      taumid (im,jm)  
       real      tauhi  (im,jm)  
       real      tautype(im,jm,lm,3)  
       real      tau    (im,jm,lm)  
       real      albedo(im,jm)      
   
       real          PK(ISTRIP,lm)  
       real         qzl(ISTRIP,lm),  CLRO(ISTRIP,lm)  
       real         TZL(ISTRIP,lm)  
       real         OZL(ISTRIP,lm)  
       real         PLE(ISTRIP,lm+1)  
       real        COSZ(ISTRIP)  
   
       real      albuvdr(ISTRIP),albuvdf(ISTRIP)  
       real      albirdr(ISTRIP),albirdf(ISTRIP)  
       real      difpar (ISTRIP),dirpar (ISTRIP)  
   
       real      fdirir(istrip),fdifir(istrip)  
       real      fdiruv(istrip),fdifuv(istrip)  
   
       real      flux   (istrip,lm+1)  
       real      fluxclr(istrip,lm+1)  
       real      dtsw   (istrip,lm)  
       real      dtswc  (istrip,lm)  
   
       real      taul   (istrip,lm)  
       real      reff   (istrip,lm,2)  
       real      tauc   (istrip,lm,2)  
       real      taua   (istrip,lm)  
       real      tstrip (istrip)  
   
       logical   first  
       data      first /.true./  
   
       integer   koz, kh2o  
       data      KOZ  /20/  
       data      kh2o /18/  
100    
101  C **********************************************************************  C **********************************************************************
102  C ****                       INITIALIZATION                         ****  C ****                       INITIALIZATION                         ****
# Line 136  C ************************************** Line 109  C **************************************
109        NTMSTP = nsecf(NDSWR)        NTMSTP = nsecf(NDSWR)
110        TMINV  = 1./float(ntmstp)        TMINV  = 1./float(ntmstp)
111    
       do j = 1,jm  
       do i = 1,im  
       PLZE(I,j,1) = SIGE(1)*PZ(I,j) + PTOP  
       enddo  
       enddo  
         
       DO L = 1,lm  
       do j = 1,jm  
       DO I = 1,im  
       PLZ (I,j,L  ) = SIG (L)  *PZ(I,j) + PTOP  
       PLZE(I,j,L+1) = SIGE(L+1)*PZ(I,j) + PTOP  
       ENDDO  
       ENDDO  
       ENDDO  
   
       call pkappa ( pz,pkht,pkz,ptop,sige,dsig,im,jm,lm )  
   
112  C Compute Temperature from Theta  C Compute Temperature from Theta
113  C ------------------------------  C ------------------------------
114        do L=1,lm        do L=1,lm
# Line 163  C ------------------------------ Line 119  C ------------------------------
119        enddo        enddo
120        enddo        enddo
121    
122  c Determine Level Indices for Low-Mid-High Cloud Regions        if (first .and. myid.eq.1 ) then
 c ------------------------------------------------------  
       low_level = lm  
       mid_level = lm  
       do L = lm-1,1,-1  
       pcheck = (1000.-ptop)*sig(l) + ptop  
       if (pcheck.gt.700.0) low_level = L  
       if (pcheck.gt.400.0) mid_level = L  
       enddo  
   
       if (first .and. myid.eq.0 ) then  
123        print *        print *
124        print *,'Low-Level Clouds are Grouped between levels: ',        print *,'Low-Level Clouds are Grouped between levels: ',
125       .         lm,' and ',low_level       .         lm,' and ',low_level
# Line 187  C ************************************** Line 133  C **************************************
133  C ****             CALCULATE COSINE OF THE ZENITH ANGLE             ****  C ****             CALCULATE COSINE OF THE ZENITH ANGLE             ****
134  C **********************************************************************  C **********************************************************************
135    
136        CALL ASTRO ( NYMD,   NHMS,  ALAT,ALON, im*jm, TEMP1,RA )        CALL ASTRO ( NYMD,   NHMS,  XLATS,XLONS, im*jm, TEMP1,RA )
137                     NYMD2 = NYMD                     NYMD2 = NYMD
138                     NHMS2 = NHMS                     NHMS2 = NHMS
139        CALL TICK  ( NYMD2,  NHMS2, NTMSTP )        CALL TICK  ( NYMD2,  NHMS2, NTMSTP )
140        CALL ASTRO ( NYMD2,  NHMS2, ALAT,ALON, im*jm, TEMP2,RA )        CALL ASTRO ( NYMD2,  NHMS2, XLATS,XLONS, im*jm, TEMP2,RA )
141    
142        do j = 1,jm        do j = 1,jm
143        do i = 1,im        do i = 1,im
# Line 203  C ************************************** Line 149  C **************************************
149        ENDDO        ENDDO
150        ENDDO        ENDDO
151    
   
152  C **********************************************************************  C **********************************************************************
153  c ****        Compute Two-Dimension Total Cloud Fraction (0-1)      ****  c ****        Compute Two-Dimension Total Cloud Fraction (0-1)      ****
154  C **********************************************************************  C **********************************************************************
# Line 281  c ------------------------- Line 226  c -------------------------
226        qdiag(i,j,icldfrc,bi,bj) =  qdiag(i,j,icldfrc,bi,bj) + totcld(i,j)        qdiag(i,j,icldfrc,bi,bj) =  qdiag(i,j,icldfrc,bi,bj) + totcld(i,j)
227        enddo        enddo
228        enddo        enddo
229        ncldfrc = ncldfrc + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) ncldfrc = ncldfrc + 1
230        endif        endif
231    
232        if( icldras.gt.0 ) then        if( icldras.gt.0 ) then
# Line 293  c ------------------------- Line 238  c -------------------------
238        enddo        enddo
239        enddo        enddo
240        enddo        enddo
241        ncldras = ncldras + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) ncldras = ncldras + 1
242        endif        endif
243    
244        if( icldtot.gt.0 ) then        if( icldtot.gt.0 ) then
# Line 305  c ------------------------- Line 250  c -------------------------
250        enddo        enddo
251        enddo        enddo
252        enddo        enddo
253        ncldtot = ncldtot + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) ncldtot = ncldtot + 1
254        endif        endif
255    
256        if( icldlow.gt.0 ) then        if( icldlow.gt.0 ) then
# Line 314  c ------------------------- Line 259  c -------------------------
259        qdiag(i,j,icldlow,bi,bj) = qdiag(i,j,icldlow,bi,bj) + cldlow(i,j)        qdiag(i,j,icldlow,bi,bj) = qdiag(i,j,icldlow,bi,bj) + cldlow(i,j)
260        enddo        enddo
261        enddo        enddo
262        ncldlow = ncldlow + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) ncldlow = ncldlow + 1
263        endif        endif
264    
265        if( icldmid.gt.0 ) then        if( icldmid.gt.0 ) then
# Line 323  c ------------------------- Line 268  c -------------------------
268        qdiag(i,j,icldmid,bi,bj) = qdiag(i,j,icldmid,bi,bj) + cldmid(i,j)        qdiag(i,j,icldmid,bi,bj) = qdiag(i,j,icldmid,bi,bj) + cldmid(i,j)
269        enddo        enddo
270        enddo        enddo
271        ncldmid = ncldmid + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) ncldmid = ncldmid + 1
272        endif        endif
273    
274        if( icldhi.gt.0 ) then        if( icldhi.gt.0 ) then
# Line 332  c ------------------------- Line 277  c -------------------------
277        qdiag(i,j,icldhi,bi,bj) = qdiag(i,j,icldhi,bi,bj) + cldhi(i,j)        qdiag(i,j,icldhi,bi,bj) = qdiag(i,j,icldhi,bi,bj) + cldhi(i,j)
278        enddo        enddo
279        enddo        enddo
280        ncldhi = ncldhi + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) ncldhi = ncldhi + 1
281        endif        endif
282    
283        if( ilzrad.gt.0 ) then        if( ilzrad.gt.0 ) then
# Line 344  c ------------------------- Line 289  c -------------------------
289        enddo        enddo
290        enddo        enddo
291        enddo        enddo
292        nlzrad = nlzrad + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) nlzrad = nlzrad + 1
293        endif        endif
294    
295  c Albedo Diagnostics  c Albedo Diagnostics
# Line 356  c ------------------ Line 301  c ------------------
301       .                                                     albvisdr(i,j)       .                                                     albvisdr(i,j)
302        enddo        enddo
303        enddo        enddo
304        nalbvisdr = nalbvisdr + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) nalbvisdr = nalbvisdr + 1
305        endif        endif
306    
307        if( ialbvisdf.gt.0 ) then        if( ialbvisdf.gt.0 ) then
# Line 366  c ------------------ Line 311  c ------------------
311       .                                                     albvisdf(i,j)       .                                                     albvisdf(i,j)
312        enddo        enddo
313        enddo        enddo
314        nalbvisdf = nalbvisdf + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) nalbvisdf = nalbvisdf + 1
315        endif        endif
316    
317        if( ialbnirdr.gt.0 ) then        if( ialbnirdr.gt.0 ) then
# Line 376  c ------------------ Line 321  c ------------------
321       .                                                     albnirdr(i,j)       .                                                     albnirdr(i,j)
322        enddo        enddo
323        enddo        enddo
324        nalbnirdr = nalbnirdr + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) nalbnirdr = nalbnirdr + 1
325        endif        endif
326    
327        if( ialbnirdf.gt.0 ) then        if( ialbnirdf.gt.0 ) then
# Line 386  c ------------------ Line 331  c ------------------
331       .                                                     albnirdf(i,j)       .                                                     albnirdf(i,j)
332        enddo        enddo
333        enddo        enddo
334        nalbnirdf = nalbnirdf + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) nalbnirdf = nalbnirdf + 1
335        endif        endif
336    
337  C Compute Optical Thicknesses and Diagnostics  C Compute Optical Thicknesses and Diagnostics
# Line 411  C -------------------------------------- Line 356  C --------------------------------------
356        enddo        enddo
357        enddo        enddo
358        enddo        enddo
359        ntauave = ntauave + 1        if ( (bi.eq.1) .and. (bj.eq.1) ) ntauave = ntauave + 1
360        endif        endif
361    
362        if( itaucld.gt.0 ) then        if( itaucld.gt.0 ) then
# Line 491  C ************************************** Line 436  C **************************************
436    
437        CALL STRIP ( zenith,COSZ,im*jm,ISTRIP,1,NN )        CALL STRIP ( zenith,COSZ,im*jm,ISTRIP,1,NN )
438    
439        CALL STRIP ( plze, ple   ,im*jm,ISTRIP,lm+1,NN)        CALL STRIP ( plze,  ple   ,im*jm,ISTRIP,lm+1,NN)
440        CALL STRIP ( pkz , pk    ,im*jm,ISTRIP,lm  ,NN)        CALL STRIP ( pkz ,  pk    ,im*jm,ISTRIP,lm  ,NN)
441        CALL STRIP ( tdry, tzl   ,im*jm,ISTRIP,lm  ,NN)        CALL STRIP ( dpres,dpstrip,im*jm,ISTRIP,lm  ,NN)
442        CALL STRIP ( qz  , qzl   ,im*jm,ISTRIP,lm  ,NN)        CALL STRIP ( tdry,  tzl   ,im*jm,ISTRIP,lm  ,NN)
443        CALL STRIP ( oz  , ozl   ,im*jm,ISTRIP,lm  ,NN)        CALL STRIP ( qz  ,  qzl   ,im*jm,ISTRIP,lm  ,NN)
444        CALL STRIP ( tau , taul  ,im*jm,ISTRIP,lm  ,NN)        CALL STRIP ( oz  ,  ozl   ,im*jm,ISTRIP,lm  ,NN)
445          CALL STRIP ( tau ,  taul  ,im*jm,ISTRIP,lm  ,NN)
446    
447        CALL STRIP ( albvisdr,albuvdr,im*jm,ISTRIP,1,NN )        CALL STRIP ( albvisdr,albuvdr,im*jm,ISTRIP,1,NN )
448        CALL STRIP ( albvisdf,albuvdf,im*jm,ISTRIP,1,NN )        CALL STRIP ( albvisdf,albuvdf,im*jm,ISTRIP,1,NN )
# Line 543  C ****     Compute Mass-Weighted Theta T Line 489  C ****     Compute Mass-Weighted Theta T
489  C **********************************************************************  C **********************************************************************
490    
491        do l=1,lm        do l=1,lm
       alf = grav/(cp*dsig(L)*100)  
492        do i=1,istrip        do i=1,istrip
493          alf = grav*(ple(i,Lm+1)-ptop)/(cp*dpstrip(i,L)*100)
494        dtsw (i,L) = alf*( flux   (i,L)-flux   (i,L+1) )/pk(i,L)        dtsw (i,L) = alf*( flux   (i,L)-flux   (i,L+1) )/pk(i,L)
495        dtswc(i,L) = alf*( fluxclr(i,L)-fluxclr(i,L+1) )/pk(i,L)        dtswc(i,L) = alf*( fluxclr(i,L)-fluxclr(i,L+1) )/pk(i,L)
496        enddo        enddo
# Line 638  C                  tau(im,jm,lm,2):  Sus Line 584  C                  tau(im,jm,lm,2):  Sus
584  C                  tau(im,jm,lm,3):  Raindrops  C                  tau(im,jm,lm,3):  Raindrops
585  C  C
586  C***********************************************************************  C***********************************************************************
 C*                  GODDARD LABORATORY FOR ATMOSPHERES                 *  
 C***********************************************************************  
587    
588        implicit none        implicit none
589    
590        integer  im,jm,lm,i,j,L        integer  im,jm,lm,i,j,L
591    
592        real  tl(im,jm,lm)        _RL  tl(im,jm,lm)
593        real  pl(im,jm,lm)        _RL  pl(im,jm,lm)
594        real ple(im,jm,lm+1)        _RL ple(im,jm,lm+1)
595        real  lz(im,jm,lm)        _RL  lz(im,jm,lm)
596        real  cf(im,jm,lm)        _RL  cf(im,jm,lm)
597        real cfm(im,jm,lm)        _RL cfm(im,jm,lm)
598        real tau(im,jm,lm,3)        _RL tau(im,jm,lm,3)
599        integer lwi(im,jm)        integer lwi(im,jm)
600    
601        real dp, alf, fracls, fraccu        _RL dp, alf, fracls, fraccu
602        real tauice, tauh2o, tauras        _RL tauice, tauh2o, tauras
603    
604  c Compute Cloud Optical Depths  c Compute Cloud Optical Depths
605  c ----------------------------  c ----------------------------
# Line 823  c*************************************** Line 767  c***************************************
767    
768  c-----Explicit Inline Directives  c-----Explicit Inline Directives
769    
770  #if CRAY  #ifdef CRAY
771  #if f77  #ifdef f77
772  cfpp$ expand (expmn)  cfpp$ expand (expmn)
773  #endif  #endif
774  #endif  #endif
775        real expmn        _RL expmn
776    
777  c-----input parameters  c-----input parameters
778    
779        integer m,n,ndim,np,ict,icb        integer m,n,ndim,np,ict,icb
780        real pl(m,ndim,np+1),ta(m,ndim,np),wa(m,ndim,np),oa(m,ndim,np)        _RL pl(m,ndim,np+1),ta(m,ndim,np),wa(m,ndim,np),oa(m,ndim,np)
781        real  taucld(m,ndim,np,2),reff(m,ndim,np,2)        _RL  taucld(m,ndim,np,2),reff(m,ndim,np,2)
782        real  fcld(m,ndim,np),taual(m,ndim,np)        _RL  fcld(m,ndim,np),taual(m,ndim,np)
783        real  rsirbm(m,ndim),rsirdf(m,ndim),        _RL  rsirbm(m,ndim),rsirdf(m,ndim),
784       *     rsuvbm(m,ndim),rsuvdf(m,ndim),cosz(m,ndim),co2       *     rsuvbm(m,ndim),rsuvdf(m,ndim),cosz(m,ndim),co2
785    
786  c-----output parameters  c-----output parameters
787    
788        real  flx(m,ndim,np+1),flc(m,ndim,np+1)        _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
789        real  fdirir(m,ndim),fdifir(m,ndim)        _RL  fdirir(m,ndim),fdifir(m,ndim)
790        real  fdirpar(m,ndim),fdifpar(m,ndim)        _RL  fdirpar(m,ndim),fdifpar(m,ndim)
791        real  fdiruv(m,ndim),fdifuv(m,ndim)        _RL  fdiruv(m,ndim),fdifuv(m,ndim)
792    
793  c-----temporary array  c-----temporary array
794    
795        integer i,j,k,ik        integer i,j,k
796        real  cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np)        _RL  cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np)
797        real  dp(m,n,np),wh(m,n,np),oh(m,n,np),scal(m,n,np)        _RL  dp(m,n,np),wh(m,n,np),oh(m,n,np),scal(m,n,np)
798        real  swh(m,n,np+1),so2(m,n,np+1),df(m,n,np+1)        _RL  swh(m,n,np+1),so2(m,n,np+1),df(m,n,np+1)
799        real  sdf(m,n),sclr(m,n),csm(m,n),taux,x        _RL  sdf(m,n),sclr(m,n),csm(m,n),x
800    
801  c-----------------------------------------------------------------  c-----------------------------------------------------------------
802    
# Line 1065  c*************************************** Line 1009  c***************************************
1009  c-----input parameters  c-----input parameters
1010    
1011        integer m,n,ndim,np,ict,icb        integer m,n,ndim,np,ict,icb
1012        real  cosz(m,ndim),fcld(m,ndim,np),taucld(m,ndim,np,2)        _RL  cosz(m,ndim),fcld(m,ndim,np),taucld(m,ndim,np,2)
1013    
1014  c-----output parameters  c-----output parameters
1015    
1016        real  cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np)        _RL  cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np)
1017    
1018  c-----temporary variables  c-----temporary variables
1019    
1020        integer i,j,k,im,it,ia,kk        integer i,j,k,im,it,ia,kk
1021        real   fm,ft,fa,xai,taucl,taux        _RL   fm,ft,fa,xai,taux
1022    
1023  c-----pre-computed table  c-----pre-computed table
1024    
1025        integer   nm,nt,na        integer   nm,nt,na
1026        parameter (nm=11,nt=9,na=11)        parameter (nm=11,nt=9,na=11)
1027        real   dm,dt,da,t1,caib(nm,nt,na),caif(nt,na)        _RL   dm,dt,da,t1,caib(nm,nt,na),caif(nt,na)
1028        parameter (dm=0.1,dt=0.30103,da=0.1,t1=-0.9031)        parameter (dm=0.1,dt=0.30103,da=0.1,t1=-0.9031)
1029    
1030  c-----include the pre-computed table for cai  c-----include the pre-computed table for cai
1031    
1032        include 'cai.dat'  #include "cai-dat.h"
1033        save caib,caif        save caib,caif
1034    
1035    
# Line 1283  c*************************************** Line 1227  c***************************************
1227    
1228  c-----Explicit Inline Directives  c-----Explicit Inline Directives
1229    
1230  #if CRAY  #ifdef CRAY
1231  #if f77  #ifdef f77
1232  cfpp$ expand (deledd)  cfpp$ expand (deledd)
1233  cfpp$ expand (sagpol)  cfpp$ expand (sagpol)
1234  cfpp$ expand (expmn)  cfpp$ expand (expmn)
1235  #endif  #endif
1236  #endif  #endif
1237        real expmn        _RL expmn
1238    
1239  c-----input parameters  c-----input parameters
1240    
1241        integer m,n,ndim,np,ict,icb        integer m,n,ndim,np,ict,icb
1242        real  taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np)        _RL  taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np)
1243        real  tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3)        _RL  tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3)
1244        real  rsirbm(m,ndim),rsirdf(m,ndim)        _RL  rsirbm(m,ndim),rsirdf(m,ndim)
1245        real  wh(m,n,np),taual(m,ndim,np),csm(m,n)        _RL  wh(m,n,np),taual(m,ndim,np),csm(m,n)
1246    
1247  c-----output (updated) parameters  c-----output (updated) parameters
1248    
1249        real  flx(m,ndim,np+1),flc(m,ndim,np+1)        _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
1250        real  fdirir(m,ndim),fdifir(m,ndim)        _RL  fdirir(m,ndim),fdifir(m,ndim)
1251    
1252  c-----static parameters  c-----static parameters
1253    
1254        integer nk,nband        integer nk,nband
1255        parameter (nk=10,nband=3)        parameter (nk=10,nband=3)
1256        real  xk(nk),hk(nband,nk),ssaal(nband),asyal(nband)        _RL  xk(nk),hk(nband,nk),ssaal(nband),asyal(nband)
1257        real  aia(nband,3),awa(nband,3),aig(nband,3),awg(nband,3)        _RL  aia(nband,3),awa(nband,3),aig(nband,3),awg(nband,3)
1258    
1259  c-----temporary array  c-----temporary array
1260    
1261        integer ib,ik,i,j,k        integer ib,ik,i,j,k
1262        real  ssacl(m,n,np),asycl(m,n,np)        _RL  ssacl(m,n,np),asycl(m,n,np)
1263        real  rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2),        _RL  rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2),
1264       *       rs(m,n,np+1,2),ts(m,n,np+1,2)       *       rs(m,n,np+1,2),ts(m,n,np+1,2)
1265        real  rssab(m,n,np+1),rabx(m,n,np+1),rsabx(m,n,np+1)        _RL  fall(m,n,np+1),fclr(m,n,np+1)
1266        real  fall(m,n,np+1),fclr(m,n,np+1)        _RL  fsdir(m,n),fsdif(m,n)
1267        real  fsdir(m,n),fsdif(m,n)  
1268          _RL  tauwv,tausto,ssatau,asysto,tauto,ssato,asyto
1269        real  tauwv,tausto,ssatau,asysto,tauto,ssato,asyto        _RL  taux,reff1,reff2,w1,w2,g1,g2
1270        real  taux,reff1,reff2,w1,w2,g1,g2        _RL  ssaclt(m,n),asyclt(m,n)
1271        real  ssaclt(m,n),asyclt(m,n)        _RL  rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)
1272        real  rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)        _RL  rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)
       real  rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)  
1273    
1274  c-----water vapor absorption coefficient for 10 k-intervals.  c-----water vapor absorption coefficient for 10 k-intervals.
1275  c     unit: cm^2/gm  c     unit: cm^2/gm
# Line 1459  c-----integration over the k-distributio Line 1402  c-----integration over the k-distributio
1402              do i= 1, m              do i= 1, m
1403    
1404               tauwv=xk(ik)*wh(i,j,k)               tauwv=xk(ik)*wh(i,j,k)
1405    
1406  c-----compute total optical thickness, single scattering albedo,  c-----compute total optical thickness, single scattering albedo,
1407  c     and asymmetry factor.  c     and asymmetry factor.
1408    
# Line 1585  c     in certain parallel processors. Line 1528  c     in certain parallel processors.
1528    
1529  c-----flux calculations  c-----flux calculations
1530    
1531           do k= 1, np+1
1532            do j= 1, n
1533             do i= 1, m
1534              fclr(i,j,k) = 0.
1535              fall(i,j,k) = 0.
1536             enddo
1537            enddo
1538           enddo
1539           do j= 1, n
1540            do i= 1, m
1541             fsdir(i,j) = 0.
1542             fsdif(i,j) = 0.
1543            enddo
1544           enddo
1545    
1546          call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,          call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,
1547       *               fclr,fall,fsdir,fsdif)       *               fclr,fall,fsdir,fsdif)
1548    
# Line 1693  c*************************************** Line 1651  c***************************************
1651    
1652  c-----Explicit Inline Directives    c-----Explicit Inline Directives  
1653        
1654  #if CRAY  #ifdef CRAY
1655  #if f77    #ifdef f77  
1656  cfpp$ expand (deledd)  cfpp$ expand (deledd)
1657  cfpp$ expand (sagpol)  cfpp$ expand (sagpol)
1658  #endif    #endif  
# Line 1703  cfpp$ expand (sagpol) Line 1661  cfpp$ expand (sagpol)
1661  c-----input parameters  c-----input parameters
1662    
1663        integer m,n,ndim,np,ict,icb        integer m,n,ndim,np,ict,icb
1664        real  taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np)        _RL  taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np)
1665        real  tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3)        _RL  tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3)
1666        real  oh(m,n,np),dp(m,n,np),taual(m,ndim,np)        _RL  oh(m,n,np),dp(m,n,np),taual(m,ndim,np)
1667        real  rsuvbm(m,ndim),rsuvdf(m,ndim),csm(m,n)        _RL  rsuvbm(m,ndim),rsuvdf(m,ndim),csm(m,n)
1668    
1669  c-----output (updated) parameter  c-----output (updated) parameter
1670    
1671        real  flx(m,ndim,np+1),flc(m,ndim,np+1)        _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
1672        real  fdirpar(m,ndim),fdifpar(m,ndim)        _RL  fdirpar(m,ndim),fdifpar(m,ndim)
1673        real  fdiruv(m,ndim),fdifuv(m,ndim)        _RL  fdiruv(m,ndim),fdifuv(m,ndim)
1674    
1675  c-----static parameters  c-----static parameters
1676    
1677        integer nband        integer nband
1678        parameter (nband=8)        parameter (nband=8)
1679        real  hk(nband),xk(nband),ry(nband)        _RL  hk(nband),xk(nband),ry(nband)
1680        real  asyal(nband),ssaal(nband),aig(3),awg(3)        _RL  asyal(nband),ssaal(nband),aig(3),awg(3)
1681    
1682  c-----temporary array  c-----temporary array
1683    
1684        integer i,j,k,ib        integer i,j,k,ib
1685        real  taurs,tauoz,tausto,ssatau,asysto,tauto,ssato,asyto        _RL  taurs,tauoz,tausto,ssatau,asysto,tauto,ssato,asyto
1686        real  taux,reff1,reff2,g1,g2,asycl(m,n,np)        _RL  taux,reff1,reff2,g1,g2,asycl(m,n,np)
1687        real  td(m,n,np+1,2),rr(m,n,np+1,2),tt(m,n,np+1,2),        _RL  td(m,n,np+1,2),rr(m,n,np+1,2),tt(m,n,np+1,2),
1688       *       rs(m,n,np+1,2),ts(m,n,np+1,2)       *       rs(m,n,np+1,2),ts(m,n,np+1,2)
1689        real  upflux(m,n,np+1),dwflux(m,n,np+1),        _RL  fall(m,n,np+1),fclr(m,n,np+1),fsdir(m,n),fsdif(m,n)
1690       *     rssab(m,n,np+1),rabx(m,n,np+1),rsabx(m,n,np+1)        _RL  asyclt(m,n)
1691        real  fall(m,n,np+1),fclr(m,n,np+1),fsdir(m,n),fsdif(m,n)        _RL  rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)
1692        real  asyclt(m,n)        _RL  rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)
       real  rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)  
       real  rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)  
1693    
1694  c-----hk is the fractional extra-terrestrial solar flux.  c-----hk is the fractional extra-terrestrial solar flux.
1695  c     the sum of hk is 0.47074.  c     the sum of hk is 0.47074.
# Line 1944  c-----compute reflectance and transmitta Line 1900  c-----compute reflectance and transmitta
1900    
1901  c-----flux calculations  c-----flux calculations
1902    
1903           do k= 1, np+1
1904            do j= 1, n
1905             do i= 1, m
1906              fclr(i,j,k) = 0.
1907              fall(i,j,k) = 0.
1908             enddo
1909            enddo
1910           enddo
1911           do j= 1, n
1912            do i= 1, m
1913             fsdir(i,j) = 0.
1914             fsdif(i,j) = 0.
1915            enddo
1916           enddo
1917          call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,          call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts,
1918       *               fclr,fall,fsdir,fsdif)       *               fclr,fall,fsdir,fsdif)
1919    
# Line 2010  c*************************************** Line 1980  c***************************************
1980    
1981  c-----Explicit Inline Directives    c-----Explicit Inline Directives  
1982        
1983  #if CRAY  #ifdef CRAY
1984  #if f77    #ifdef f77  
1985  cfpp$ expand (expmn)  cfpp$ expand (expmn)
1986  #endif    #endif  
1987  #endif  #endif
1988        real expmn        _RL expmn
1989    
1990        real  zero,one,two,three,four,fourth,seven,tumin        _RL  zero,one,two,three,four,fourth,seven,tumin
1991        parameter (one=1., three=3.)        parameter (one=1., three=3.)
1992        parameter (seven=7., two=2.)        parameter (seven=7., two=2.)
1993        parameter (four=4., fourth=.25)        parameter (four=4., fourth=.25)
1994        parameter (zero=0., tumin=1.e-20)        parameter (zero=0., tumin=1.e-20)
1995    
1996  c-----input parameters  c-----input parameters
1997        real  tau,ssc,g0,csm        _RL  tau,ssc,g0,csm
1998    
1999  c-----output parameters  c-----output parameters
2000        real  rr,tt,td        _RL  rr,tt,td
2001    
2002  c-----temporary parameters  c-----temporary parameters
2003    
2004        real  zth,ff,xx,taup,sscp,gp,gm1,gm2,gm3,akk,alf1,alf2,        _RL  zth,ff,xx,taup,sscp,gp,gm1,gm2,gm3,akk,alf1,alf2,
2005       *     all,bll,st7,st8,cll,dll,fll,ell,st1,st2,st3,st4       *     all1,bll,st7,st8,cll,dll,fll,ell,st1,st2,st3,st4
2006  c  c
2007                  zth = one / csm                  zth = one / csm
2008    
# Line 2069  c   alf1 and alf2 are alpha1 and alpha2 Line 2039  c   alf1 and alf2 are alpha1 and alpha2
2039                  alf1 = gm1 - gm3 * xx                  alf1 = gm1 - gm3 * xx
2040                  alf2 = gm2 + gm3 * xx                  alf2 = gm2 + gm3 * xx
2041    
2042  c  all is last term in eq(21) of K & H  c  all1 is last term in eq(21) of K & H
2043  c  bll is last term in eq(22) of K & H  c  bll is last term in eq(22) of K & H
2044    
2045                  xx  = akk * two                  xx  = akk * two
2046                  all = (gm3 - alf2 * zth    )*xx*td                  all1 = (gm3 - alf2 * zth    )*xx*td
2047                  bll = (one - gm3 + alf1*zth)*xx                  bll = (one - gm3 + alf1*zth)*xx
2048    
2049                  xx  = akk * zth                  xx  = akk * zth
# Line 2099  c  bll is last term in eq(22) of K & H Line 2069  c  bll is last term in eq(22) of K & H
2069  c   rr is r-hat of eq(21) of K & H  c   rr is r-hat of eq(21) of K & H
2070  c   tt is diffuse part of t-hat of eq(22) of K & H  c   tt is diffuse part of t-hat of eq(22) of K & H
2071    
2072                  rr =   ( cll-dll*st4    -all*st2)*st1                  rr =   ( cll-dll*st4    -all1*st2)*st1
2073                  tt = - ((fll-ell*st4)*td-bll*st2)*st1                  tt = - ((fll-ell*st4)*td-bll*st2)*st1
2074    
2075                  rr = max(rr,zero)                  rr = max(rr,zero)
# Line 2134  c*************************************** Line 2104  c***************************************
2104    
2105  c-----Explicit Inline Directives    c-----Explicit Inline Directives  
2106        
2107  #if CRAY  #ifdef CRAY
2108  #if f77    #ifdef f77  
2109  cfpp$ expand (expmn)  cfpp$ expand (expmn)
2110  #endif    #endif  
2111  #endif  #endif
2112        real expmn        _RL expmn
2113    
2114        real  one,three,four        _RL  one,three,four
2115        parameter (one=1., three=3., four=4.)        parameter (one=1., three=3., four=4.)
2116    
2117  c-----output parameters:  c-----output parameters:
2118    
2119        real  tau,ssc,g0        _RL  tau,ssc,g0
2120    
2121  c-----output parameters:  c-----output parameters:
2122    
2123        real  rll,tll        _RL  rll,tll
2124    
2125  c-----temporary arrays  c-----temporary arrays
2126    
2127        real  xx,uuu,ttt,emt,up1,um1,st1        _RL  xx,uuu,ttt,emt,up1,um1,st1
2128    
2129               xx  = one-ssc*g0               xx  = one-ssc*g0
2130               uuu = sqrt( xx/(one-ssc))               uuu = sqrt( xx/(one-ssc))
# Line 2176  c*************************************** Line 2146  c***************************************
2146    
2147  c*******************************************************************  c*******************************************************************
2148  c compute exponential for arguments in the range 0> fin > -10.  c compute exponential for arguments in the range 0> fin > -10.
2149    c*******************************************************************
2150          implicit none
2151          _RL  fin,expmn
2152    
2153          _RL one,expmin,e1,e2,e3,e4
2154        parameter (one=1.0, expmin=-10.0)        parameter (one=1.0, expmin=-10.0)
2155        parameter (e1=1.0,        e2=-2.507213e-1)        parameter (e1=1.0,        e2=-2.507213e-1)
2156        parameter (e3=2.92732e-2, e4=-3.827800e-3)        parameter (e3=2.92732e-2, e4=-3.827800e-3)
       real  fin,expmn  
2157    
2158        if (fin .lt. expmin) fin = expmin        if (fin .lt. expmin) fin = expmin
2159        expmn = ((e4*fin + e3)*fin+e2)*fin+e1        expmn = ((e4*fin + e3)*fin+e2)*fin+e1
# Line 2230  c-----input parameters Line 2203  c-----input parameters
2203    
2204        integer m,n,np,ict,icb        integer m,n,np,ict,icb
2205    
2206        real  rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2)        _RL  rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2)
2207        real  rs(m,n,np+1,2),ts(m,n,np+1,2)        _RL  rs(m,n,np+1,2),ts(m,n,np+1,2)
2208        real  cc(m,n,3)        _RL  cc(m,n,3)
2209    
2210  c-----temporary array  c-----temporary array
2211    
2212        integer i,j,k,ih,im,is        integer i,j,k,ih,im,is
2213        real  rra(m,n,np+1,2,2),tta(m,n,np+1,2,2),tda(m,n,np+1,2,2)        _RL  rra(m,n,np+1,2,2),tta(m,n,np+1,2,2),tda(m,n,np+1,2,2)
2214        real  rsa(m,n,np+1,2,2),rxa(m,n,np+1,2,2)        _RL  rsa(m,n,np+1,2,2),rxa(m,n,np+1,2,2)
2215        real  ch(m,n),cm(m,n),ct(m,n),flxdn(m,n,np+1)        _RL  ch(m,n),cm(m,n),ct(m,n),flxdn(m,n,np+1)
2216        real  fdndir(m,n),fdndif(m,n),fupdif        _RL  fdndir(m,n),fdndif(m,n),fupdif
2217        real  denm,xx        _RL  denm,xx
2218    
2219  c-----output parameters  c-----output parameters
2220    
2221        real  fclr(m,n,np+1),fall(m,n,np+1)        _RL  fclr(m,n,np+1),fall(m,n,np+1)
2222        real  fsdir(m,n),fsdif(m,n)        _RL  fsdir(m,n),fsdif(m,n)
2223    
2224  c-----initialize all-sky flux (fall) and surface downward fluxes  c-----initialize all-sky flux (fall) and surface downward fluxes
2225    
# Line 2546  c     due to co2 absorption. Line 2519  c     due to co2 absorption.
2519  c-----input parameters  c-----input parameters
2520    
2521        integer m,n,np        integer m,n,np
2522        real  csm(m,n),swc(m,n,np+1),swh(m,n,np+1),cah(22,19)        _RL  csm(m,n),swc(m,n,np+1),swh(m,n,np+1),cah(22,19)
2523    
2524  c-----output (undated) parameter  c-----output (undated) parameter
2525    
2526        real  df(m,n,np+1)        _RL  df(m,n,np+1)
2527    
2528  c-----temporary array  c-----temporary array
2529    
2530        integer i,j,k,ic,iw        integer i,j,k,ic,iw
2531        real  xx,clog,wlog,dc,dw,x1,x2,y2        _RL  xx,clog1,wlog,dc,dw,x1,x2,y2
2532    
2533  c********************************************************************  c********************************************************************
2534  c-----include co2 look-up table  c-----include co2 look-up table
2535    
2536        include 'cah.dat'  #include "cah-dat.h"
2537        save cah  c     save cah
2538    
2539  c********************************************************************  c********************************************************************
2540  c-----table look-up for the reduction of clear-sky solar  c-----table look-up for the reduction of clear-sky solar
# Line 2572  c     extraterrestrial solar flux in the Line 2545  c     extraterrestrial solar flux in the
2545         do j= 1, n         do j= 1, n
2546          do i= 1, m          do i= 1, m
2547            xx=1./.3            xx=1./.3
2548            clog=log10(swc(i,j,k)*csm(i,j))            clog1=log10(swc(i,j,k)*csm(i,j))
2549            wlog=log10(swh(i,j,k)*csm(i,j))            wlog=log10(swh(i,j,k)*csm(i,j))
2550            ic=int( (clog+3.15)*xx+1.)            ic=int( (clog1+3.15)*xx+1.)
2551            iw=int( (wlog+4.15)*xx+1.)            iw=int( (wlog+4.15)*xx+1.)
2552            if(ic.lt.2)ic=2            if(ic.lt.2)ic=2
2553            if(iw.lt.2)iw=2            if(iw.lt.2)iw=2
2554            if(ic.gt.22)ic=22            if(ic.gt.22)ic=22
2555            if(iw.gt.19)iw=19                if(iw.gt.19)iw=19    
2556            dc=clog-float(ic-2)*.3+3.            dc=clog1-float(ic-2)*.3+3.
2557            dw=wlog-float(iw-2)*.3+4.              dw=wlog-float(iw-2)*.3+4.  
2558            x1=cah(1,iw-1)+(cah(1,iw)-cah(1,iw-1))*xx*dw            x1=cah(1,iw-1)+(cah(1,iw)-cah(1,iw-1))*xx*dw
2559            x2=cah(ic-1,iw-1)+(cah(ic-1,iw)-cah(ic-1,iw-1))*xx*dw            x2=cah(ic-1,iw-1)+(cah(ic-1,iw)-cah(ic-1,iw-1))*xx*dw
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22