/[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.24 by ce107, Thu Jun 16 16:46:12 2005 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
 #ifdef ALLOW_DIAGNOSTICS  
 #include "diagnostics.h"  
 #endif  
15    
16  c Input Variables  c Input Variables
17  c ---------------  c ---------------
18        integer nymd,nhms,ndswr,istrip,npcs,bi,bj        integer nymd,nhms,bi,bj,ndswr,myid,istrip,npcs
19          integer mid_level,low_level
20        integer im,jm,lm                integer im,jm,lm        
21        real  ptop                      _RL  ptop
22        real  sige(lm+1)                _RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1),dpres(im,jm,lm)
23        real   sig(lm)                  _RL pkht(im,jm,lm+1),pkz(im,jm,lm)
24        real  dsig(lm)                  _RL tz(im,jm,lm),qz(im,jm,lm)
25          _RL oz(im,jm,lm)
26        real    pz(im,jm)              _RL co2
27        real    tz(im,jm,lm)            _RL albvisdr(im,jm),albvisdf(im,jm),albnirdr(im,jm)
28        real  pkht(im,jm,lm)            _RL albnirdf(im,jm)
29          _RL radswg(im,jm),swgclr(im,jm),fdifpar(im,jm),fdirpar(im,jm)
30        real    co2                    _RL osr(im,jm),osrclr(im,jm),dtradsw(im,jm,lm)
31        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)  
   
32        integer nswcld,nswlz            integer nswcld,nswlz    
33        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)    
   
34        logical lpnt                    logical lpnt            
35        integer imstturb                integer imstturb        
36        real qliqave(im,jm,lm)          _RL qliqave(im,jm,lm),fccave(im,jm,lm)  
       real  fccave(im,jm,lm)    
   
37        integer landtype(im,jm)        integer landtype(im,jm)
38          _RL xlats(im,jm),xlons(im,jm)
39    
40  c Local Variables  c Local Variables
41  c ---------------  c ---------------
42        integer   i,j,L,nn,nsecf,mid_level, low_level        integer   i,j,L,nn,nsecf
43        integer   nb2,ntmstp,nymd2,nhms2        integer   ntmstp,nymd2,nhms2
44        real      getcon,grav,cp,undef,pcheck        _RL      getcon,grav,cp,undef
45        real      ra,alf,reffw,reffi,tminv        _RL      ra,alf,reffw,reffi,tminv
46    
47        parameter ( reffw = 10.0 )          parameter ( reffw = 10.0 )  
48        parameter ( reffi = 65.0 )          parameter ( reffi = 65.0 )  
49    
50        real      alat(im,jm)        _RL tdry(im,jm,lm)
51        real      alon(im,jm)        _RL TEMP1(im,jm)
52          _RL TEMP2(im,jm)
53          _RL zenith (im,jm)
54          _RL cldtot (im,jm,lm)
55          _RL cldmxo (im,jm,lm)
56          _RL totcld (im,jm)
57          _RL cldlow (im,jm)
58          _RL cldmid (im,jm)
59          _RL cldhi  (im,jm)
60          _RL taulow (im,jm)
61          _RL taumid (im,jm)
62          _RL tauhi  (im,jm)
63          _RL tautype(im,jm,lm,3)
64          _RL tau(im,jm,lm)
65          _RL albedo(im,jm)    
66    
67          _RL PK(ISTRIP,lm)
68          _RL qzl(ISTRIP,lm),CLRO(ISTRIP,lm)
69          _RL TZL(ISTRIP,lm)
70          _RL OZL(ISTRIP,lm)
71          _RL PLE(ISTRIP,lm+1)
72          _RL COSZ(ISTRIP)
73          _RL dpstrip(ISTRIP,lm)
74    
75          _RL albuvdr(ISTRIP),albuvdf(ISTRIP)
76          _RL albirdr(ISTRIP),albirdf(ISTRIP)
77          _RL difpar (ISTRIP),dirpar (ISTRIP)
78    
79          _RL fdirir(istrip),fdifir(istrip)
80          _RL fdiruv(istrip),fdifuv(istrip)
81    
82          _RL flux(istrip,lm+1)
83          _RL fluxclr(istrip,lm+1)
84          _RL dtsw(istrip,lm)
85          _RL dtswc(istrip,lm)
86    
87          _RL taul   (istrip,lm)
88          _RL reff   (istrip,lm,2)
89          _RL tauc   (istrip,lm,2)
90          _RL taua   (istrip,lm)
91          _RL tstrip (istrip)
92    #ifdef ALLOW_DIAGNOSTICS
93          logical  diagnostics_is_on
94          external diagnostics_is_on
95          _RL tmpdiag(im,jm),tmpdiag2(im,jm)
96    #endif
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 224  c -------------------------------------- Line 169  c --------------------------------------
169          do L =1,lm          do L =1,lm
170          do j =1,jm          do j =1,jm
171          do i =1,im          do i =1,im
172           cldtot(i,j,L)=min(1.0,max(cldsw(i,j,L),fccave(i,j,L)/imstturb))           cldtot(i,j,L)=min(1.0 _d 0,max(cldsw(i,j,L),fccave(i,j,L)/imstturb))
173           cldmxo(i,j,L)=min(1.0,cswmo(i,j,L))           cldmxo(i,j,L)=min(1.0 _d 0,cswmo(i,j,L))
174             swlz(i,j,L)=swlz(i,j,L)+qliqave(i,j,L)/imstturb             swlz(i,j,L)=swlz(i,j,L)+qliqave(i,j,L)/imstturb
175          enddo          enddo
176          enddo          enddo
# Line 234  c -------------------------------------- Line 179  c --------------------------------------
179          do L =1,lm          do L =1,lm
180          do j =1,jm          do j =1,jm
181          do i =1,im          do i =1,im
182           cldtot(i,j,L) =  min( 1.0,cldsw(i,j,L) )           cldtot(i,j,L) =  min( 1.0 _d 0,cldsw(i,j,L) )
183           cldmxo(i,j,L) =  min( 1.0,cswmo(i,j,L) )           cldmxo(i,j,L) =  min( 1.0 _d 0,cswmo(i,j,L) )
184          enddo          enddo
185          enddo          enddo
186          enddo          enddo
# Line 273  c -------------------------------------- Line 218  c --------------------------------------
218        enddo        enddo
219        enddo        enddo
220    
221    #ifdef ALLOW_DIAGNOSTICS
222    
223  c Compute Cloud Diagnostics  c Compute Cloud Diagnostics
224  c -------------------------  c -------------------------
225        if(icldfrc.gt.0) then        if(diagnostics_is_on('CLDFRC  ',myid) ) then
226        do j=1,jm         call diagnostics_fill(totcld,'CLDFRC  ',0,1,3,bi,bj,myid)
       do i=1,im  
       qdiag(i,j,icldfrc,bi,bj) =  qdiag(i,j,icldfrc,bi,bj) + totcld(i,j)  
       enddo  
       enddo  
       ncldfrc = ncldfrc + 1  
227        endif        endif
228    
229        if( icldras.gt.0 ) then        if(diagnostics_is_on('CLDRAS  ',myid) ) then
230        do L=1,lm         do L=1,lm
231        do j=1,jm         do j=1,jm
232        do i=1,im         do i=1,im
233        qdiag(i,j,icldras+L-1,bi,bj) = qdiag(i,j,icldras+L-1,bi,bj) +          tmpdiag(i,j) = cswmo(i,j,L)
234       .                                                     cswmo(i,j,L)         enddo
235        enddo         enddo
236        enddo         call diagnostics_fill(tmpdiag,'CLDRAS  ',L,1,3,bi,bj,myid)
237        enddo        enddo
       ncldras = ncldras + 1  
238        endif        endif
239    
240        if( icldtot.gt.0 ) then        if(diagnostics_is_on('CLDTOT  ',myid) ) then
241        do L=1,lm         do L=1,lm
242        do j=1,jm         do j=1,jm
243        do i=1,im         do i=1,im
244        qdiag(i,j,icldtot+L-1,bi,bj) = qdiag(i,j,icldtot+L-1,bi,bj) +          tmpdiag(i,j) = cldtot(i,j,L)
245       .                                                     cldtot(i,j,L)         enddo
246        enddo         enddo
247        enddo         call diagnostics_fill(tmpdiag,'CLDTOT  ',L,1,3,bi,bj,myid)
248        enddo        enddo
       ncldtot = ncldtot + 1  
249        endif        endif
250    
251        if( icldlow.gt.0 ) then        if(diagnostics_is_on('CLDLOW  ',myid) ) then
252        do j=1,jm         call diagnostics_fill(cldlow,'CLDLOW  ',0,1,3,bi,bj,myid)
       do i=1,im  
       qdiag(i,j,icldlow,bi,bj) = qdiag(i,j,icldlow,bi,bj) + cldlow(i,j)  
       enddo  
       enddo  
       ncldlow = ncldlow + 1  
253        endif        endif
254    
255        if( icldmid.gt.0 ) then        if(diagnostics_is_on('CLDMID  ',myid) ) then
256        do j=1,jm         call diagnostics_fill(cldmid,'CLDMID  ',0,1,3,bi,bj,myid)
       do i=1,im  
       qdiag(i,j,icldmid,bi,bj) = qdiag(i,j,icldmid,bi,bj) + cldmid(i,j)  
       enddo  
       enddo  
       ncldmid = ncldmid + 1  
257        endif        endif
258    
259        if( icldhi.gt.0 ) then        if(diagnostics_is_on('CLDHI   ',myid) ) then
260        do j=1,jm         call diagnostics_fill(cldhi,'CLDHI   ',0,1,3,bi,bj,myid)
       do i=1,im  
       qdiag(i,j,icldhi,bi,bj) = qdiag(i,j,icldhi,bi,bj) + cldhi(i,j)  
       enddo  
       enddo  
       ncldhi = ncldhi + 1  
261        endif        endif
262    
263        if( ilzrad.gt.0 ) then        if(diagnostics_is_on('LZRAD   ',myid) ) then
264        do L=1,lm         do L=1,lm
265        do j=1,jm         do j=1,jm
266        do i=1,im         do i=1,im
267        qdiag(i,j,ilzrad+L-1,bi,bj) = qdiag(i,j,ilzrad+L-1,bi,bj) +          tmpdiag(i,j) = swlz(i,j,L) * 1.0e6
268       .                                                     swlz(i,j,L)*1.0e6         enddo
269        enddo         enddo
270        enddo         call diagnostics_fill(tmpdiag,'LZRAD   ',L,1,3,bi,bj,myid)
271        enddo         enddo
       nlzrad = nlzrad + 1  
272        endif        endif
273    
274  c Albedo Diagnostics  c Albedo Diagnostics
275  c ------------------  c ------------------
276        if( ialbvisdr.gt.0 ) then        if(diagnostics_is_on('ALBVISDR',myid) ) then
277        do j=1,jm         call diagnostics_fill(albvisdr,'ALBVISDR',0,1,3,bi,bj,myid)
       do i=1,im  
       qdiag(i,j,ialbvisdr,bi,bj) = qdiag(i,j,ialbvisdr,bi,bj) +  
      .                                                     albvisdr(i,j)  
       enddo  
       enddo  
       nalbvisdr = nalbvisdr + 1  
278        endif        endif
279    
280        if( ialbvisdf.gt.0 ) then        if(diagnostics_is_on('ALBVISDF',myid) ) then
281        do j=1,jm         call diagnostics_fill(albvisdf,'ALBVISDF',0,1,3,bi,bj,myid)
       do i=1,im  
       qdiag(i,j,ialbvisdf,bi,bj) = qdiag(i,j,ialbvisdf,bi,bj) +  
      .                                                     albvisdf(i,j)  
       enddo  
       enddo  
       nalbvisdf = nalbvisdf + 1  
282        endif        endif
283    
284        if( ialbnirdr.gt.0 ) then        if(diagnostics_is_on('ALBNIRDR',myid) ) then
285        do j=1,jm         call diagnostics_fill(albnirdr,'ALBNIRDR',0,1,3,bi,bj,myid)
       do i=1,im  
       qdiag(i,j,ialbnirdr,bi,bj) = qdiag(i,j,ialbnirdr,bi,bj) +  
      .                                                     albnirdr(i,j)  
       enddo  
       enddo  
       nalbnirdr = nalbnirdr + 1  
286        endif        endif
287    
288        if( ialbnirdf.gt.0 ) then        if(diagnostics_is_on('ALBNIRDF',myid) ) then
289        do j=1,jm         call diagnostics_fill(albnirdf,'ALBNIRDF',0,1,3,bi,bj,myid)
       do i=1,im  
       qdiag(i,j,ialbnirdf,bi,bj) = qdiag(i,j,ialbnirdf,bi,bj) +  
      .                                                     albnirdf(i,j)  
       enddo  
       enddo  
       nalbnirdf = nalbnirdf + 1  
290        endif        endif
291    
292    #endif
293    
294  C Compute Optical Thicknesses and Diagnostics  C Compute Optical Thicknesses and Diagnostics
295  C -------------------------------------------  C -------------------------------------------
296        call opthk(tdry,plz,plze,swlz,cldtot,cldmxo,landtype,im,jm,lm,        call opthk(tdry,plz,plze,swlz,cldtot,cldmxo,landtype,im,jm,lm,
# Line 402  C -------------------------------------- Line 304  C --------------------------------------
304        enddo        enddo
305        enddo        enddo
306    
307        if( itauave.gt.0 ) then  #ifdef ALLOW_DIAGNOSTICS
308        do L=1,lm        if(diagnostics_is_on('TAUAVE  ',myid) ) then
309        do j=1,jm         do L=1,lm
310        do i=1,im         do j=1,jm
311        qdiag(i,j,itauave+L-1,bi,bj) = qdiag(i,j,itauave+L-1,bi,bj) +         do i=1,im
312       .                        tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L))          tmpdiag(i,j) = tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L))
313        enddo         enddo
314        enddo         enddo
315           call diagnostics_fill(tmpdiag,'TAUAVE  ',L,1,3,bi,bj,myid)
316        enddo        enddo
       ntauave = ntauave + 1  
317        endif        endif
318    
319        if( itaucld.gt.0 ) then        if(diagnostics_is_on('TAUCLD  ',myid) .and.
320        do L=1,lm       .                 diagnostics_is_on('TAUCLDC ',myid) ) then
321        do j=1,jm         do L=1,lm
322        do i=1,im         do j=1,jm
323         if( cldtot(i,j,L).ne.0.0 ) then         do i=1,im
324          qdiag(i,j,itaucld +L-1,bi,bj) = qdiag(i,j,itaucld +L-1,bi,bj) +          if( cldtot(i,j,L).ne.0.0 ) then
325       .                        tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L))          tmpdiag(i,j) = tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L))
326          qdiag(i,j,itaucldc+L-1,bi,bj) =          tmpdiag2(i,j) = 1.
327       .                             qdiag(i,j,itaucldc+L-1,bi,bj) + 1.0          else
328         endif          tmpdiag(i,j) = 0.
329        enddo          tmpdiag2(i,j) = 0.
330        enddo          endif
331           enddo
332           enddo
333           call diagnostics_fill(tmpdiag,'TAUCLD  ',L,1,3,bi,bj,myid)
334           call diagnostics_fill(tmpdiag,'TAUCLDC ',L,1,3,bi,bj,myid)
335        enddo        enddo
336        endif        endif
337    
338  c Compute Low, Mid, and High Cloud Optical Depth Diagnostics  c Compute Low, Mid, and High Cloud Optical Depth Diagnostics
339  c ----------------------------------------------------------  c ----------------------------------------------------------
340        if( itaulow.ne.0 ) then        if(diagnostics_is_on('TAULOW  ',myid) .and.
341         .                 diagnostics_is_on('TAULOWC ',myid) ) then
342         do j = 1,jm         do j = 1,jm
343         do i = 1,im         do i = 1,im
344            taulow(i,j) =  0.0
345          if( cldlow(i,j).ne.0.0 ) then          if( cldlow(i,j).ne.0.0 ) then
          taulow(i,j) =  0.0  
346           do L = low_level,lm           do L = low_level,lm
347            taulow(i,j) = taulow(i,j) + tau(i,j,L)            taulow(i,j) = taulow(i,j) + tau(i,j,L)
348           enddo           enddo
349           qdiag(i,j,itaulow,bi,bj ) = qdiag(i,j,itaulow,bi,bj ) +           tmpdiag2(i,j) = 1.
350       .                                                    taulow(i,j)          else
351           qdiag(i,j,itaulowc,bi,bj) = qdiag(i,j,itaulowc,bi,bj) + 1.0           tmpdiag(i,j) = 0.
352          endif          endif
353         enddo         enddo
354         enddo         enddo
355           call diagnostics_fill(taulow,'TAULOW  ',0,1,3,bi,bj,myid)
356           call diagnostics_fill(tmpdiag2,'TAULOWC ',0,1,3,bi,bj,myid)
357        endif        endif
358    
359        if( itaumid.ne.0 ) then        if(diagnostics_is_on('TAUMID  ',myid) .and.
360         .                 diagnostics_is_on('TAUMIDC ',myid) ) then
361         do j = 1,jm         do j = 1,jm
362         do i = 1,im         do i = 1,im
363            taumid(i,j) =  0.0
364          if( cldmid(i,j).ne.0.0 ) then          if( cldmid(i,j).ne.0.0 ) then
          taumid(i,j) =  0.0  
365           do L = mid_level,low_level+1           do L = mid_level,low_level+1
366            taumid(i,j) = taumid(i,j) + tau(i,j,L)            taumid(i,j) = taumid(i,j) + tau(i,j,L)
367           enddo           enddo
368           qdiag(i,j,itaumid,bi,bj ) = qdiag(i,j,itaumid,bi,bj ) +           tmpdiag2(i,j) = 1.
369       .                                                    taumid(i,j)          else
370           qdiag(i,j,itaumidc,bi,bj) = qdiag(i,j,itaumidc,bi,bj) + 1.0           tmpdiag(i,j) = 0.
371          endif          endif
372         enddo         enddo
373         enddo         enddo
374           call diagnostics_fill(taumid,'TAUMID  ',0,1,3,bi,bj,myid)
375           call diagnostics_fill(tmpdiag2,'TAUMIDC ',0,1,3,bi,bj,myid)
376        endif        endif
377    
378        if( itauhi.ne.0 ) then        if(diagnostics_is_on('TAUHI   ',myid) .and.
379         .                 diagnostics_is_on('TAUHIC  ',myid) ) then
380         do j = 1,jm         do j = 1,jm
381         do i = 1,im         do i = 1,im
382            tauhi(i,j) =  0.0
383          if( cldhi(i,j).ne.0.0 ) then          if( cldhi(i,j).ne.0.0 ) then
          tauhi(i,j) =  0.0  
384           do L = 1,mid_level+1           do L = 1,mid_level+1
385            tauhi(i,j) = tauhi(i,j) + tau(i,j,L)            tauhi(i,j) = tauhi(i,j) + tau(i,j,L)
386           enddo           enddo
387           qdiag(i,j,itauhi,bi,bj ) = qdiag(i,j,itauhi,bi,bj ) +           tmpdiag2(i,j) = 1.
388       .                                                   tauhi(i,j)          else
389           qdiag(i,j,itauhic,bi,bj) = qdiag(i,j,itauhic,bi,bj) + 1.0           tmpdiag(i,j) = 0.
390          endif          endif
391         enddo         enddo
392         enddo         enddo
393           call diagnostics_fill(tauhi,'TAUHI   ',0,1,3,bi,bj,myid)
394           call diagnostics_fill(tmpdiag2,'TAUHIC  ',0,1,3,bi,bj,myid)
395        endif        endif
396    #endif
397    
398  C***********************************************************************  C***********************************************************************
399  C ****                     LOOP OVER GLOBAL REGIONS                 ****  C ****                     LOOP OVER GLOBAL REGIONS                 ****
# Line 491  C ************************************** Line 407  C **************************************
407    
408        CALL STRIP ( zenith,COSZ,im*jm,ISTRIP,1,NN )        CALL STRIP ( zenith,COSZ,im*jm,ISTRIP,1,NN )
409    
410        CALL STRIP ( plze, ple   ,im*jm,ISTRIP,lm+1,NN)        CALL STRIP ( plze,  ple   ,im*jm,ISTRIP,lm+1,NN)
411        CALL STRIP ( pkz , pk    ,im*jm,ISTRIP,lm  ,NN)        CALL STRIP ( pkz ,  pk    ,im*jm,ISTRIP,lm  ,NN)
412        CALL STRIP ( tdry, tzl   ,im*jm,ISTRIP,lm  ,NN)        CALL STRIP ( dpres,dpstrip,im*jm,ISTRIP,lm  ,NN)
413        CALL STRIP ( qz  , qzl   ,im*jm,ISTRIP,lm  ,NN)        CALL STRIP ( tdry,  tzl   ,im*jm,ISTRIP,lm  ,NN)
414        CALL STRIP ( oz  , ozl   ,im*jm,ISTRIP,lm  ,NN)        CALL STRIP ( qz  ,  qzl   ,im*jm,ISTRIP,lm  ,NN)
415        CALL STRIP ( tau , taul  ,im*jm,ISTRIP,lm  ,NN)        CALL STRIP ( oz  ,  ozl   ,im*jm,ISTRIP,lm  ,NN)
416          CALL STRIP ( tau ,  taul  ,im*jm,ISTRIP,lm  ,NN)
417    
418        CALL STRIP ( albvisdr,albuvdr,im*jm,ISTRIP,1,NN )        CALL STRIP ( albvisdr,albuvdr,im*jm,ISTRIP,1,NN )
419        CALL STRIP ( albvisdf,albuvdf,im*jm,ISTRIP,1,NN )        CALL STRIP ( albvisdf,albuvdf,im*jm,ISTRIP,1,NN )
# Line 509  c Partition Tau between Water and Ice Pa Line 426  c Partition Tau between Water and Ice Pa
426  c ---------------------------------------------  c ---------------------------------------------
427        do L= 1,lm        do L= 1,lm
428        do i= 1,istrip        do i= 1,istrip
429                alf = min( max((tzl(i,l)-253.15)/20.,0.) ,1.)                alf = min( max((tzl(i,l)-253.15)/20.,0. _d 0) ,1. _d 0)
430        taua(i,L)   = 0.        taua(i,L)   = 0.
431    
432        if( alf.ne.0.0 .and. alf.ne.1.0 ) then        if( alf.ne.0.0 .and. alf.ne.1.0 ) then
# Line 543  C ****     Compute Mass-Weighted Theta T Line 460  C ****     Compute Mass-Weighted Theta T
460  C **********************************************************************  C **********************************************************************
461    
462        do l=1,lm        do l=1,lm
       alf = grav/(cp*dsig(L)*100)  
463        do i=1,istrip        do i=1,istrip
464          alf = grav*(ple(i,Lm+1)-ptop)/(cp*dpstrip(i,L)*100)
465        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)
466        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)
467        enddo        enddo
# Line 577  c --------------------- Line 494  c ---------------------
494    
495   1000 CONTINUE   1000 CONTINUE
496    
497    #ifdef ALLOW_DIAGNOSTICS
498    
499  c Mean Albedo Diagnostic  c Mean Albedo Diagnostic
500  c ----------------------  c ----------------------
501        if( ialbedo.gt.0 ) then        if(diagnostics_is_on('ALBEDO  ',myid) .and.
502        do j=1,jm       .                 diagnostics_is_on('ALBEDOC ',myid) ) then
503        do i=1,im         do j=1,jm
504        if( albedo(i,j).ne.undef ) then         do i=1,im
505        qdiag(i,j,ialbedo,bi,bj ) = qdiag(i,j,ialbedo,bi,bj )+albedo(i,j)          if( albedo(i,j).ne.undef ) then
506        qdiag(i,j,ialbedoc,bi,bj) = qdiag(i,j,ialbedoc,bi,bj) + 1.0           tmpdiag(i,j) = albedo(i,j)
507        endif           tmpdiag2(i,j) = 1.
508        enddo          else
509        enddo           tmpdiag(i,j) = 0.
510             tmpdiag2(i,j) = 0.
511            endif
512           enddo
513           enddo
514           call diagnostics_fill(tmpdiag,'ALBEDO  ',0,1,3,bi,bj,myid)
515           call diagnostics_fill(tmpdiag2,'ALBEDOC ',0,1,3,bi,bj,myid)
516        endif        endif
517    #endif
518    
519  C **********************************************************************  C **********************************************************************
520  C ****                 ZERO-OUT OR BUMP COUNTERS                    ****  C ****                 ZERO-OUT OR BUMP COUNTERS                    ****
# Line 638  C                  tau(im,jm,lm,2):  Sus Line 564  C                  tau(im,jm,lm,2):  Sus
564  C                  tau(im,jm,lm,3):  Raindrops  C                  tau(im,jm,lm,3):  Raindrops
565  C  C
566  C***********************************************************************  C***********************************************************************
 C*                  GODDARD LABORATORY FOR ATMOSPHERES                 *  
 C***********************************************************************  
567    
568        implicit none        implicit none
569    
570        integer  im,jm,lm,i,j,L        integer  im,jm,lm,i,j,L
571    
572        real  tl(im,jm,lm)        _RL  tl(im,jm,lm)
573        real  pl(im,jm,lm)        _RL  pl(im,jm,lm)
574        real ple(im,jm,lm+1)        _RL ple(im,jm,lm+1)
575        real  lz(im,jm,lm)        _RL  lz(im,jm,lm)
576        real  cf(im,jm,lm)        _RL  cf(im,jm,lm)
577        real cfm(im,jm,lm)        _RL cfm(im,jm,lm)
578        real tau(im,jm,lm,3)        _RL tau(im,jm,lm,3)
579        integer lwi(im,jm)        integer lwi(im,jm)
580    
581        real dp, alf, fracls, fraccu        _RL dp, alf, fracls, fraccu
582        real tauice, tauh2o, tauras        _RL tauice, tauh2o, tauras
583    
584  c Compute Cloud Optical Depths  c Compute Cloud Optical Depths
585  c ----------------------------  c ----------------------------
586        do L=1,lm        do L=1,lm
587        do j=1,jm        do j=1,jm
588        do i=1,im        do i=1,im
589                   alf   =  min( max(( tl(i,j,L)-233.15)/20.,0.) ,1.)                   alf   =  min( max(( tl(i,j,L)-233.15)/20.,0. _d 0) ,1. _d 0)
590                    dp   =  ple(i,j,L+1)-ple(i,j,L)                    dp   =  ple(i,j,L+1)-ple(i,j,L)
591           tau(i,j,L,1)  = 0.0           tau(i,j,L,1)  = 0.0
592           tau(i,j,L,2)  = 0.0           tau(i,j,L,2)  = 0.0
# Line 682  c -------------------------------------- Line 606  c --------------------------------------
606    
607  c Large-Scale Ice  c Large-Scale Ice
608  c ---------------  c ---------------
609           tauice = max( 0.0002, 0.002*min(500*lz(i,j,L)*1000,1.0) )           tauice = max( 0.0002 _d 0, 0.002*min(500*lz(i,j,L)*1000,1.0 _d 0) )
610           tau(i,j,L,1) = fracls*(1-alf)*tauice*dp           tau(i,j,L,1) = fracls*(1-alf)*tauice*dp
611    
612  c Large-Scale Water  c Large-Scale Water
613  c -----------------  c -----------------
614  C Over Land  C Over Land
615           if( lwi(i,j).le.10 ) then           if( lwi(i,j).le.10 ) then
616            tauh2o = max( 0.0020, 0.200*min(200*lz(i,j,L)*1000,1.0) )              tauh2o = max( 0.0020 _d 0, 0.200*min(200*lz(i,j,L)*1000,1.0 _d 0) )  
617            tau(i,j,L,3) = fracls*alf*tauh2o*dp            tau(i,j,L,3) = fracls*alf*tauh2o*dp
618           else           else
619  C Non-Precipitation Clouds Over Ocean  C Non-Precipitation Clouds Over Ocean
# Line 698  C Non-Precipitation Clouds Over Ocean Line 622  C Non-Precipitation Clouds Over Ocean
622             tau(i,j,L,2) = fracls*alf*tauh2o*dp             tau(i,j,L,2) = fracls*alf*tauh2o*dp
623            else            else
624  C Over Ocean  C Over Ocean
625             tauh2o = max( 0.0020, 0.120*min( 20*lz(i,j,L)*1000,1.0) )               tauh2o = max( 0.0020 _d 0, 0.120*min( 20*lz(i,j,L)*1000,1.0 _d 0) )  
626             tau(i,j,L,3) = fracls*alf*tauh2o*dp             tau(i,j,L,3) = fracls*alf*tauh2o*dp
627            endif            endif
628           endif           endif
# Line 823  c*************************************** Line 747  c***************************************
747    
748  c-----Explicit Inline Directives  c-----Explicit Inline Directives
749    
750  #if CRAY  #ifdef CRAY
751  #if f77  #ifdef f77
752  cfpp$ expand (expmn)  cfpp$ expand (expmn)
753  #endif  #endif
754  #endif  #endif
755        real expmn        _RL expmn
756    
757  c-----input parameters  c-----input parameters
758    
759        integer m,n,ndim,np,ict,icb        integer m,n,ndim,np,ict,icb
760        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)
761        real  taucld(m,ndim,np,2),reff(m,ndim,np,2)        _RL  taucld(m,ndim,np,2),reff(m,ndim,np,2)
762        real  fcld(m,ndim,np),taual(m,ndim,np)        _RL  fcld(m,ndim,np),taual(m,ndim,np)
763        real  rsirbm(m,ndim),rsirdf(m,ndim),        _RL  rsirbm(m,ndim),rsirdf(m,ndim),
764       *     rsuvbm(m,ndim),rsuvdf(m,ndim),cosz(m,ndim),co2       *     rsuvbm(m,ndim),rsuvdf(m,ndim),cosz(m,ndim),co2
765    
766  c-----output parameters  c-----output parameters
767    
768        real  flx(m,ndim,np+1),flc(m,ndim,np+1)        _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
769        real  fdirir(m,ndim),fdifir(m,ndim)        _RL  fdirir(m,ndim),fdifir(m,ndim)
770        real  fdirpar(m,ndim),fdifpar(m,ndim)        _RL  fdirpar(m,ndim),fdifpar(m,ndim)
771        real  fdiruv(m,ndim),fdifuv(m,ndim)        _RL  fdiruv(m,ndim),fdifuv(m,ndim)
772    
773  c-----temporary array  c-----temporary array
774    
775        integer i,j,k,ik        integer i,j,k
776        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)
777        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)
778        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)
779        real  sdf(m,n),sclr(m,n),csm(m,n),taux,x        _RL  sdf(m,n),sclr(m,n),csm(m,n),x
780    
781  c-----------------------------------------------------------------  c-----------------------------------------------------------------
782    
# Line 1065  c*************************************** Line 989  c***************************************
989  c-----input parameters  c-----input parameters
990    
991        integer m,n,ndim,np,ict,icb        integer m,n,ndim,np,ict,icb
992        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)
993    
994  c-----output parameters  c-----output parameters
995    
996        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)
997    
998  c-----temporary variables  c-----temporary variables
999    
1000        integer i,j,k,im,it,ia,kk        integer i,j,k,im,it,ia,kk
1001        real   fm,ft,fa,xai,taucl,taux        _RL   fm,ft,fa,xai,taux
1002    
1003  c-----pre-computed table  c-----pre-computed table
1004    
1005        integer   nm,nt,na        integer   nm,nt,na
1006        parameter (nm=11,nt=9,na=11)        parameter (nm=11,nt=9,na=11)
1007        real   dm,dt,da,t1,caib(nm,nt,na),caif(nt,na)        _RL   dm,dt,da,t1,caib(nm,nt,na),caif(nt,na)
1008        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)
1009    
1010  c-----include the pre-computed table for cai  c-----include the pre-computed table for cai
1011    
1012        include 'cai.dat'  #include "cai-dat.h"
1013        save caib,caif        save caib,caif
1014    
1015    
# Line 1155  c-----normalize cloud cover Line 1079  c-----normalize cloud cover
1079    
1080  c-----table look-up  c-----table look-up
1081    
1082             taux=min(taux,32.)             taux=min(taux,32. _d 0)
1083    
1084             fm=cosz(i,j)/dm             fm=cosz(i,j)/dm
1085             ft=(log10(taux)-t1)/dt             ft=(log10(taux)-t1)/dt
# Line 1191  c     angle, optical thickness, and clou Line 1115  c     angle, optical thickness, and clou
1115       *     caib(im,it,ia+1)*(1.+fa))*fa*.5+caib(im,it,ia)*(1.-fa*fa)       *     caib(im,it,ia+1)*(1.+fa))*fa*.5+caib(im,it,ia)*(1.-fa*fa)
1116    
1117             xai= xai-2.*caib(im,it,ia)             xai= xai-2.*caib(im,it,ia)
1118             xai=max(xai,0.0)             xai=max(xai,0.0 _d 0)
1119            
1120             tauclb(i,j,k) = taux*xai             tauclb(i,j,k) = taux*xai
1121    
# Line 1206  c     thickness and cover but not the so Line 1130  c     thickness and cover but not the so
1130       *      caif(it,ia+1)*(1.+fa))*fa*.5+caif(it,ia)*(1.-fa*fa)       *      caif(it,ia+1)*(1.+fa))*fa*.5+caif(it,ia)*(1.-fa*fa)
1131    
1132             xai= xai-caif(it,ia)             xai= xai-caif(it,ia)
1133             xai=max(xai,0.0)             xai=max(xai,0.0 _d 0)
1134            
1135             tauclf(i,j,k) = taux*xai             tauclf(i,j,k) = taux*xai
1136    
# Line 1283  c*************************************** Line 1207  c***************************************
1207    
1208  c-----Explicit Inline Directives  c-----Explicit Inline Directives
1209    
1210  #if CRAY  #ifdef CRAY
1211  #if f77  #ifdef f77
1212  cfpp$ expand (deledd)  cfpp$ expand (deledd)
1213  cfpp$ expand (sagpol)  cfpp$ expand (sagpol)
1214  cfpp$ expand (expmn)  cfpp$ expand (expmn)
1215  #endif  #endif
1216  #endif  #endif
1217        real expmn        _RL expmn
1218    
1219  c-----input parameters  c-----input parameters
1220    
1221        integer m,n,ndim,np,ict,icb        integer m,n,ndim,np,ict,icb
1222        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)
1223        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)
1224        real  rsirbm(m,ndim),rsirdf(m,ndim)        _RL  rsirbm(m,ndim),rsirdf(m,ndim)
1225        real  wh(m,n,np),taual(m,ndim,np),csm(m,n)        _RL  wh(m,n,np),taual(m,ndim,np),csm(m,n)
1226    
1227  c-----output (updated) parameters  c-----output (updated) parameters
1228    
1229        real  flx(m,ndim,np+1),flc(m,ndim,np+1)        _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
1230        real  fdirir(m,ndim),fdifir(m,ndim)        _RL  fdirir(m,ndim),fdifir(m,ndim)
1231    
1232  c-----static parameters  c-----static parameters
1233    
1234        integer nk,nband        integer nk,nband
1235        parameter (nk=10,nband=3)        parameter (nk=10,nband=3)
1236        real  xk(nk),hk(nband,nk),ssaal(nband),asyal(nband)        _RL  xk(nk),hk(nband,nk),ssaal(nband),asyal(nband)
1237        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)
1238    
1239  c-----temporary array  c-----temporary array
1240    
1241        integer ib,ik,i,j,k        integer ib,ik,i,j,k
1242        real  ssacl(m,n,np),asycl(m,n,np)        _RL  ssacl(m,n,np),asycl(m,n,np)
1243        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),
1244       *       rs(m,n,np+1,2),ts(m,n,np+1,2)       *       rs(m,n,np+1,2),ts(m,n,np+1,2)
1245        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)
1246        real  fall(m,n,np+1),fclr(m,n,np+1)        _RL  fsdir(m,n),fsdif(m,n)
1247        real  fsdir(m,n),fsdif(m,n)  
1248          _RL  tauwv,tausto,ssatau,asysto,tauto,ssato,asyto
1249        real  tauwv,tausto,ssatau,asysto,tauto,ssato,asyto        _RL  taux,reff1,reff2,w1,w2,g1,g2
1250        real  taux,reff1,reff2,w1,w2,g1,g2        _RL  ssaclt(m,n),asyclt(m,n)
1251        real  ssaclt(m,n),asyclt(m,n)        _RL  rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)
1252        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)  
1253    
1254  c-----water vapor absorption coefficient for 10 k-intervals.  c-----water vapor absorption coefficient for 10 k-intervals.
1255  c     unit: cm^2/gm  c     unit: cm^2/gm
# Line 1418  c     for a mixture of ice and liquid pa Line 1341  c     for a mixture of ice and liquid pa
1341             taux=taucld(i,j,k,1)+taucld(i,j,k,2)             taux=taucld(i,j,k,1)+taucld(i,j,k,2)
1342            if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then            if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then
1343    
1344             reff1=min(reff(i,j,k,1),130.)             reff1=min(reff(i,j,k,1),130. _d 0)
1345             reff2=min(reff(i,j,k,2),20.0)             reff2=min(reff(i,j,k,2),20.0 _d 0)
1346    
1347             w1=(1.-(aia(ib,1)+(aia(ib,2)+             w1=(1.-(aia(ib,1)+(aia(ib,2)+
1348       *         aia(ib,3)*reff1)*reff1))*taucld(i,j,k,1)       *         aia(ib,3)*reff1)*reff1))*taucld(i,j,k,1)
# Line 1459  c-----integration over the k-distributio Line 1382  c-----integration over the k-distributio
1382              do i= 1, m              do i= 1, m
1383    
1384               tauwv=xk(ik)*wh(i,j,k)               tauwv=xk(ik)*wh(i,j,k)
1385    
1386  c-----compute total optical thickness, single scattering albedo,  c-----compute total optical thickness, single scattering albedo,
1387  c     and asymmetry factor.  c     and asymmetry factor.
1388    
# Line 1474  c-----compute reflectance and transmitta Line 1397  c-----compute reflectance and transmitta
1397    
1398  c            if (ssato .gt. 0.001) then  c            if (ssato .gt. 0.001) then
1399    
1400  c             ssato=min(ssato,0.999999)  c             ssato=min(ssato,0.999999 _d 0)
1401  c             asyto=asysto/(ssato*tauto)  c             asyto=asysto/(ssato*tauto)
1402    
1403  c             call deledd(tauto,ssato,asyto,csm(i,j),  c             call deledd(tauto,ssato,asyto,csm(i,j),
# Line 1506  c-----compute reflectance and transmitta Line 1429  c-----compute reflectance and transmitta
1429    
1430                tauto=tausto+tauclb(i,j,k)                tauto=tausto+tauclb(i,j,k)
1431                ssato=(ssatau+ssacl(i,j,k)*tauclb(i,j,k))/tauto+1.0e-8                ssato=(ssatau+ssacl(i,j,k)*tauclb(i,j,k))/tauto+1.0e-8
1432                ssato=min(ssato,0.999999)                ssato=min(ssato,0.999999 _d 0)
1433                asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclb(i,j,k))/                asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclb(i,j,k))/
1434       *              (ssato*tauto)       *              (ssato*tauto)
1435    
# Line 1515  c-----compute reflectance and transmitta Line 1438  c-----compute reflectance and transmitta
1438    
1439                tauto=tausto+tauclf(i,j,k)                tauto=tausto+tauclf(i,j,k)
1440                ssato=(ssatau+ssacl(i,j,k)*tauclf(i,j,k))/tauto+1.0e-8                ssato=(ssatau+ssacl(i,j,k)*tauclf(i,j,k))/tauto+1.0e-8
1441                ssato=min(ssato,0.999999)                ssato=min(ssato,0.999999 _d 0)
1442                asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclf(i,j,k))/                asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclf(i,j,k))/
1443       *              (ssato*tauto)       *              (ssato*tauto)
1444    
# Line 1585  c     in certain parallel processors. Line 1508  c     in certain parallel processors.
1508    
1509  c-----flux calculations  c-----flux calculations
1510    
1511           do k= 1, np+1
1512            do j= 1, n
1513             do i= 1, m
1514              fclr(i,j,k) = 0.
1515              fall(i,j,k) = 0.
1516             enddo
1517            enddo
1518           enddo
1519           do j= 1, n
1520            do i= 1, m
1521             fsdir(i,j) = 0.
1522             fsdif(i,j) = 0.
1523            enddo
1524           enddo
1525    
1526          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,
1527       *               fclr,fall,fsdir,fsdif)       *               fclr,fall,fsdir,fsdif)
1528    
# Line 1693  c*************************************** Line 1631  c***************************************
1631    
1632  c-----Explicit Inline Directives    c-----Explicit Inline Directives  
1633        
1634  #if CRAY  #ifdef CRAY
1635  #if f77    #ifdef f77  
1636  cfpp$ expand (deledd)  cfpp$ expand (deledd)
1637  cfpp$ expand (sagpol)  cfpp$ expand (sagpol)
1638  #endif    #endif  
# Line 1703  cfpp$ expand (sagpol) Line 1641  cfpp$ expand (sagpol)
1641  c-----input parameters  c-----input parameters
1642    
1643        integer m,n,ndim,np,ict,icb        integer m,n,ndim,np,ict,icb
1644        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)
1645        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)
1646        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)
1647        real  rsuvbm(m,ndim),rsuvdf(m,ndim),csm(m,n)        _RL  rsuvbm(m,ndim),rsuvdf(m,ndim),csm(m,n)
1648    
1649  c-----output (updated) parameter  c-----output (updated) parameter
1650    
1651        real  flx(m,ndim,np+1),flc(m,ndim,np+1)        _RL  flx(m,ndim,np+1),flc(m,ndim,np+1)
1652        real  fdirpar(m,ndim),fdifpar(m,ndim)        _RL  fdirpar(m,ndim),fdifpar(m,ndim)
1653        real  fdiruv(m,ndim),fdifuv(m,ndim)        _RL  fdiruv(m,ndim),fdifuv(m,ndim)
1654    
1655  c-----static parameters  c-----static parameters
1656    
1657        integer nband        integer nband
1658        parameter (nband=8)        parameter (nband=8)
1659        real  hk(nband),xk(nband),ry(nband)        _RL  hk(nband),xk(nband),ry(nband)
1660        real  asyal(nband),ssaal(nband),aig(3),awg(3)        _RL  asyal(nband),ssaal(nband),aig(3),awg(3)
1661    
1662  c-----temporary array  c-----temporary array
1663    
1664        integer i,j,k,ib        integer i,j,k,ib
1665        real  taurs,tauoz,tausto,ssatau,asysto,tauto,ssato,asyto        _RL  taurs,tauoz,tausto,ssatau,asysto,tauto,ssato,asyto
1666        real  taux,reff1,reff2,g1,g2,asycl(m,n,np)        _RL  taux,reff1,reff2,g1,g2,asycl(m,n,np)
1667        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),
1668       *       rs(m,n,np+1,2),ts(m,n,np+1,2)       *       rs(m,n,np+1,2),ts(m,n,np+1,2)
1669        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)
1670       *     rssab(m,n,np+1),rabx(m,n,np+1),rsabx(m,n,np+1)        _RL  asyclt(m,n)
1671        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)
1672        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)  
1673    
1674  c-----hk is the fractional extra-terrestrial solar flux.  c-----hk is the fractional extra-terrestrial solar flux.
1675  c     the sum of hk is 0.47074.  c     the sum of hk is 0.47074.
# Line 1796  c     liquid and ice particles.  unit of Line 1732  c     liquid and ice particles.  unit of
1732             taux=taucld(i,j,k,1)+taucld(i,j,k,2)             taux=taucld(i,j,k,1)+taucld(i,j,k,2)
1733            if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then            if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then
1734    
1735             reff1=min(reff(i,j,k,1),130.)             reff1=min(reff(i,j,k,1),130. _d 0)
1736             reff2=min(reff(i,j,k,2),20.0)             reff2=min(reff(i,j,k,2),20.0 _d 0)
1737    
1738             g1=(aig(1)+(aig(2)+aig(3)*reff1)*reff1)*taucld(i,j,k,1)             g1=(aig(1)+(aig(2)+aig(3)*reff1)*reff1)*taucld(i,j,k,1)
1739             g2=(awg(1)+(awg(2)+awg(3)*reff2)*reff2)*taucld(i,j,k,2)             g2=(awg(1)+(awg(2)+awg(3)*reff2)*reff2)*taucld(i,j,k,2)
# Line 1848  c-----compute reflectance and transmitta Line 1784  c-----compute reflectance and transmitta
1784    
1785            tauto=tausto            tauto=tausto
1786            ssato=ssatau/tauto+1.0e-8            ssato=ssatau/tauto+1.0e-8
1787            ssato=min(ssato,0.999999)            ssato=min(ssato,0.999999 _d 0)
1788            asyto=asysto/(ssato*tauto)            asyto=asysto/(ssato*tauto)
1789    
1790            call deledd(tauto,ssato,asyto,csm(i,j),            call deledd(tauto,ssato,asyto,csm(i,j),
# Line 1870  c-----compute reflectance and transmitta Line 1806  c-----compute reflectance and transmitta
1806    
1807             tauto=tausto+tauclb(i,j,k)             tauto=tausto+tauclb(i,j,k)
1808             ssato=(ssatau+tauclb(i,j,k))/tauto+1.0e-8             ssato=(ssatau+tauclb(i,j,k))/tauto+1.0e-8
1809             ssato=min(ssato,0.999999)             ssato=min(ssato,0.999999 _d 0)
1810             asyto=(asysto+asycl(i,j,k)*tauclb(i,j,k))/(ssato*tauto)             asyto=(asysto+asycl(i,j,k)*tauclb(i,j,k))/(ssato*tauto)
1811    
1812             call deledd(tauto,ssato,asyto,csm(i,j),             call deledd(tauto,ssato,asyto,csm(i,j),
# Line 1878  c-----compute reflectance and transmitta Line 1814  c-----compute reflectance and transmitta
1814    
1815             tauto=tausto+tauclf(i,j,k)             tauto=tausto+tauclf(i,j,k)
1816             ssato=(ssatau+tauclf(i,j,k))/tauto+1.0e-8             ssato=(ssatau+tauclf(i,j,k))/tauto+1.0e-8
1817             ssato=min(ssato,0.999999)             ssato=min(ssato,0.999999 _d 0)
1818             asyto=(asysto+asycl(i,j,k)*tauclf(i,j,k))/(ssato*tauto)             asyto=(asysto+asycl(i,j,k)*tauclf(i,j,k))/(ssato*tauto)
1819    
1820             call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j))             call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j))
# Line 1944  c-----compute reflectance and transmitta Line 1880  c-----compute reflectance and transmitta
1880    
1881  c-----flux calculations  c-----flux calculations
1882    
1883           do k= 1, np+1
1884            do j= 1, n
1885             do i= 1, m
1886              fclr(i,j,k) = 0.
1887              fall(i,j,k) = 0.
1888             enddo
1889            enddo
1890           enddo
1891           do j= 1, n
1892            do i= 1, m
1893             fsdir(i,j) = 0.
1894             fsdif(i,j) = 0.
1895            enddo
1896           enddo
1897          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,
1898       *               fclr,fall,fsdir,fsdif)       *               fclr,fall,fsdir,fsdif)
1899    
# Line 2010  c*************************************** Line 1960  c***************************************
1960    
1961  c-----Explicit Inline Directives    c-----Explicit Inline Directives  
1962        
1963  #if CRAY  #ifdef CRAY
1964  #if f77    #ifdef f77  
1965  cfpp$ expand (expmn)  cfpp$ expand (expmn)
1966  #endif    #endif  
1967  #endif  #endif
1968        real expmn        _RL expmn
1969    
1970        real  zero,one,two,three,four,fourth,seven,tumin        _RL  zero,one,two,three,four,fourth,seven,tumin
1971        parameter (one=1., three=3.)        parameter (one=1., three=3.)
1972        parameter (seven=7., two=2.)        parameter (seven=7., two=2.)
1973        parameter (four=4., fourth=.25)        parameter (four=4., fourth=.25)
1974        parameter (zero=0., tumin=1.e-20)        parameter (zero=0., tumin=1.e-20)
1975    
1976  c-----input parameters  c-----input parameters
1977        real  tau,ssc,g0,csm        _RL  tau,ssc,g0,csm
1978    
1979  c-----output parameters  c-----output parameters
1980        real  rr,tt,td        _RL  rr,tt,td
1981    
1982  c-----temporary parameters  c-----temporary parameters
1983    
1984        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,
1985       *     all,bll,st7,st8,cll,dll,fll,ell,st1,st2,st3,st4       *     all1,bll,st7,st8,cll,dll,fll,ell,st1,st2,st3,st4
1986  c  c
1987                  zth = one / csm                  zth = one / csm
1988    
# Line 2069  c   alf1 and alf2 are alpha1 and alpha2 Line 2019  c   alf1 and alf2 are alpha1 and alpha2
2019                  alf1 = gm1 - gm3 * xx                  alf1 = gm1 - gm3 * xx
2020                  alf2 = gm2 + gm3 * xx                  alf2 = gm2 + gm3 * xx
2021    
2022  c  all is last term in eq(21) of K & H  c  all1 is last term in eq(21) of K & H
2023  c  bll is last term in eq(22) of K & H  c  bll is last term in eq(22) of K & H
2024    
2025                  xx  = akk * two                  xx  = akk * two
2026                  all = (gm3 - alf2 * zth    )*xx*td                  all1 = (gm3 - alf2 * zth    )*xx*td
2027                  bll = (one - gm3 + alf1*zth)*xx                  bll = (one - gm3 + alf1*zth)*xx
2028    
2029                  xx  = akk * zth                  xx  = akk * zth
# Line 2099  c  bll is last term in eq(22) of K & H Line 2049  c  bll is last term in eq(22) of K & H
2049  c   rr is r-hat of eq(21) of K & H  c   rr is r-hat of eq(21) of K & H
2050  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
2051    
2052                  rr =   ( cll-dll*st4    -all*st2)*st1                  rr =   ( cll-dll*st4    -all1*st2)*st1
2053                  tt = - ((fll-ell*st4)*td-bll*st2)*st1                  tt = - ((fll-ell*st4)*td-bll*st2)*st1
2054    
2055                  rr = max(rr,zero)                  rr = max(rr,zero)
# Line 2134  c*************************************** Line 2084  c***************************************
2084    
2085  c-----Explicit Inline Directives    c-----Explicit Inline Directives  
2086        
2087  #if CRAY  #ifdef CRAY
2088  #if f77    #ifdef f77  
2089  cfpp$ expand (expmn)  cfpp$ expand (expmn)
2090  #endif    #endif  
2091  #endif  #endif
2092        real expmn        _RL expmn
2093    
2094        real  one,three,four        _RL  one,three,four
2095        parameter (one=1., three=3., four=4.)        parameter (one=1., three=3., four=4.)
2096    
2097  c-----output parameters:  c-----output parameters:
2098    
2099        real  tau,ssc,g0        _RL  tau,ssc,g0
2100    
2101  c-----output parameters:  c-----output parameters:
2102    
2103        real  rll,tll        _RL  rll,tll
2104    
2105  c-----temporary arrays  c-----temporary arrays
2106    
2107        real  xx,uuu,ttt,emt,up1,um1,st1        _RL  xx,uuu,ttt,emt,up1,um1,st1
2108    
2109               xx  = one-ssc*g0               xx  = one-ssc*g0
2110               uuu = sqrt( xx/(one-ssc))               uuu = sqrt( xx/(one-ssc))
# Line 2176  c*************************************** Line 2126  c***************************************
2126    
2127  c*******************************************************************  c*******************************************************************
2128  c compute exponential for arguments in the range 0> fin > -10.  c compute exponential for arguments in the range 0> fin > -10.
2129    c*******************************************************************
2130          implicit none
2131          _RL  fin,expmn
2132    
2133          _RL one,expmin,e1,e2,e3,e4
2134        parameter (one=1.0, expmin=-10.0)        parameter (one=1.0, expmin=-10.0)
2135        parameter (e1=1.0,        e2=-2.507213e-1)        parameter (e1=1.0,        e2=-2.507213e-1)
2136        parameter (e3=2.92732e-2, e4=-3.827800e-3)        parameter (e3=2.92732e-2, e4=-3.827800e-3)
       real  fin,expmn  
2137    
2138        if (fin .lt. expmin) fin = expmin        if (fin .lt. expmin) fin = expmin
2139        expmn = ((e4*fin + e3)*fin+e2)*fin+e1        expmn = ((e4*fin + e3)*fin+e2)*fin+e1
# Line 2230  c-----input parameters Line 2183  c-----input parameters
2183    
2184        integer m,n,np,ict,icb        integer m,n,np,ict,icb
2185    
2186        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)
2187        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)
2188        real  cc(m,n,3)        _RL  cc(m,n,3)
2189    
2190  c-----temporary array  c-----temporary array
2191    
2192        integer i,j,k,ih,im,is        integer i,j,k,ih,im,is
2193        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)
2194        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)
2195        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)
2196        real  fdndir(m,n),fdndif(m,n),fupdif        _RL  fdndir(m,n),fdndif(m,n),fupdif
2197        real  denm,xx        _RL  denm,xx
2198    
2199  c-----output parameters  c-----output parameters
2200    
2201        real  fclr(m,n,np+1),fall(m,n,np+1)        _RL  fclr(m,n,np+1),fall(m,n,np+1)
2202        real  fsdir(m,n),fsdif(m,n)        _RL  fsdir(m,n),fsdif(m,n)
2203    
2204  c-----initialize all-sky flux (fall) and surface downward fluxes  c-----initialize all-sky flux (fall) and surface downward fluxes
2205    
# Line 2546  c     due to co2 absorption. Line 2499  c     due to co2 absorption.
2499  c-----input parameters  c-----input parameters
2500    
2501        integer m,n,np        integer m,n,np
2502        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)
2503    
2504  c-----output (undated) parameter  c-----output (undated) parameter
2505    
2506        real  df(m,n,np+1)        _RL  df(m,n,np+1)
2507    
2508  c-----temporary array  c-----temporary array
2509    
2510        integer i,j,k,ic,iw        integer i,j,k,ic,iw
2511        real  xx,clog,wlog,dc,dw,x1,x2,y2        _RL  xx,clog1,wlog,dc,dw,x1,x2,y2
2512    
2513  c********************************************************************  c********************************************************************
2514  c-----include co2 look-up table  c-----include co2 look-up table
2515    
2516        include 'cah.dat'  #include "cah-dat.h"
2517        save cah  c     save cah
2518    
2519  c********************************************************************  c********************************************************************
2520  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 2525  c     extraterrestrial solar flux in the
2525         do j= 1, n         do j= 1, n
2526          do i= 1, m          do i= 1, m
2527            xx=1./.3            xx=1./.3
2528            clog=log10(swc(i,j,k)*csm(i,j))            clog1=log10(swc(i,j,k)*csm(i,j))
2529            wlog=log10(swh(i,j,k)*csm(i,j))            wlog=log10(swh(i,j,k)*csm(i,j))
2530            ic=int( (clog+3.15)*xx+1.)            ic=int( (clog1+3.15)*xx+1.)
2531            iw=int( (wlog+4.15)*xx+1.)            iw=int( (wlog+4.15)*xx+1.)
2532            if(ic.lt.2)ic=2            if(ic.lt.2)ic=2
2533            if(iw.lt.2)iw=2            if(iw.lt.2)iw=2
2534            if(ic.gt.22)ic=22            if(ic.gt.22)ic=22
2535            if(iw.gt.19)iw=19                if(iw.gt.19)iw=19    
2536            dc=clog-float(ic-2)*.3+3.            dc=clog1-float(ic-2)*.3+3.
2537            dw=wlog-float(iw-2)*.3+4.              dw=wlog-float(iw-2)*.3+4.  
2538            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
2539            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.24
  ViewVC Help
Powered by ViewVC 1.1.22