/[lmdze]/trunk/Sources/phylmd/Radlwsw/swclr.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/Radlwsw/swclr.f

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

revision 217 by guez, Thu Mar 30 14:25:18 2017 UTC revision 218 by guez, Thu Mar 30 15:37:51 2017 UTC
# Line 6  contains Line 6  contains
6    
7    SUBROUTINE swclr(knu, flag_aer, palbp, pdsig, prayl, psec, pcgaz, ppizaz, &    SUBROUTINE swclr(knu, flag_aer, palbp, pdsig, prayl, psec, pcgaz, ppizaz, &
8         pray1, pray2, prefz, prj, prk, prmu0, ptauaz, ptra1, ptra2)         pray1, pray2, prefz, prj, prk, prmu0, ptauaz, ptra1, ptra2)
       
     USE raddim, only: kdlon, kflev  
     USE radepsi, only: repsct, zepsec  
     USE radopt, only: novlp  
9    
     ! ------------------------------------------------------------------  
10      ! PURPOSE.      ! PURPOSE.
     ! --------  
11      ! COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY IN CASE OF      ! COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY IN CASE OF
12      ! CLEAR-SKY COLUMN      ! CLEAR-SKY COLUMN
13    
14      ! REFERENCE.      ! REFERENCE.
     ! ----------  
   
15      ! SEE RADIATION'S PART OF THE ECMWF RESEARCH DEPARTMENT      ! SEE RADIATION'S PART OF THE ECMWF RESEARCH DEPARTMENT
16      ! DOCUMENTATION, AND FOUQUART AND BONNEL (1980)      ! DOCUMENTATION, AND FOUQUART AND BONNEL (1980)
17    
18      ! AUTHOR.      ! AUTHOR.
19      ! -------      ! JEAN-JACQUES MORCRETTE *ECMWF*
     ! JEAN-JACQUES MORCRETTE  *ECMWF*  
20    
21      ! MODIFICATIONS.      ! MODIFICATIONS.
     ! --------------  
22      ! ORIGINAL : 94-11-15      ! ORIGINAL : 94-11-15
23      ! ------------------------------------------------------------------  
24      ! * ARGUMENTS:      USE raddim, only: kdlon, kflev
25        USE radepsi, only: zepsec
26        USE radopt, only: novlp
27    
28        ! ARGUMENTS:
29    
30      INTEGER knu      INTEGER knu
31      ! -OB      ! -OB
32      logical, intent(in):: flag_aer      logical, intent(in):: flag_aer
33      DOUBLE PRECISION palbp(kdlon, 2)      DOUBLE PRECISION palbp(kdlon, 2)
34      DOUBLE PRECISION pdsig(kdlon, kflev)      DOUBLE PRECISION, intent(in):: pdsig(kdlon, kflev)
35      DOUBLE PRECISION prayl(kdlon)      DOUBLE PRECISION, intent(in):: prayl(kdlon)
36      DOUBLE PRECISION psec(kdlon)      DOUBLE PRECISION psec(kdlon)
37    
38      DOUBLE PRECISION pcgaz(kdlon, kflev)      DOUBLE PRECISION, intent(out):: pcgaz(kdlon, kflev)
39      DOUBLE PRECISION ppizaz(kdlon, kflev)      DOUBLE PRECISION, intent(out):: ppizaz(kdlon, kflev)
40      DOUBLE PRECISION pray1(kdlon, kflev+1)      DOUBLE PRECISION pray1(kdlon, kflev + 1)
41      DOUBLE PRECISION pray2(kdlon, kflev+1)      DOUBLE PRECISION pray2(kdlon, kflev + 1)
42      DOUBLE PRECISION prefz(kdlon, 2, kflev+1)      DOUBLE PRECISION prefz(kdlon, 2, kflev + 1)
43      DOUBLE PRECISION prj(kdlon, 6, kflev+1)      DOUBLE PRECISION prj(kdlon, 6, kflev + 1)
44      DOUBLE PRECISION prk(kdlon, 6, kflev+1)      DOUBLE PRECISION prk(kdlon, 6, kflev + 1)
45      DOUBLE PRECISION prmu0(kdlon, kflev+1)      DOUBLE PRECISION prmu0(kdlon, kflev + 1)
46      DOUBLE PRECISION ptauaz(kdlon, kflev)      DOUBLE PRECISION, intent(out):: ptauaz(kdlon, kflev)
47      DOUBLE PRECISION ptra1(kdlon, kflev+1)      DOUBLE PRECISION ptra1(kdlon, kflev + 1)
48      DOUBLE PRECISION ptra2(kdlon, kflev+1)      DOUBLE PRECISION ptra2(kdlon, kflev + 1)
   
     ! * LOCAL VARIABLES:  
49    
50      DOUBLE PRECISION zc0i(kdlon, kflev+1)      ! LOCAL VARIABLES:
51        DOUBLE PRECISION zc0i(kdlon, kflev + 1)
52      DOUBLE PRECISION zclear(kdlon)      DOUBLE PRECISION zclear(kdlon)
53      DOUBLE PRECISION zr21(kdlon)      DOUBLE PRECISION zr21(kdlon)
54      DOUBLE PRECISION zss0(kdlon)      DOUBLE PRECISION zss0(kdlon)
55      DOUBLE PRECISION zscat(kdlon)      DOUBLE PRECISION zscat(kdlon)
56      DOUBLE PRECISION ztr(kdlon, 2, kflev+1)      DOUBLE PRECISION ztr(kdlon, 2, kflev + 1)
   
57      INTEGER jl, jk, ja, jkl, jklp1, jaj, jkm1      INTEGER jl, jk, ja, jkl, jklp1, jaj, jkm1
58      DOUBLE PRECISION ztray, zgar, zratio, zff, zfacoa, zcorae      DOUBLE PRECISION zfacoa, zcorae
59      DOUBLE PRECISION zmue, zgap, zww, zto, zden, zmu1, zden1      DOUBLE PRECISION zmue, zgap, zww, zto, zden, zmu1, zden1
60      DOUBLE PRECISION zbmu0, zbmu1, zre11      DOUBLE PRECISION zbmu0, zbmu1, zre11
61        double precision, parameter:: REPSCT = 1d-10
62    
63      ! ------------------------------------------------------------------      !------------------------------------------------------------------
64        
65      ! *         1.    OPTICAL PARAMETERS FOR AEROSOLS AND RAYLEIGH      ! 1. OPTICAL PARAMETERS FOR AEROSOLS AND RAYLEIGH
     ! --------------------------------------------  
   
66    
67      DO jk = 1, kflev + 1      DO jk = 1, kflev + 1
68         DO ja = 1, 6         DO ja = 1, 6
69            DO jl = 1, kdlon            DO jl = 1, kdlon
70               prj(jl, ja, jk) = 0.               prj(jl, ja, jk) = 0d0
71               prk(jl, ja, jk) = 0.               prk(jl, ja, jk) = 0d0
72            END DO            END DO
73         END DO         END DO
74      END DO      END DO
75    
76      DO jk = 1, kflev      DO jk = 1, kflev
77         DO jl = 1, kdlon         DO jl = 1, kdlon
           ptauaz(jl, jk) = 0d0  
           ppizaz(jl, jk) = 0d0  
78            pcgaz(jl, jk) = 0d0            pcgaz(jl, jk) = 0d0
79              ptauaz(jl, jk) = prayl(jl) * pdsig(jl, jk)
80         END DO         END DO
81    
82         IF (flag_aer) THEN         IF (flag_aer) THEN
83            ! -OB            ! -OB
84            DO jl = 1, kdlon            DO jl = 1, kdlon
85               ztray = prayl(jl)*pdsig(jl, jk)               ppizaz(jl, jk) = 1d0
              zratio = ztray/(ztray+ptauaz(jl,jk))  
              zgar = pcgaz(jl, jk)  
              zff = zgar*zgar  
              ptauaz(jl, jk) = ztray + ptauaz(jl, jk)*(1.-ppizaz(jl,jk)*zff)  
              pcgaz(jl, jk) = zgar*(1.-zratio)/(1.+zgar)  
              ppizaz(jl, jk) = zratio + (1.-zratio)*ppizaz(jl, jk)*(1.-zff)/(1.- &  
                   ppizaz(jl,jk)*zff)  
86            END DO            END DO
87         ELSE         ELSE
88            DO jl = 1, kdlon            DO jl = 1, kdlon
89               ztray = prayl(jl)*pdsig(jl, jk)               ppizaz(jl, jk) = 1d0 - repsct
              ptauaz(jl, jk) = ztray  
              pcgaz(jl, jk) = 0.  
              ppizaz(jl, jk) = 1. - repsct  
90            END DO            END DO
91         END IF         END IF
92      END DO      END DO
93    
94      ! ------------------------------------------------------------------      ! 2. TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL
   
     ! *         2.    TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL  
     ! ----------------------------------------------  
   
95    
96      DO jl = 1, kdlon      DO jl = 1, kdlon
97         zc0i(jl, kflev+1) = 0.         zc0i(jl, kflev + 1) = 0d0
98         zclear(jl) = 1.         zclear(jl) = 1d0
99         zscat(jl) = 0.         zscat(jl) = 0d0
100      END DO      END DO
101    
102      jk = 1      jk = 1
103      jkl = kflev + 1 - jk      jkl = kflev + 1 - jk
104      jklp1 = jkl + 1      jklp1 = jkl + 1
105      DO jl = 1, kdlon      DO jl = 1, kdlon
106         zfacoa = 1. - ppizaz(jl, jkl)*pcgaz(jl, jkl)*pcgaz(jl, jkl)         zfacoa = 1d0
107         zcorae = zfacoa*ptauaz(jl, jkl)*psec(jl)         zcorae = zfacoa * ptauaz(jl, jkl) * psec(jl)
108         zr21(jl) = exp(-zcorae)         zr21(jl) = exp(- zcorae)
109         zss0(jl) = 1. - zr21(jl)         zss0(jl) = 1d0 - zr21(jl)
110    
111         IF (novlp==1) THEN         IF (novlp == 1) THEN
112            ! * maximum-random            ! maximum-random
113            zclear(jl) = zclear(jl)*(1.0-max(zss0(jl),zscat(jl)))/ &            zclear(jl) = zclear(jl) * (1d0 - max(zss0(jl), zscat(jl))) / (1d0 &
114                 (1.0-min(zscat(jl),1.-zepsec))                 - min(zscat(jl), 1d0 - zepsec))
115            zc0i(jl, jkl) = 1.0 - zclear(jl)            zc0i(jl, jkl) = 1d0 - zclear(jl)
116            zscat(jl) = zss0(jl)            zscat(jl) = zss0(jl)
117         ELSE IF (novlp==2) THEN         ELSE IF (novlp == 2) THEN
118            ! * maximum            ! maximum
119            zscat(jl) = max(zss0(jl), zscat(jl))            zscat(jl) = max(zss0(jl), zscat(jl))
120            zc0i(jl, jkl) = zscat(jl)            zc0i(jl, jkl) = zscat(jl)
121         ELSE IF (novlp==3) THEN         ELSE IF (novlp == 3) THEN
122            ! * random            ! random
123            zclear(jl) = zclear(jl)*(1.0-zss0(jl))            zclear(jl) = zclear(jl) * (1d0 - zss0(jl))
124            zscat(jl) = 1.0 - zclear(jl)            zscat(jl) = 1d0 - zclear(jl)
125            zc0i(jl, jkl) = zscat(jl)            zc0i(jl, jkl) = zscat(jl)
126         END IF         END IF
127      END DO      END DO
# Line 154  contains Line 130  contains
130         jkl = kflev + 1 - jk         jkl = kflev + 1 - jk
131         jklp1 = jkl + 1         jklp1 = jkl + 1
132         DO jl = 1, kdlon         DO jl = 1, kdlon
133            zfacoa = 1. - ppizaz(jl, jkl)*pcgaz(jl, jkl)*pcgaz(jl, jkl)            zfacoa = 1d0
134            zcorae = zfacoa*ptauaz(jl, jkl)*psec(jl)            zcorae = zfacoa * ptauaz(jl, jkl) * psec(jl)
135            zr21(jl) = exp(-zcorae)            zr21(jl) = exp(- zcorae)
136            zss0(jl) = 1. - zr21(jl)            zss0(jl) = 1d0 - zr21(jl)
137    
138            IF (novlp==1) THEN            IF (novlp == 1) THEN
139               ! * maximum-random               ! maximum-random
140               zclear(jl) = zclear(jl)*(1.0-max(zss0(jl),zscat(jl)))/ &               zclear(jl) = zclear(jl) * (1d0 - max(zss0(jl), zscat(jl))) &
141                    (1.0-min(zscat(jl),1.-zepsec))                    / (1d0 - min(zscat(jl), 1d0 - zepsec))
142               zc0i(jl, jkl) = 1.0 - zclear(jl)               zc0i(jl, jkl) = 1d0 - zclear(jl)
143               zscat(jl) = zss0(jl)               zscat(jl) = zss0(jl)
144            ELSE IF (novlp==2) THEN            ELSE IF (novlp == 2) THEN
145               ! * maximum               ! maximum
146               zscat(jl) = max(zss0(jl), zscat(jl))               zscat(jl) = max(zss0(jl), zscat(jl))
147               zc0i(jl, jkl) = zscat(jl)               zc0i(jl, jkl) = zscat(jl)
148            ELSE IF (novlp==3) THEN            ELSE IF (novlp == 3) THEN
149               ! * random               ! random
150               zclear(jl) = zclear(jl)*(1.0-zss0(jl))               zclear(jl) = zclear(jl) * (1d0 - zss0(jl))
151               zscat(jl) = 1.0 - zclear(jl)               zscat(jl) = 1d0 - zclear(jl)
152               zc0i(jl, jkl) = zscat(jl)               zc0i(jl, jkl) = zscat(jl)
153            END IF            END IF
154         END DO         END DO
155      END DO      END DO
156    
157      ! ------------------------------------------------------------------      ! 3. REFLECTIVITY/TRANSMISSIVITY FOR PURE SCATTERING
   
     ! *         3.    REFLECTIVITY/TRANSMISSIVITY FOR PURE SCATTERING  
     ! -----------------------------------------------  
   
158    
159      DO jl = 1, kdlon      DO jl = 1, kdlon
160         pray1(jl, kflev+1) = 0.         pray1(jl, kflev + 1) = 0d0
161         pray2(jl, kflev+1) = 0.         pray2(jl, kflev + 1) = 0d0
162         prefz(jl, 2, 1) = palbp(jl, knu)         prefz(jl, 2, 1) = palbp(jl, knu)
163         prefz(jl, 1, 1) = palbp(jl, knu)         prefz(jl, 1, 1) = palbp(jl, knu)
164         ptra1(jl, kflev+1) = 1.         ptra1(jl, kflev + 1) = 1d0
165         ptra2(jl, kflev+1) = 1.         ptra2(jl, kflev + 1) = 1d0
166      END DO      END DO
167    
168      DO jk = 2, kflev + 1      DO jk = 2, kflev + 1
169         jkm1 = jk - 1         jkm1 = jk - 1
170         DO jl = 1, kdlon         DO jl = 1, kdlon
171    
172              ! 3.1 EQUIVALENT ZENITH ANGLE
173    
174            ! ------------------------------------------------------------------            zmue = (1d0 - zc0i(jl, jk)) * psec(jl) + zc0i(jl, jk) * 1.66d0
175              prmu0(jl, jk) = 1d0 / zmue
176    
177            ! *         3.1  EQUIVALENT ZENITH ANGLE            ! 3.2 REFLECT./TRANSMISSIVITY DUE TO RAYLEIGH AND AEROSOLS
           ! -----------------------  
178    
179              zgap = 0d0
180            zmue = (1.-zc0i(jl,jk))*psec(jl) + zc0i(jl, jk)*1.66            zbmu0 = 0.5d0 - 0.75d0 * zgap / zmue
           prmu0(jl, jk) = 1./zmue  
   
   
           ! ------------------------------------------------------------------  
   
           ! *         3.2  REFLECT./TRANSMISSIVITY DUE TO RAYLEIGH AND AEROSOLS  
           ! ----------------------------------------------------  
   
   
           zgap = pcgaz(jl, jkm1)  
           zbmu0 = 0.5 - 0.75*zgap/zmue  
181            zww = ppizaz(jl, jkm1)            zww = ppizaz(jl, jkm1)
182            zto = ptauaz(jl, jkm1)            zto = ptauaz(jl, jkm1)
183            zden = 1. + (1.-zww+zbmu0*zww)*zto*zmue + (1-zww)*(1.-zww+2.*zbmu0*zww) &            zden = 1d0 + (1d0 - zww + zbmu0 * zww) * zto * zmue + (1d0 - zww) &
184                 *zto*zto*zmue*zmue                 * (1d0 - zww + 2d0 * zbmu0 * zww) * zto * zto * zmue * zmue
185            pray1(jl, jkm1) = zbmu0*zww*zto*zmue/zden            pray1(jl, jkm1) = zbmu0 * zww * zto * zmue / zden
186            ptra1(jl, jkm1) = 1./zden            ptra1(jl, jkm1) = 1d0 / zden
187    
188            zmu1 = 0.5            zmu1 = 0.5d0
189            zbmu1 = 0.5 - 0.75*zgap*zmu1            zbmu1 = 0.5d0 - 0.75d0 * zgap * zmu1
190            zden1 = 1. + (1.-zww+zbmu1*zww)*zto/zmu1 + (1-zww)*(1.-zww+2.*zbmu1*zww &            zden1 = 1d0 + (1d0 - zww + zbmu1 * zww) * zto / zmu1 + (1d0 - zww) &
191                 )*zto*zto/zmu1/zmu1                 * (1d0 - zww + 2d0 * zbmu1 * zww &
192            pray2(jl, jkm1) = zbmu1*zww*zto/zmu1/zden1                 ) * zto * zto / zmu1 / zmu1
193            ptra2(jl, jkm1) = 1./zden1            pray2(jl, jkm1) = zbmu1 * zww * zto / zmu1 / zden1
194              ptra2(jl, jkm1) = 1d0 / zden1
195    
196              prefz(jl, 1, jk) = (pray1(jl, jkm1) + prefz(jl, 1, jkm1) &
197                   * ptra1(jl, jkm1)* ptra2(jl, jkm1) / (1d0 - pray2(jl, jkm1) &
198                   * prefz(jl, 1, jkm1)))
199    
200              ztr(jl, 1, jkm1) = (ptra1(jl, jkm1) / (1d0 - pray2(jl, jkm1) &
201                   * prefz(jl, 1, jkm1)))
202    
203            prefz(jl, 1, jk) = (pray1(jl,jkm1)+prefz(jl,1,jkm1)*ptra1(jl,jkm1)* &            prefz(jl, 2, jk) = (pray1(jl, jkm1) + prefz(jl, 2, jkm1) &
204                 ptra2(jl,jkm1)/(1.-pray2(jl,jkm1)*prefz(jl,1,jkm1)))                 * ptra1(jl, jkm1) * ptra2(jl, jkm1))
   
           ztr(jl, 1, jkm1) = (ptra1(jl,jkm1)/(1.-pray2(jl,jkm1)*prefz(jl,1, &  
                jkm1)))  
   
           prefz(jl, 2, jk) = (pray1(jl,jkm1)+prefz(jl,2,jkm1)*ptra1(jl,jkm1)* &  
                ptra2(jl,jkm1))  
205    
206            ztr(jl, 2, jkm1) = ptra1(jl, jkm1)            ztr(jl, 2, jkm1) = ptra1(jl, jkm1)
207    
208         END DO         END DO
209      END DO      END DO
210      DO jl = 1, kdlon      DO jl = 1, kdlon
211         zmue = (1.-zc0i(jl,1))*psec(jl) + zc0i(jl, 1)*1.66         zmue = (1d0 - zc0i(jl, 1)) * psec(jl) + zc0i(jl, 1) * 1.66d0
212         prmu0(jl, 1) = 1./zmue         prmu0(jl, 1) = 1d0 / zmue
213      END DO      END DO
214    
215        ! 3.5 REFLECT./TRANSMISSIVITY BETWEEN SURFACE AND LEVEL
216    
217      ! ------------------------------------------------------------------      IF (knu == 1) THEN
   
     ! *         3.5    REFLECT./TRANSMISSIVITY BETWEEN SURFACE AND LEVEL  
     ! -------------------------------------------------  
   
   
     IF (knu==1) THEN  
218         jaj = 2         jaj = 2
219         DO jl = 1, kdlon         DO jl = 1, kdlon
220            prj(jl, jaj, kflev+1) = 1.            prj(jl, jaj, kflev + 1) = 1d0
221            prk(jl, jaj, kflev+1) = prefz(jl, 1, kflev+1)            prk(jl, jaj, kflev + 1) = prefz(jl, 1, kflev + 1)
222         END DO         END DO
223    
224         DO jk = 1, kflev         DO jk = 1, kflev
225            jkl = kflev + 1 - jk            jkl = kflev + 1 - jk
226            jklp1 = jkl + 1            jklp1 = jkl + 1
227            DO jl = 1, kdlon            DO jl = 1, kdlon
228               zre11 = prj(jl, jaj, jklp1)*ztr(jl, 1, jkl)               zre11 = prj(jl, jaj, jklp1) * ztr(jl, 1, jkl)
229               prj(jl, jaj, jkl) = zre11               prj(jl, jaj, jkl) = zre11
230               prk(jl, jaj, jkl) = zre11*prefz(jl, 1, jkl)               prk(jl, jaj, jkl) = zre11 * prefz(jl, 1, jkl)
231            END DO            END DO
232         END DO         END DO
   
233      ELSE      ELSE
   
234         DO jaj = 1, 2         DO jaj = 1, 2
235            DO jl = 1, kdlon            DO jl = 1, kdlon
236               prj(jl, jaj, kflev+1) = 1.               prj(jl, jaj, kflev + 1) = 1d0
237               prk(jl, jaj, kflev+1) = prefz(jl, jaj, kflev+1)               prk(jl, jaj, kflev + 1) = prefz(jl, jaj, kflev + 1)
238            END DO            END DO
239    
240            DO jk = 1, kflev            DO jk = 1, kflev
241               jkl = kflev + 1 - jk               jkl = kflev + 1 - jk
242               jklp1 = jkl + 1               jklp1 = jkl + 1
243               DO jl = 1, kdlon               DO jl = 1, kdlon
244                  zre11 = prj(jl, jaj, jklp1)*ztr(jl, jaj, jkl)                  zre11 = prj(jl, jaj, jklp1) * ztr(jl, jaj, jkl)
245                  prj(jl, jaj, jkl) = zre11                  prj(jl, jaj, jkl) = zre11
246                  prk(jl, jaj, jkl) = zre11*prefz(jl, jaj, jkl)                  prk(jl, jaj, jkl) = zre11 * prefz(jl, jaj, jkl)
247               END DO               END DO
248            END DO            END DO
249         END DO         END DO
   
250      END IF      END IF
251    
252    END SUBROUTINE swclr    END SUBROUTINE swclr
Legend:
Removed from v.217  
changed lines
  Added in v.218
  ViewVC Help
Powered by ViewVC 1.1.21