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

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