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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 157 - (show annotations)
Mon Jul 20 16:01:49 2015 UTC (8 years, 9 months ago) by guez
File size: 8497 byte(s)
Just encapsulated SUBROUTINE vlsplt in a module and cleaned it.

In procedure vlx, local variables dxqu and adxqu only need indices
iip2:ip1jm. Otherwise, just cleaned vlx.

Procedures dynredem0 and dynredem1 no longer have argument fichnom,
they just operate on a file named "restart.nc". The programming
guideline here is that gcm should not be more complex than it needs by
itself, other programs (ce0l etc.) just have to adapt to gcm. So ce0l
now creates files "restart.nc" and "restartphy.nc".

In order to facilitate decentralizing the writing of "restartphy.nc",
created a procedure phyredem0 out of phyredem. phyredem0 creates the
NetCDF header of "restartphy.nc" while phyredem writes the NetCDF
variables. As the global attribute itau_phy needs to be filled in
phyredem0, at the beginnig of the run, we must compute its value
instead of just using itap. So we have a dummy argument lmt_pas of
phyredem0. Also, the ncid of "startphy.nc" is upgraded from local
variable of phyetat0 to dummy argument. phyetat0 no longer closes
"startphy.nc".

Following the same decentralizing objective, the ncid of "restart.nc"
is upgraded from local variable of dynredem0 to module variable of
dynredem0_m. "restart.nc" is not closed at the end of dynredem0 nor
opened at the beginning of dynredem1.

In procedure etat0, instead of creating many vectors of size klon
which will be filled with zeroes, just create one array null_array.

In procedure phytrac, instead of writing trs(: 1) to a text file,
write it to "restartphy.nc" (following LMDZ). This is better because
now trs(: 1) is next to its coordinates. We can write to
"restartphy.nc" from phytrac directly, and not add trs(: 1) to the
long list of variables in physiq, thanks to the decentralizing of
"restartphy.nc".

In procedure phyetat0, we no longer write to standard output the
minimum and maximum values of read arrays. It is ok to check input and
abort on invalid values but just printing statistics on input seems too
much useless computation and out of place clutter.

1 SUBROUTINE swclr(knu, paer, flag_aer, tauae, pizae, cgae, palbp, pdsig, &
2 prayl, psec, pcgaz, ppizaz, pray1, pray2, prefz, prj, prk, prmu0, ptauaz, &
3 ptra1, 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 ! CLEAR-SKY COLUMN
16
17 ! REFERENCE.
18 ! ----------
19
20 ! SEE RADIATION'S PART OF THE ECMWF RESEARCH DEPARTMENT
21 ! DOCUMENTATION, AND FOUQUART AND BONNEL (1980)
22
23 ! AUTHOR.
24 ! -------
25 ! JEAN-JACQUES MORCRETTE *ECMWF*
26
27 ! MODIFICATIONS.
28 ! --------------
29 ! ORIGINAL : 94-11-15
30 ! ------------------------------------------------------------------
31 ! * ARGUMENTS:
32
33 INTEGER knu
34 ! -OB
35 DOUBLE PRECISION flag_aer
36 DOUBLE PRECISION tauae(kdlon, kflev, 2)
37 DOUBLE PRECISION pizae(kdlon, kflev, 2)
38 DOUBLE PRECISION cgae(kdlon, kflev, 2)
39 DOUBLE PRECISION paer(kdlon, kflev, 5)
40 DOUBLE PRECISION palbp(kdlon, 2)
41 DOUBLE PRECISION pdsig(kdlon, kflev)
42 DOUBLE PRECISION prayl(kdlon)
43 DOUBLE PRECISION psec(kdlon)
44
45 DOUBLE PRECISION pcgaz(kdlon, kflev)
46 DOUBLE PRECISION ppizaz(kdlon, kflev)
47 DOUBLE PRECISION pray1(kdlon, kflev+1)
48 DOUBLE PRECISION pray2(kdlon, kflev+1)
49 DOUBLE PRECISION prefz(kdlon, 2, kflev+1)
50 DOUBLE PRECISION prj(kdlon, 6, kflev+1)
51 DOUBLE PRECISION prk(kdlon, 6, kflev+1)
52 DOUBLE PRECISION prmu0(kdlon, kflev+1)
53 DOUBLE PRECISION ptauaz(kdlon, kflev)
54 DOUBLE PRECISION ptra1(kdlon, kflev+1)
55 DOUBLE PRECISION ptra2(kdlon, kflev+1)
56
57 ! * LOCAL VARIABLES:
58
59 DOUBLE PRECISION zc0i(kdlon, kflev+1)
60 DOUBLE PRECISION zcle0(kdlon, kflev)
61 DOUBLE PRECISION zclear(kdlon)
62 DOUBLE PRECISION zr21(kdlon)
63 DOUBLE PRECISION zr23(kdlon)
64 DOUBLE PRECISION zss0(kdlon)
65 DOUBLE PRECISION zscat(kdlon)
66 DOUBLE PRECISION ztr(kdlon, 2, kflev+1)
67
68 INTEGER jl, jk, ja, jkl, jklp1, jaj, jkm1
69 DOUBLE PRECISION ztray, zgar, zratio, zff, zfacoa, zcorae
70 DOUBLE PRECISION zmue, zgap, zww, zto, zden, zmu1, zden1
71 DOUBLE PRECISION zbmu0, zbmu1, zre11
72
73 ! ------------------------------------------------------------------
74
75 ! * 1. OPTICAL PARAMETERS FOR AEROSOLS AND RAYLEIGH
76 ! --------------------------------------------
77
78
79 DO jk = 1, kflev + 1
80 DO ja = 1, 6
81 DO jl = 1, kdlon
82 prj(jl, ja, jk) = 0.
83 prk(jl, ja, jk) = 0.
84 END DO
85 END DO
86 END DO
87
88 DO jk = 1, kflev
89 DO jl = 1, kdlon
90 ptauaz(jl, jk) = flag_aer*tauae(jl, jk, knu)
91 ppizaz(jl, jk) = flag_aer*pizae(jl, jk, knu)
92 pcgaz(jl, jk) = flag_aer*cgae(jl, jk, knu)
93 END DO
94
95 IF (flag_aer>0) THEN
96 ! -OB
97 DO jl = 1, kdlon
98 ! PCGAZ(JL,JK)=PCGAZ(JL,JK)/PPIZAZ(JL,JK)
99 ! PPIZAZ(JL,JK)=PPIZAZ(JL,JK)/PTAUAZ(JL,JK)
100 ztray = prayl(jl)*pdsig(jl, jk)
101 zratio = ztray/(ztray+ptauaz(jl,jk))
102 zgar = pcgaz(jl, jk)
103 zff = zgar*zgar
104 ptauaz(jl, jk) = ztray + ptauaz(jl, jk)*(1.-ppizaz(jl,jk)*zff)
105 pcgaz(jl, jk) = zgar*(1.-zratio)/(1.+zgar)
106 ppizaz(jl, jk) = zratio + (1.-zratio)*ppizaz(jl, jk)*(1.-zff)/(1.- &
107 ppizaz(jl,jk)*zff)
108 END DO
109 ELSE
110 DO jl = 1, kdlon
111 ztray = prayl(jl)*pdsig(jl, jk)
112 ptauaz(jl, jk) = ztray
113 pcgaz(jl, jk) = 0.
114 ppizaz(jl, jk) = 1. - repsct
115 END DO
116 END IF ! check flag_aer
117 END DO
118
119 ! ------------------------------------------------------------------
120
121 ! * 2. TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL
122 ! ----------------------------------------------
123
124
125 DO jl = 1, kdlon
126 zr23(jl) = 0.
127 zc0i(jl, kflev+1) = 0.
128 zclear(jl) = 1.
129 zscat(jl) = 0.
130 END DO
131
132 jk = 1
133 jkl = kflev + 1 - jk
134 jklp1 = jkl + 1
135 DO jl = 1, kdlon
136 zfacoa = 1. - ppizaz(jl, jkl)*pcgaz(jl, jkl)*pcgaz(jl, jkl)
137 zcorae = zfacoa*ptauaz(jl, jkl)*psec(jl)
138 zr21(jl) = exp(-zcorae)
139 zss0(jl) = 1. - zr21(jl)
140 zcle0(jl, jkl) = zss0(jl)
141
142 IF (novlp==1) THEN
143 ! * maximum-random
144 zclear(jl) = zclear(jl)*(1.0-max(zss0(jl),zscat(jl)))/ &
145 (1.0-min(zscat(jl),1.-zepsec))
146 zc0i(jl, jkl) = 1.0 - zclear(jl)
147 zscat(jl) = zss0(jl)
148 ELSE IF (novlp==2) THEN
149 ! * maximum
150 zscat(jl) = max(zss0(jl), zscat(jl))
151 zc0i(jl, jkl) = zscat(jl)
152 ELSE IF (novlp==3) THEN
153 ! * random
154 zclear(jl) = zclear(jl)*(1.0-zss0(jl))
155 zscat(jl) = 1.0 - zclear(jl)
156 zc0i(jl, jkl) = zscat(jl)
157 END IF
158 END DO
159
160 DO jk = 2, kflev
161 jkl = kflev + 1 - jk
162 jklp1 = jkl + 1
163 DO jl = 1, kdlon
164 zfacoa = 1. - ppizaz(jl, jkl)*pcgaz(jl, jkl)*pcgaz(jl, jkl)
165 zcorae = zfacoa*ptauaz(jl, jkl)*psec(jl)
166 zr21(jl) = exp(-zcorae)
167 zss0(jl) = 1. - zr21(jl)
168 zcle0(jl, jkl) = zss0(jl)
169
170 IF (novlp==1) THEN
171 ! * maximum-random
172 zclear(jl) = zclear(jl)*(1.0-max(zss0(jl),zscat(jl)))/ &
173 (1.0-min(zscat(jl),1.-zepsec))
174 zc0i(jl, jkl) = 1.0 - zclear(jl)
175 zscat(jl) = zss0(jl)
176 ELSE IF (novlp==2) THEN
177 ! * maximum
178 zscat(jl) = max(zss0(jl), zscat(jl))
179 zc0i(jl, jkl) = zscat(jl)
180 ELSE IF (novlp==3) THEN
181 ! * random
182 zclear(jl) = zclear(jl)*(1.0-zss0(jl))
183 zscat(jl) = 1.0 - zclear(jl)
184 zc0i(jl, jkl) = zscat(jl)
185 END IF
186 END DO
187 END DO
188
189 ! ------------------------------------------------------------------
190
191 ! * 3. REFLECTIVITY/TRANSMISSIVITY FOR PURE SCATTERING
192 ! -----------------------------------------------
193
194
195 DO jl = 1, kdlon
196 pray1(jl, kflev+1) = 0.
197 pray2(jl, kflev+1) = 0.
198 prefz(jl, 2, 1) = palbp(jl, knu)
199 prefz(jl, 1, 1) = palbp(jl, knu)
200 ptra1(jl, kflev+1) = 1.
201 ptra2(jl, kflev+1) = 1.
202 END DO
203
204 DO jk = 2, kflev + 1
205 jkm1 = jk - 1
206 DO jl = 1, kdlon
207
208
209 ! ------------------------------------------------------------------
210
211 ! * 3.1 EQUIVALENT ZENITH ANGLE
212 ! -----------------------
213
214
215 zmue = (1.-zc0i(jl,jk))*psec(jl) + zc0i(jl, jk)*1.66
216 prmu0(jl, jk) = 1./zmue
217
218
219 ! ------------------------------------------------------------------
220
221 ! * 3.2 REFLECT./TRANSMISSIVITY DUE TO RAYLEIGH AND AEROSOLS
222 ! ----------------------------------------------------
223
224
225 zgap = pcgaz(jl, jkm1)
226 zbmu0 = 0.5 - 0.75*zgap/zmue
227 zww = ppizaz(jl, jkm1)
228 zto = ptauaz(jl, jkm1)
229 zden = 1. + (1.-zww+zbmu0*zww)*zto*zmue + (1-zww)*(1.-zww+2.*zbmu0*zww) &
230 *zto*zto*zmue*zmue
231 pray1(jl, jkm1) = zbmu0*zww*zto*zmue/zden
232 ptra1(jl, jkm1) = 1./zden
233
234 zmu1 = 0.5
235 zbmu1 = 0.5 - 0.75*zgap*zmu1
236 zden1 = 1. + (1.-zww+zbmu1*zww)*zto/zmu1 + (1-zww)*(1.-zww+2.*zbmu1*zww &
237 )*zto*zto/zmu1/zmu1
238 pray2(jl, jkm1) = zbmu1*zww*zto/zmu1/zden1
239 ptra2(jl, jkm1) = 1./zden1
240
241
242
243 prefz(jl, 1, jk) = (pray1(jl,jkm1)+prefz(jl,1,jkm1)*ptra1(jl,jkm1)* &
244 ptra2(jl,jkm1)/(1.-pray2(jl,jkm1)*prefz(jl,1,jkm1)))
245
246 ztr(jl, 1, jkm1) = (ptra1(jl,jkm1)/(1.-pray2(jl,jkm1)*prefz(jl,1, &
247 jkm1)))
248
249 prefz(jl, 2, jk) = (pray1(jl,jkm1)+prefz(jl,2,jkm1)*ptra1(jl,jkm1)* &
250 ptra2(jl,jkm1))
251
252 ztr(jl, 2, jkm1) = ptra1(jl, jkm1)
253
254 END DO
255 END DO
256 DO jl = 1, kdlon
257 zmue = (1.-zc0i(jl,1))*psec(jl) + zc0i(jl, 1)*1.66
258 prmu0(jl, 1) = 1./zmue
259 END DO
260
261
262 ! ------------------------------------------------------------------
263
264 ! * 3.5 REFLECT./TRANSMISSIVITY BETWEEN SURFACE AND LEVEL
265 ! -------------------------------------------------
266
267
268 IF (knu==1) THEN
269 jaj = 2
270 DO jl = 1, kdlon
271 prj(jl, jaj, kflev+1) = 1.
272 prk(jl, jaj, kflev+1) = prefz(jl, 1, kflev+1)
273 END DO
274
275 DO jk = 1, kflev
276 jkl = kflev + 1 - jk
277 jklp1 = jkl + 1
278 DO jl = 1, kdlon
279 zre11 = prj(jl, jaj, jklp1)*ztr(jl, 1, jkl)
280 prj(jl, jaj, jkl) = zre11
281 prk(jl, jaj, jkl) = zre11*prefz(jl, 1, jkl)
282 END DO
283 END DO
284
285 ELSE
286
287 DO jaj = 1, 2
288 DO jl = 1, kdlon
289 prj(jl, jaj, kflev+1) = 1.
290 prk(jl, jaj, kflev+1) = prefz(jl, jaj, kflev+1)
291 END DO
292
293 DO jk = 1, kflev
294 jkl = kflev + 1 - jk
295 jklp1 = jkl + 1
296 DO jl = 1, kdlon
297 zre11 = prj(jl, jaj, jklp1)*ztr(jl, jaj, jkl)
298 prj(jl, jaj, jkl) = zre11
299 prk(jl, jaj, jkl) = zre11*prefz(jl, jaj, jkl)
300 END DO
301 END DO
302 END DO
303
304 END IF
305
306 ! ------------------------------------------------------------------
307
308 RETURN
309 END SUBROUTINE swclr

  ViewVC Help
Powered by ViewVC 1.1.21