/[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 207 - (show annotations)
Thu Sep 1 10:30:53 2016 UTC (7 years, 8 months ago) by guez
File size: 10204 byte(s)
New philosophy on compiler options.

Removed source code for thermcep = f. (Not used in LMDZ either.)

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

  ViewVC Help
Powered by ViewVC 1.1.21