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

Contents of /trunk/Sources/phylmd/Radlwsw/swr.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (show annotations)
Wed Mar 5 14:57:53 2014 UTC (10 years, 2 months ago) by guez
Original Path: trunk/phylmd/Radlwsw/swr.f
File size: 9586 byte(s)
Changed all ".f90" suffixes to ".f".
1 SUBROUTINE swr(knu, palbd, pcg, pcld, pdsig, pomega, prayl, psec, ptau, &
2 pcgaz, ppizaz, pray1, pray2, prefz, prj, prk, prmue, ptauaz, ptra1, &
3 ptra2)
4 USE dimens_m
5 USE dimphy
6 USE raddim
7 USE radepsi
8 USE radopt
9 IMPLICIT NONE
10
11 ! ------------------------------------------------------------------
12 ! PURPOSE.
13 ! --------
14 ! COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY IN CASE OF
15 ! CONTINUUM SCATTERING
16
17 ! METHOD.
18 ! -------
19
20 ! 1. COMPUTES CONTINUUM FLUXES CORRESPONDING TO AEROSOL
21 ! OR/AND RAYLEIGH SCATTERING (NO MOLECULAR GAS ABSORPTION)
22
23 ! REFERENCE.
24 ! ----------
25
26 ! SEE RADIATION'S PART OF THE ECMWF RESEARCH DEPARTMENT
27 ! DOCUMENTATION, AND FOUQUART AND BONNEL (1980)
28
29 ! AUTHOR.
30 ! -------
31 ! JEAN-JACQUES MORCRETTE *ECMWF*
32
33 ! MODIFICATIONS.
34 ! --------------
35 ! ORIGINAL : 89-07-14
36 ! ------------------------------------------------------------------
37 ! * ARGUMENTS:
38
39 INTEGER knu
40 DOUBLE PRECISION palbd(kdlon, 2)
41 DOUBLE PRECISION pcg(kdlon, 2, kflev)
42 DOUBLE PRECISION pcld(kdlon, kflev)
43 DOUBLE PRECISION pdsig(kdlon, kflev)
44 DOUBLE PRECISION pomega(kdlon, 2, kflev)
45 DOUBLE PRECISION prayl(kdlon)
46 DOUBLE PRECISION psec(kdlon)
47 DOUBLE PRECISION ptau(kdlon, 2, kflev)
48
49 DOUBLE PRECISION pray1(kdlon, kflev+1)
50 DOUBLE PRECISION pray2(kdlon, kflev+1)
51 DOUBLE PRECISION prefz(kdlon, 2, kflev+1)
52 DOUBLE PRECISION prj(kdlon, 6, kflev+1)
53 DOUBLE PRECISION prk(kdlon, 6, kflev+1)
54 DOUBLE PRECISION prmue(kdlon, kflev+1)
55 DOUBLE PRECISION pcgaz(kdlon, kflev)
56 DOUBLE PRECISION ppizaz(kdlon, kflev)
57 DOUBLE PRECISION ptauaz(kdlon, kflev)
58 DOUBLE PRECISION ptra1(kdlon, kflev+1)
59 DOUBLE PRECISION ptra2(kdlon, kflev+1)
60
61 ! * LOCAL VARIABLES:
62
63 DOUBLE PRECISION zc1i(kdlon, kflev+1)
64 DOUBLE PRECISION zcleq(kdlon, kflev)
65 DOUBLE PRECISION zclear(kdlon)
66 DOUBLE PRECISION zcloud(kdlon)
67 DOUBLE PRECISION zgg(kdlon)
68 DOUBLE PRECISION zref(kdlon)
69 DOUBLE PRECISION zre1(kdlon)
70 DOUBLE PRECISION zre2(kdlon)
71 DOUBLE PRECISION zrmuz(kdlon)
72 DOUBLE PRECISION zrneb(kdlon)
73 DOUBLE PRECISION zr21(kdlon)
74 DOUBLE PRECISION zr22(kdlon)
75 DOUBLE PRECISION zr23(kdlon)
76 DOUBLE PRECISION zss1(kdlon)
77 DOUBLE PRECISION zto1(kdlon)
78 DOUBLE PRECISION ztr(kdlon, 2, kflev+1)
79 DOUBLE PRECISION ztr1(kdlon)
80 DOUBLE PRECISION ztr2(kdlon)
81 DOUBLE PRECISION zw(kdlon)
82
83 INTEGER jk, jl, ja, jkl, jklp1, jkm1, jaj
84 DOUBLE PRECISION zfacoa, zfacoc, zcorae, zcorcd
85 DOUBLE PRECISION zmue, zgap, zww, zto, zden, zden1
86 DOUBLE PRECISION zmu1, zre11, zbmu0, zbmu1
87
88 ! ------------------------------------------------------------------
89
90 ! * 1. INITIALIZATION
91 ! --------------
92
93
94 DO jk = 1, kflev + 1
95 DO ja = 1, 6
96 DO jl = 1, kdlon
97 prj(jl, ja, jk) = 0.
98 prk(jl, ja, jk) = 0.
99 END DO
100 END DO
101 END DO
102
103
104 ! ------------------------------------------------------------------
105
106 ! * 2. TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL
107 ! ----------------------------------------------
108
109
110 DO jl = 1, kdlon
111 zr23(jl) = 0.
112 zc1i(jl, kflev+1) = 0.
113 zclear(jl) = 1.
114 zcloud(jl) = 0.
115 END DO
116
117 jk = 1
118 jkl = kflev + 1 - jk
119 jklp1 = jkl + 1
120 DO jl = 1, kdlon
121 zfacoa = 1. - ppizaz(jl, jkl)*pcgaz(jl, jkl)*pcgaz(jl, jkl)
122 zfacoc = 1. - pomega(jl, knu, jkl)*pcg(jl, knu, jkl)*pcg(jl, knu, jkl)
123 zcorae = zfacoa*ptauaz(jl, jkl)*psec(jl)
124 zcorcd = zfacoc*ptau(jl, knu, jkl)*psec(jl)
125 zr21(jl) = exp(-zcorae)
126 zr22(jl) = exp(-zcorcd)
127 zss1(jl) = pcld(jl, jkl)*(1.0-zr21(jl)*zr22(jl)) + &
128 (1.0-pcld(jl,jkl))*(1.0-zr21(jl))
129 zcleq(jl, jkl) = zss1(jl)
130
131 IF (novlp==1) THEN
132 ! * maximum-random
133 zclear(jl) = zclear(jl)*(1.0-max(zss1(jl),zcloud(jl)))/ &
134 (1.0-min(zcloud(jl),1.-zepsec))
135 zc1i(jl, jkl) = 1.0 - zclear(jl)
136 zcloud(jl) = zss1(jl)
137 ELSE IF (novlp==2) THEN
138 ! * maximum
139 zcloud(jl) = max(zss1(jl), zcloud(jl))
140 zc1i(jl, jkl) = zcloud(jl)
141 ELSE IF (novlp==3) THEN
142 ! * random
143 zclear(jl) = zclear(jl)*(1.0-zss1(jl))
144 zcloud(jl) = 1.0 - zclear(jl)
145 zc1i(jl, jkl) = zcloud(jl)
146 END IF
147 END DO
148
149 DO jk = 2, kflev
150 jkl = kflev + 1 - jk
151 jklp1 = jkl + 1
152 DO jl = 1, kdlon
153 zfacoa = 1. - ppizaz(jl, jkl)*pcgaz(jl, jkl)*pcgaz(jl, jkl)
154 zfacoc = 1. - pomega(jl, knu, jkl)*pcg(jl, knu, jkl)*pcg(jl, knu, jkl)
155 zcorae = zfacoa*ptauaz(jl, jkl)*psec(jl)
156 zcorcd = zfacoc*ptau(jl, knu, jkl)*psec(jl)
157 zr21(jl) = exp(-zcorae)
158 zr22(jl) = exp(-zcorcd)
159 zss1(jl) = pcld(jl, jkl)*(1.0-zr21(jl)*zr22(jl)) + &
160 (1.0-pcld(jl,jkl))*(1.0-zr21(jl))
161 zcleq(jl, jkl) = zss1(jl)
162
163 IF (novlp==1) THEN
164 ! * maximum-random
165 zclear(jl) = zclear(jl)*(1.0-max(zss1(jl),zcloud(jl)))/ &
166 (1.0-min(zcloud(jl),1.-zepsec))
167 zc1i(jl, jkl) = 1.0 - zclear(jl)
168 zcloud(jl) = zss1(jl)
169 ELSE IF (novlp==2) THEN
170 ! * maximum
171 zcloud(jl) = max(zss1(jl), zcloud(jl))
172 zc1i(jl, jkl) = zcloud(jl)
173 ELSE IF (novlp==3) THEN
174 ! * random
175 zclear(jl) = zclear(jl)*(1.0-zss1(jl))
176 zcloud(jl) = 1.0 - zclear(jl)
177 zc1i(jl, jkl) = zcloud(jl)
178 END IF
179 END DO
180 END DO
181
182 ! ------------------------------------------------------------------
183
184 ! * 3. REFLECTIVITY/TRANSMISSIVITY FOR PURE SCATTERING
185 ! -----------------------------------------------
186
187
188 DO jl = 1, kdlon
189 pray1(jl, kflev+1) = 0.
190 pray2(jl, kflev+1) = 0.
191 prefz(jl, 2, 1) = palbd(jl, knu)
192 prefz(jl, 1, 1) = palbd(jl, knu)
193 ptra1(jl, kflev+1) = 1.
194 ptra2(jl, kflev+1) = 1.
195 END DO
196
197 DO jk = 2, kflev + 1
198 jkm1 = jk - 1
199 DO jl = 1, kdlon
200 zrneb(jl) = pcld(jl, jkm1)
201 zre1(jl) = 0.
202 ztr1(jl) = 0.
203 zre2(jl) = 0.
204 ztr2(jl) = 0.
205
206
207 ! ------------------------------------------------------------------
208
209 ! * 3.1 EQUIVALENT ZENITH ANGLE
210 ! -----------------------
211
212
213 zmue = (1.-zc1i(jl,jk))*psec(jl) + zc1i(jl, jk)*1.66
214 prmue(jl, jk) = 1./zmue
215
216
217 ! ------------------------------------------------------------------
218
219 ! * 3.2 REFLECT./TRANSMISSIVITY DUE TO RAYLEIGH AND AEROSOLS
220 ! ----------------------------------------------------
221
222
223 zgap = pcgaz(jl, jkm1)
224 zbmu0 = 0.5 - 0.75*zgap/zmue
225 zww = ppizaz(jl, jkm1)
226 zto = ptauaz(jl, jkm1)
227 zden = 1. + (1.-zww+zbmu0*zww)*zto*zmue + (1-zww)*(1.-zww+2.*zbmu0*zww) &
228 *zto*zto*zmue*zmue
229 pray1(jl, jkm1) = zbmu0*zww*zto*zmue/zden
230 ptra1(jl, jkm1) = 1./zden
231 ! PRINT *,' LOOP 342 ** 3 ** JL=',JL,PRAY1(JL,JKM1),PTRA1(JL,JKM1)
232
233 zmu1 = 0.5
234 zbmu1 = 0.5 - 0.75*zgap*zmu1
235 zden1 = 1. + (1.-zww+zbmu1*zww)*zto/zmu1 + (1-zww)*(1.-zww+2.*zbmu1*zww &
236 )*zto*zto/zmu1/zmu1
237 pray2(jl, jkm1) = zbmu1*zww*zto/zmu1/zden1
238 ptra2(jl, jkm1) = 1./zden1
239
240
241 ! ------------------------------------------------------------------
242
243 ! * 3.3 EFFECT OF CLOUD LAYER
244 ! ---------------------
245
246
247 zw(jl) = pomega(jl, knu, jkm1)
248 zto1(jl) = ptau(jl, knu, jkm1)/zw(jl) + ptauaz(jl, jkm1)/ppizaz(jl, &
249 jkm1)
250 zr21(jl) = ptau(jl, knu, jkm1) + ptauaz(jl, jkm1)
251 zr22(jl) = ptau(jl, knu, jkm1)/zr21(jl)
252 zgg(jl) = zr22(jl)*pcg(jl, knu, jkm1) + (1.-zr22(jl))*pcgaz(jl, jkm1)
253 ! Modif PhD - JJM 19/03/96 pour erreurs arrondis
254 ! machine
255 ! PHD PROTECTION ZW(JL) = ZR21(JL) / ZTO1(JL)
256 IF (zw(jl)==1. .AND. ppizaz(jl,jkm1)==1.) THEN
257 zw(jl) = 1.
258 ELSE
259 zw(jl) = zr21(jl)/zto1(jl)
260 END IF
261 zref(jl) = prefz(jl, 1, jkm1)
262 zrmuz(jl) = prmue(jl, jk)
263 END DO
264
265 CALL swde(zgg, zref, zrmuz, zto1, zw, zre1, zre2, ztr1, ztr2)
266
267 DO jl = 1, kdlon
268
269 prefz(jl, 1, jk) = (1.-zrneb(jl))*(pray1(jl,jkm1)+prefz(jl,1,jkm1)* &
270 ptra1(jl,jkm1)*ptra2(jl,jkm1)/(1.-pray2(jl,jkm1)*prefz(jl,1, &
271 jkm1))) + zrneb(jl)*zre2(jl)
272
273 ztr(jl, 1, jkm1) = zrneb(jl)*ztr2(jl) + (ptra1(jl,jkm1)/(1.-pray2(jl, &
274 jkm1)*prefz(jl,1,jkm1)))*(1.-zrneb(jl))
275
276 prefz(jl, 2, jk) = (1.-zrneb(jl))*(pray1(jl,jkm1)+prefz(jl,2,jkm1)* &
277 ptra1(jl,jkm1)*ptra2(jl,jkm1)) + zrneb(jl)*zre1(jl)
278
279 ztr(jl, 2, jkm1) = zrneb(jl)*ztr1(jl) + ptra1(jl, jkm1)*(1.-zrneb(jl))
280
281 END DO
282 END DO
283 DO jl = 1, kdlon
284 zmue = (1.-zc1i(jl,1))*psec(jl) + zc1i(jl, 1)*1.66
285 prmue(jl, 1) = 1./zmue
286 END DO
287
288
289 ! ------------------------------------------------------------------
290
291 ! * 3.5 REFLECT./TRANSMISSIVITY BETWEEN SURFACE AND LEVEL
292 ! -------------------------------------------------
293
294
295 IF (knu==1) THEN
296 jaj = 2
297 DO jl = 1, kdlon
298 prj(jl, jaj, kflev+1) = 1.
299 prk(jl, jaj, kflev+1) = prefz(jl, 1, kflev+1)
300 END DO
301
302 DO jk = 1, kflev
303 jkl = kflev + 1 - jk
304 jklp1 = jkl + 1
305 DO jl = 1, kdlon
306 zre11 = prj(jl, jaj, jklp1)*ztr(jl, 1, jkl)
307 prj(jl, jaj, jkl) = zre11
308 prk(jl, jaj, jkl) = zre11*prefz(jl, 1, jkl)
309 END DO
310 END DO
311
312 ELSE
313
314 DO jaj = 1, 2
315 DO jl = 1, kdlon
316 prj(jl, jaj, kflev+1) = 1.
317 prk(jl, jaj, kflev+1) = prefz(jl, jaj, kflev+1)
318 END DO
319
320 DO jk = 1, kflev
321 jkl = kflev + 1 - jk
322 jklp1 = jkl + 1
323 DO jl = 1, kdlon
324 zre11 = prj(jl, jaj, jklp1)*ztr(jl, jaj, jkl)
325 prj(jl, jaj, jkl) = zre11
326 prk(jl, jaj, jkl) = zre11*prefz(jl, jaj, jkl)
327 END DO
328 END DO
329 END DO
330
331 END IF
332
333 ! ------------------------------------------------------------------
334
335 RETURN
336 END SUBROUTINE swr

  ViewVC Help
Powered by ViewVC 1.1.21