/[lmdze]/trunk/Sources/phylmd/yamada4.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/yamada4.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 227 - (show annotations)
Thu Nov 2 15:47:03 2017 UTC (6 years, 6 months ago) by guez
File size: 12041 byte(s)
Rename phisinit to phis in restart.nc: clearer, same name as Fortran variable.

In aaam_bud, use rlat and rlon from phyetat0_m instead of having these
module variables associated to actual arguments in physiq.

In clmain, too many wind variables make the procedure hard to
understand. Use yu(:knon, 1) and yv(:knon, 1) instead of u1lay(:knon)
and v1lay(:knon). Note that when yu(:knon, 1) and yv(:knon, 1) are
used as actual arguments, they are probably copied to new arrays since
the elements are not contiguous. Rename yu10m to wind10m because this
is the norm of wind vector, not its zonal component. Rename yustar to
ustar. Rename uzon and vmer to u1 and v1 since these are wind
components at first layer and u1 and v1 are the names of corresponding
dummy arguments in stdlevvar.

In clmain, rename yzlev to zlev.

In clmain, screenc, stdlevvar and coefcdrag, remove the code
corresponding to zxli true (not used in LMDZ either).

Subroutine ustarhb becomes a function. Simplifications using the fact
that zx_alf2 = 0 and zx_alf1 = 1 (discarding the possibility to change
this).

In procedure vdif_kcay, remove unused dummy argument plev. Remove
useless computations of sss and sssq.

In clouds_gno, exp(100.) would overflow in single precision. Set
maximum to exp(80.) instead.

In physiq, use u(:, 1) and v(:, 1) as arguments to phytrac instead of
creating ad hoc variables yu1 and yv1.

In stdlevvar, rename dummy argument u_10m to wind10m, following the
corresponding modification in clmain. Simplifications using the fact
that ok_pred = 0 and ok_corr = 1 (discarding the possibility to change
this).

1 module yamada4_m
2
3 IMPLICIT NONE
4
5 private
6 public yamada4
7 real, parameter:: kap = 0.4
8
9 contains
10
11 SUBROUTINE yamada4(dt, g, zlev, zlay, u, v, teta, cd, q2, km, kn, kq, &
12 ustar, iflag_pbl)
13
14 ! From LMDZ4/libf/phylmd/yamada4.F, version 1.1 2004/06/22 11:45:36
15
16 use nr_util, only: assert, assert_eq
17 USE dimphy, ONLY: klev
18
19 REAL, intent(in):: dt ! pas de temps
20 real, intent(in):: g
21
22 REAL zlev(:, :) ! (knon, klev + 1)
23 ! altitude \`a chaque niveau (interface inf\'erieure de la couche de
24 ! m\^eme indice)
25
26 REAL, intent(in):: zlay(:, :) ! (knon, klev) altitude au centre de
27 ! chaque couche
28
29 REAL, intent(in):: u(:, :), v(:, :) ! (knon, klev)
30 ! vitesse au centre de chaque couche (en entr\'ee : la valeur au
31 ! d\'ebut du pas de temps)
32
33 REAL, intent(in):: teta(:, :) ! (knon, klev)
34 ! temp\'erature potentielle au centre de chaque couche (en entr\'ee :
35 ! la valeur au d\'ebut du pas de temps)
36
37 REAL, intent(in):: cd(:) ! (knon) cdrag, valeur au d\'ebut du pas de temps
38
39 REAL, intent(inout):: q2(:, :) ! (knon, klev + 1)
40 ! $q^2$ au bas de chaque couche
41 ! En entr\'ee : la valeur au d\'ebut du pas de temps ; en sortie : la
42 ! valeur \`a la fin du pas de temps.
43
44 REAL km(:, :) ! (knon, klev + 1)
45 ! diffusivit\'e turbulente de quantit\'e de mouvement (au bas de
46 ! chaque couche) (en sortie : la valeur \`a la fin du pas de temps)
47
48 REAL kn(:, :) ! (knon, klev + 1)
49 ! diffusivit\'e turbulente des scalaires (au bas de chaque couche)
50 ! (en sortie : la valeur \`a la fin du pas de temps)
51
52 REAL kq(:, :) ! (knon, klev + 1)
53 real, intent(in):: ustar(:) ! (knon)
54
55 integer, intent(in):: iflag_pbl
56 ! iflag_pbl doit valoir entre 6 et 9
57 ! l = 6, on prend syst\'ematiquement une longueur d'\'equilibre
58 ! iflag_pbl = 6 : MY 2.0
59 ! iflag_pbl = 7 : MY 2.0.Fournier
60 ! iflag_pbl = 8 : MY 2.5
61 ! iflag_pbl = 9 : un test ?
62
63 ! Local:
64 integer knon
65 real kmin, qmin
66 real pblhmin(size(cd)), coriol(size(cd)) ! (knon)
67 real qpre
68 REAL unsdz(size(zlay, 1), size(zlay, 2)) ! (knon, klev)
69 REAL unsdzdec(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
70 REAL kmpre(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
71 real tmp2
72 REAL mpre(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
73 real delta(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
74 real aa(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
75 real aa1
76 logical:: first = .true.
77 integer:: ipas = 0
78 integer ig, k
79 real ri
80 real, dimension(size(zlev, 1), size(zlev, 2)):: rif, sm ! (knon, klev + 1)
81 real alpha(size(zlay, 1), size(zlay, 2)) ! (knon, klev)
82
83 real, dimension(size(zlev, 1), size(zlev, 2)):: m2, dz, n2
84 ! (knon, klev + 1)
85
86 real zq
87 real dtetadz(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
88 real m2cstat, mcstat, kmcstat
89 real l(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
90 real l0(size(cd)) ! (knon)
91 real sq(size(cd)), sqz(size(cd)) ! (knon)
92 real zz(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
93 integer iter
94 real:: ric = 0.195, rifc = 0.191, b1 = 16.6
95
96 !-----------------------------------------------------------------------
97
98 call assert(iflag_pbl >= 6 .and. iflag_pbl <= 9, "yamada4 iflag_pbl")
99 knon = assert_eq([size(zlev, 1), size(zlay, 1), size(u, 1), size(v, 1), &
100 size(teta, 1), size(cd), size(q2, 1), size(km, 1), size(kn, 1), &
101 size(kq, 1)], "yamada4 knon")
102 call assert(klev == [size(zlev, 2) - 1, size(zlay, 2), size(u, 2), &
103 size(v, 2), size(teta, 2), size(q2, 2) - 1, size(km, 2) - 1, &
104 size(kn, 2) - 1, size(kq, 2) - 1], "yamada4 klev")
105
106 ipas = ipas + 1
107
108 ! les increments verticaux
109 DO ig = 1, knon
110 ! alerte: zlev n'est pas declare a klev + 1
111 zlev(ig, klev + 1) = zlay(ig, klev) + (zlay(ig, klev) - zlev(ig, klev))
112 ENDDO
113
114 DO k = 1, klev
115 DO ig = 1, knon
116 unsdz(ig, k) = 1.E+0/(zlev(ig, k + 1)-zlev(ig, k))
117 ENDDO
118 ENDDO
119
120 DO ig = 1, knon
121 unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1))
122 ENDDO
123
124 DO k = 2, klev
125 DO ig = 1, knon
126 unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1))
127 ENDDO
128 ENDDO
129
130 DO ig = 1, knon
131 unsdzdec(ig, klev + 1) = 1.E+0/(zlev(ig, klev + 1)-zlay(ig, klev))
132 ENDDO
133
134 do k = 2, klev
135 do ig = 1, knon
136 dz(ig, k) = zlay(ig, k)-zlay(ig, k-1)
137 m2(ig, k) = ((u(ig, k)-u(ig, k-1))**2 + (v(ig, k)-v(ig, k-1))**2) &
138 /(dz(ig, k)*dz(ig, k))
139 dtetadz(ig, k) = (teta(ig, k)-teta(ig, k-1))/dz(ig, k)
140 n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig, k-1) + teta(ig, k))
141 ri = n2(ig, k)/max(m2(ig, k), 1.e-10)
142 if (ri.lt.ric) then
143 rif(ig, k) = frif(ri)
144 else
145 rif(ig, k) = rifc
146 endif
147 if (rif(ig, k).lt.0.16) then
148 alpha(ig, k) = falpha(rif(ig, k))
149 sm(ig, k) = fsm(rif(ig, k))
150 else
151 alpha(ig, k) = 1.12
152 sm(ig, k) = 0.085
153 endif
154 zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig, k))*sm(ig, k)
155 enddo
156 enddo
157
158 ! Au premier appel, on d\'etermine l et q2 de fa\ccon it\'erative.
159 ! It\'eration pour d\'eterminer la longueur de m\'elange
160
161 if (first .or. iflag_pbl == 6) then
162 do ig = 1, knon
163 l0(ig) = 10.
164 enddo
165 do k = 2, klev-1
166 do ig = 1, knon
167 l(ig, k) = l0(ig) * kap * zlev(ig, k) &
168 / (kap * zlev(ig, k) + l0(ig))
169 enddo
170 enddo
171
172 do iter = 1, 10
173 do ig = 1, knon
174 sq(ig) = 1e-10
175 sqz(ig) = 1e-10
176 enddo
177 do k = 2, klev-1
178 do ig = 1, knon
179 q2(ig, k) = l(ig, k)**2 * zz(ig, k)
180 l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k))
181 zq = sqrt(q2(ig, k))
182 sqz(ig) = sqz(ig) + zq * zlev(ig, k) &
183 * (zlay(ig, k) - zlay(ig, k-1))
184 sq(ig) = sq(ig) + zq * (zlay(ig, k) - zlay(ig, k-1))
185 enddo
186 enddo
187 do ig = 1, knon
188 l0(ig) = 0.2 * sqz(ig) / sq(ig)
189 enddo
190 enddo
191 endif
192
193 ! Calcul de la longueur de melange.
194
195 ! Mise a jour de l0
196 do ig = 1, knon
197 sq(ig) = 1.e-10
198 sqz(ig) = 1.e-10
199 enddo
200 do k = 2, klev-1
201 do ig = 1, knon
202 zq = sqrt(q2(ig, k))
203 sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig, k)-zlay(ig, k-1))
204 sq(ig) = sq(ig) + zq*(zlay(ig, k)-zlay(ig, k-1))
205 enddo
206 enddo
207 do ig = 1, knon
208 l0(ig) = 0.2*sqz(ig)/sq(ig)
209 enddo
210 ! calcul de l(z)
211 do k = 2, klev
212 do ig = 1, knon
213 l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k))
214 if (first) then
215 q2(ig, k) = l(ig, k)**2 * zz(ig, k)
216 endif
217 enddo
218 enddo
219
220 ! Yamada 2.0
221 if (iflag_pbl == 6) then
222 do k = 2, klev
223 do ig = 1, knon
224 q2(ig, k) = l(ig, k)**2 * zz(ig, k)
225 enddo
226 enddo
227 else if (iflag_pbl == 7) then
228 ! Yamada 2.Fournier
229
230 ! Calcul de l, km, au pas precedent
231 do k = 2, klev
232 do ig = 1, knon
233 delta(ig, k) = q2(ig, k) / (l(ig, k)**2 * sm(ig, k))
234 kmpre(ig, k) = l(ig, k) * sqrt(q2(ig, k)) * sm(ig, k)
235 mpre(ig, k) = sqrt(m2(ig, k))
236 enddo
237 enddo
238
239 do k = 2, klev-1
240 do ig = 1, knon
241 m2cstat = max(alpha(ig, k)*n2(ig, k) + delta(ig, k)/b1, 1.e-12)
242 mcstat = sqrt(m2cstat)
243
244 ! puis on ecrit la valeur de q qui annule l'equation de m
245 ! supposee en q3
246
247 IF (k == 2) THEN
248 kmcstat = 1.E+0 / mcstat &
249 *(unsdz(ig, k)*kmpre(ig, k + 1) &
250 *mpre(ig, k + 1) &
251 + unsdz(ig, k-1) &
252 *cd(ig) &
253 *(sqrt(u(ig, 3)**2 + v(ig, 3)**2) &
254 -mcstat/unsdzdec(ig, k) &
255 -mpre(ig, k + 1)/unsdzdec(ig, k + 1))**2) &
256 /(unsdz(ig, k) + unsdz(ig, k-1))
257 ELSE
258 kmcstat = 1.E+0 / mcstat &
259 *(unsdz(ig, k)*kmpre(ig, k + 1) &
260 *mpre(ig, k + 1) &
261 + unsdz(ig, k-1)*kmpre(ig, k-1) &
262 *mpre(ig, k-1)) &
263 /(unsdz(ig, k) + unsdz(ig, k-1))
264 ENDIF
265 tmp2 = kmcstat / (sm(ig, k) / q2(ig, k)) /l(ig, k)
266 q2(ig, k) = max(tmp2, 1.e-12)**(2./3.)
267 enddo
268 enddo
269 else if (iflag_pbl >= 8) then
270 ! Yamada 2.5 a la Didi
271
272 ! Calcul de l, km, au pas precedent
273 do k = 2, klev
274 do ig = 1, knon
275 delta(ig, k) = q2(ig, k)/(l(ig, k)**2*sm(ig, k))
276 if (delta(ig, k).lt.1.e-20) then
277 delta(ig, k) = 1.e-20
278 endif
279 km(ig, k) = l(ig, k)*sqrt(q2(ig, k))*sm(ig, k)
280 aa1 = (m2(ig, k)*(1.-rif(ig, k))-delta(ig, k)/b1)
281 aa(ig, k) = aa1*dt/(delta(ig, k)*l(ig, k))
282 qpre = sqrt(q2(ig, k))
283 if (iflag_pbl == 8) then
284 if (aa(ig, k).gt.0.) then
285 q2(ig, k) = (qpre + aa(ig, k)*qpre*qpre)**2
286 else
287 q2(ig, k) = (qpre/(1.-aa(ig, k)*qpre))**2
288 endif
289 else
290 ! iflag_pbl = 9
291 if (aa(ig, k)*qpre.gt.0.9) then
292 q2(ig, k) = (qpre*10.)**2
293 else
294 q2(ig, k) = (qpre/(1.-aa(ig, k)*qpre))**2
295 endif
296 endif
297 q2(ig, k) = min(max(q2(ig, k), 1.e-10), 1.e4)
298 enddo
299 enddo
300 endif
301
302 ! Calcul des coefficients de m\'elange
303 do k = 2, klev
304 do ig = 1, knon
305 zq = sqrt(q2(ig, k))
306 km(ig, k) = l(ig, k)*zq*sm(ig, k)
307 kn(ig, k) = km(ig, k)*alpha(ig, k)
308 kq(ig, k) = l(ig, k)*zq*0.2
309 enddo
310 enddo
311
312 ! Traitement des cas noctrunes avec l'introduction d'une longueur
313 ! minilale.
314
315 ! Traitement particulier pour les cas tres stables.
316 ! D'apres Holtslag Boville.
317
318 do ig = 1, knon
319 coriol(ig) = 1.e-4
320 pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546e-5)
321 enddo
322
323 do k = 2, klev
324 do ig = 1, knon
325 if (teta(ig, 2).gt.teta(ig, 1)) then
326 qmin = ustar(ig)*(max(1.-zlev(ig, k)/pblhmin(ig), 0.))**2
327 kmin = kap*zlev(ig, k)*qmin
328 else
329 kmin = -1. ! kmin n'est utilise que pour les SL stables.
330 endif
331 if (kn(ig, k).lt.kmin.or.km(ig, k).lt.kmin) then
332 kn(ig, k) = kmin
333 km(ig, k) = kmin
334 kq(ig, k) = kmin
335 ! la longueur de melange est suposee etre l = kap z
336 ! K = l q Sm d'ou q2 = (K/l Sm)**2
337 q2(ig, k) = (qmin/sm(ig, k))**2
338 endif
339 enddo
340 enddo
341
342 first = .false.
343
344 end SUBROUTINE yamada4
345
346 !*******************************************************************
347
348 real function frif(ri)
349
350 real, intent(in):: ri
351
352 frif = 0.6588*(ri + 0.1776-sqrt(ri*ri-0.3221*ri + 0.03156))
353
354 end function frif
355
356 !*******************************************************************
357
358 real function falpha(ri)
359
360 real, intent(in):: ri
361
362 falpha = 1.318*(0.2231-ri)/(0.2341-ri)
363
364 end function falpha
365
366 !*******************************************************************
367
368 real function fsm(ri)
369
370 real, intent(in):: ri
371
372 fsm = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
373
374 end function fsm
375
376 !*******************************************************************
377
378 real function fl(zzz, zl0, zq2, zn2)
379
380 real, intent(in):: zzz, zl0, zq2, zn2
381
382 fl = max(min(zl0 * kap * zzz / (kap * zzz + zl0), &
383 0.5 * sqrt(zq2) / sqrt(max(zn2, 1e-10))), 1.)
384
385 end function fl
386
387 end module yamada4_m

  ViewVC Help
Powered by ViewVC 1.1.21