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

Annotation of /trunk/phylmd/yamada4.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 227 - (hide annotations)
Thu Nov 2 15:47:03 2017 UTC (6 years, 6 months ago) by guez
Original Path: trunk/Sources/phylmd/yamada4.f
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 guez 47 module yamada4_m
2 guez 3
3 guez 47 IMPLICIT NONE
4 guez 3
5 guez 118 private
6     public yamada4
7     real, parameter:: kap = 0.4
8 guez 3
9 guez 47 contains
10 guez 3
11 guez 227 SUBROUTINE yamada4(dt, g, zlev, zlay, u, v, teta, cd, q2, km, kn, kq, &
12 guez 118 ustar, iflag_pbl)
13 guez 3
14 guez 47 ! From LMDZ4/libf/phylmd/yamada4.F, version 1.1 2004/06/22 11:45:36
15 guez 3
16 guez 227 use nr_util, only: assert, assert_eq
17 guez 118 USE dimphy, ONLY: klev
18 guez 3
19 guez 47 REAL, intent(in):: dt ! pas de temps
20     real, intent(in):: g
21 guez 3
22 guez 227 REAL zlev(:, :) ! (knon, klev + 1)
23 guez 178 ! altitude \`a chaque niveau (interface inf\'erieure de la couche de
24     ! m\^eme indice)
25 guez 3
26 guez 227 REAL, intent(in):: zlay(:, :) ! (knon, klev) altitude au centre de
27     ! chaque couche
28 guez 3
29 guez 227 REAL, intent(in):: u(:, :), v(:, :) ! (knon, klev)
30 guez 178 ! vitesse au centre de chaque couche (en entr\'ee : la valeur au
31     ! d\'ebut du pas de temps)
32 guez 3
33 guez 227 REAL, intent(in):: teta(:, :) ! (knon, klev)
34 guez 178 ! temp\'erature potentielle au centre de chaque couche (en entr\'ee :
35     ! la valeur au d\'ebut du pas de temps)
36 guez 3
37 guez 227 REAL, intent(in):: cd(:) ! (knon) cdrag, valeur au d\'ebut du pas de temps
38 guez 3
39 guez 227 REAL, intent(inout):: q2(:, :) ! (knon, klev + 1)
40 guez 47 ! $q^2$ au bas de chaque couche
41 guez 178 ! 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 guez 3
44 guez 227 REAL km(:, :) ! (knon, klev + 1)
45 guez 178 ! 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 guez 3
48 guez 227 REAL kn(:, :) ! (knon, klev + 1)
49 guez 178 ! diffusivit\'e turbulente des scalaires (au bas de chaque couche)
50     ! (en sortie : la valeur \`a la fin du pas de temps)
51 guez 3
52 guez 227 REAL kq(:, :) ! (knon, klev + 1)
53     real, intent(in):: ustar(:) ! (knon)
54 guez 3
55 guez 119 integer, intent(in):: iflag_pbl
56 guez 47 ! iflag_pbl doit valoir entre 6 et 9
57 guez 178 ! l = 6, on prend syst\'ematiquement une longueur d'\'equilibre
58 guez 47 ! 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 guez 3
63 guez 47 ! Local:
64 guez 227 integer knon
65     real kmin, qmin
66     real pblhmin(size(cd)), coriol(size(cd)) ! (knon)
67 guez 47 real qpre
68 guez 227 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 guez 47 logical:: first = .true.
77     integer:: ipas = 0
78     integer ig, k
79     real ri
80 guez 227 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 guez 47 real m2cstat, mcstat, kmcstat
89 guez 227 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 guez 47 integer iter
94 guez 196 real:: ric = 0.195, rifc = 0.191, b1 = 16.6
95 guez 3
96 guez 47 !-----------------------------------------------------------------------
97 guez 3
98 guez 227 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 guez 3
106 guez 227 ipas = ipas + 1
107 guez 3
108 guez 47 ! les increments verticaux
109 guez 227 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 guez 47 ENDDO
113 guez 3
114 guez 118 DO k = 1, klev
115 guez 227 DO ig = 1, knon
116     unsdz(ig, k) = 1.E+0/(zlev(ig, k + 1)-zlev(ig, k))
117 guez 47 ENDDO
118     ENDDO
119 guez 118
120 guez 227 DO ig = 1, knon
121 guez 47 unsdzdec(ig, 1) = 1.E+0/(zlay(ig, 1)-zlev(ig, 1))
122     ENDDO
123 guez 118
124     DO k = 2, klev
125 guez 227 DO ig = 1, knon
126 guez 47 unsdzdec(ig, k) = 1.E+0/(zlay(ig, k)-zlay(ig, k-1))
127     ENDDO
128     ENDDO
129 guez 118
130 guez 227 DO ig = 1, knon
131     unsdzdec(ig, klev + 1) = 1.E+0/(zlev(ig, klev + 1)-zlay(ig, klev))
132 guez 47 ENDDO
133 guez 3
134 guez 47 do k = 2, klev
135 guez 227 do ig = 1, knon
136 guez 47 dz(ig, k) = zlay(ig, k)-zlay(ig, k-1)
137 guez 227 m2(ig, k) = ((u(ig, k)-u(ig, k-1))**2 + (v(ig, k)-v(ig, k-1))**2) &
138 guez 47 /(dz(ig, k)*dz(ig, k))
139     dtetadz(ig, k) = (teta(ig, k)-teta(ig, k-1))/dz(ig, k)
140 guez 227 n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig, k-1) + teta(ig, k))
141 guez 47 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 guez 119 if (rif(ig, k).lt.0.16) then
148 guez 47 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 guez 3
158 guez 196 ! Au premier appel, on d\'etermine l et q2 de fa\ccon it\'erative.
159 guez 178 ! It\'eration pour d\'eterminer la longueur de m\'elange
160 guez 3
161 guez 47 if (first .or. iflag_pbl == 6) then
162 guez 227 do ig = 1, knon
163 guez 47 l0(ig) = 10.
164     enddo
165     do k = 2, klev-1
166 guez 227 do ig = 1, knon
167 guez 47 l(ig, k) = l0(ig) * kap * zlev(ig, k) &
168     / (kap * zlev(ig, k) + l0(ig))
169     enddo
170     enddo
171 guez 3
172 guez 47 do iter = 1, 10
173 guez 227 do ig = 1, knon
174 guez 47 sq(ig) = 1e-10
175     sqz(ig) = 1e-10
176     enddo
177     do k = 2, klev-1
178 guez 227 do ig = 1, knon
179 guez 47 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 guez 227 do ig = 1, knon
188 guez 47 l0(ig) = 0.2 * sqz(ig) / sq(ig)
189     enddo
190     enddo
191     endif
192 guez 3
193 guez 47 ! Calcul de la longueur de melange.
194 guez 3
195 guez 47 ! Mise a jour de l0
196 guez 227 do ig = 1, knon
197 guez 47 sq(ig) = 1.e-10
198     sqz(ig) = 1.e-10
199     enddo
200     do k = 2, klev-1
201 guez 227 do ig = 1, knon
202 guez 47 zq = sqrt(q2(ig, k))
203 guez 227 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 guez 47 enddo
206     enddo
207 guez 227 do ig = 1, knon
208 guez 47 l0(ig) = 0.2*sqz(ig)/sq(ig)
209     enddo
210     ! calcul de l(z)
211     do k = 2, klev
212 guez 227 do ig = 1, knon
213 guez 47 l(ig, k) = fl(zlev(ig, k), l0(ig), q2(ig, k), n2(ig, k))
214 guez 119 if (first) then
215 guez 47 q2(ig, k) = l(ig, k)**2 * zz(ig, k)
216     endif
217     enddo
218     enddo
219 guez 3
220 guez 47 ! Yamada 2.0
221     if (iflag_pbl == 6) then
222     do k = 2, klev
223 guez 227 do ig = 1, knon
224 guez 47 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 guez 3
230 guez 47 ! Calcul de l, km, au pas precedent
231     do k = 2, klev
232 guez 227 do ig = 1, knon
233 guez 47 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 guez 3
239 guez 47 do k = 2, klev-1
240 guez 227 do ig = 1, knon
241     m2cstat = max(alpha(ig, k)*n2(ig, k) + delta(ig, k)/b1, 1.e-12)
242 guez 47 mcstat = sqrt(m2cstat)
243 guez 3
244 guez 47 ! puis on ecrit la valeur de q qui annule l'equation de m
245     ! supposee en q3
246 guez 3
247 guez 47 IF (k == 2) THEN
248     kmcstat = 1.E+0 / mcstat &
249 guez 227 *(unsdz(ig, k)*kmpre(ig, k + 1) &
250     *mpre(ig, k + 1) &
251     + unsdz(ig, k-1) &
252 guez 47 *cd(ig) &
253 guez 227 *(sqrt(u(ig, 3)**2 + v(ig, 3)**2) &
254 guez 47 -mcstat/unsdzdec(ig, k) &
255 guez 227 -mpre(ig, k + 1)/unsdzdec(ig, k + 1))**2) &
256     /(unsdz(ig, k) + unsdz(ig, k-1))
257 guez 47 ELSE
258     kmcstat = 1.E+0 / mcstat &
259 guez 227 *(unsdz(ig, k)*kmpre(ig, k + 1) &
260     *mpre(ig, k + 1) &
261     + unsdz(ig, k-1)*kmpre(ig, k-1) &
262 guez 47 *mpre(ig, k-1)) &
263 guez 227 /(unsdz(ig, k) + unsdz(ig, k-1))
264 guez 47 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 guez 3
272 guez 47 ! Calcul de l, km, au pas precedent
273     do k = 2, klev
274 guez 227 do ig = 1, knon
275 guez 47 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 guez 227 q2(ig, k) = (qpre + aa(ig, k)*qpre*qpre)**2
286 guez 47 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 guez 3
302 guez 178 ! Calcul des coefficients de m\'elange
303 guez 47 do k = 2, klev
304 guez 227 do ig = 1, knon
305 guez 47 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 guez 3
312 guez 47 ! Traitement des cas noctrunes avec l'introduction d'une longueur
313     ! minilale.
314 guez 3
315 guez 47 ! Traitement particulier pour les cas tres stables.
316     ! D'apres Holtslag Boville.
317 guez 3
318 guez 227 do ig = 1, knon
319 guez 47 coriol(ig) = 1.e-4
320     pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546e-5)
321     enddo
322 guez 3
323 guez 47 do k = 2, klev
324 guez 227 do ig = 1, knon
325 guez 47 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 guez 3
342 guez 47 first = .false.
343 guez 3
344 guez 47 end SUBROUTINE yamada4
345 guez 3
346 guez 47 !*******************************************************************
347 guez 3
348 guez 62 real function frif(ri)
349 guez 3
350 guez 47 real, intent(in):: ri
351 guez 3
352 guez 227 frif = 0.6588*(ri + 0.1776-sqrt(ri*ri-0.3221*ri + 0.03156))
353 guez 3
354 guez 47 end function frif
355 guez 3
356 guez 47 !*******************************************************************
357 guez 3
358 guez 62 real function falpha(ri)
359 guez 3
360 guez 47 real, intent(in):: ri
361 guez 3
362 guez 47 falpha = 1.318*(0.2231-ri)/(0.2341-ri)
363 guez 3
364 guez 47 end function falpha
365    
366     !*******************************************************************
367    
368 guez 62 real function fsm(ri)
369 guez 47
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 guez 62 real function fl(zzz, zl0, zq2, zn2)
379 guez 47
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