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

Contents of /trunk/Sources/phylmd/vdif_kcay.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: 17964 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 vdif_kcay_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE vdif_kcay(knon, dt, g, zlev, zlay, u, v, teta, cd, q2, q2diag, &
8 km, kn, ustar, l_mix)
9
10 ! From LMDZ4/libf/phylmd/vdif_kcay.F, version 1.1, 2004/06/22 11:45:36
11
12 USE dimphy, ONLY: klev, klon
13 use yamada_m, only: yamada
14
15 INTEGER knon
16 ! knon : nombre de points de grille
17
18 REAL, intent(in):: dt
19 ! dt : pas de temps
20 real, intent(in):: g
21 ! g : g
22 REAL zlev(klon, klev+1)
23 ! zlev : altitude a chaque niveau (interface inferieure de la couche
24 ! de meme indice)
25 REAL zlay(klon, klev)
26 ! zlay : altitude au centre de chaque couche
27 REAL, intent(in):: u(klon, klev)
28 REAL, intent(in):: v(klon, klev)
29 ! u, v : vitesse au centre de chaque couche
30 ! (en entree : la valeur au debut du pas de temps)
31 REAL teta(klon, klev)
32 ! teta : temperature potentielle au centre de chaque couche
33 ! (en entree : la valeur au debut du pas de temps)
34 REAL, intent(in):: cd (:) ! (knon) cdrag, valeur au debut du pas de temps
35 REAL q2(klon, klev+1)
36 ! q2 : $q^2$ au bas de chaque couche
37 ! (en entree : la valeur au debut du pas de temps)
38 ! (en sortie : la valeur a la fin du pas de temps)
39 REAL q2diag(klon, klev+1)
40 REAL km(klon, klev+1)
41 ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
42 ! couche)
43 ! (en sortie : la valeur a la fin du pas de temps)
44 REAL kn(klon, klev+1)
45 ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
46 ! (en sortie : la valeur a la fin du pas de temps)
47 real, intent(in):: ustar(:) ! (knon)
48 integer l_mix
49
50 ! Local:
51
52 real snstable
53 real sq(klon), sqz(klon), zq, long0(klon)
54
55 INTEGER nlay, nlev
56 ! nlay : nombre de couches
57 ! nlev : nombre de niveaux
58 REAL unsdz(klon, klev)
59 ! unsdz : 1 sur l'epaisseur de couche
60 REAL unsdzdec(klon, klev+1)
61 ! unsdzdec : 1 sur la distance entre le centre de la couche et le
62 ! centre de la couche inferieure
63 REAL q(klon, klev+1)
64 ! q : echelle de vitesse au bas de chaque couche
65 ! (valeur a la fin du pas de temps)
66
67 REAL kmpre(klon, klev+1)
68 ! kmpre : km au debut du pas de temps
69 REAL qcstat
70 ! qcstat : q : solution stationnaire du probleme couple
71 ! (valeur a la fin du pas de temps)
72 REAL q2cstat
73 ! q2cstat : q2 : solution stationnaire du probleme couple
74 ! (valeur a la fin du pas de temps)
75
76 REAL long(klon, klev+1)
77 ! long : longueur de melange calculee selon Blackadar
78
79 ! kmq3 : terme en q^3 dans le developpement de km
80 ! (valeur au debut du pas de temps)
81 ! kmcstat : valeur de km solution stationnaire du systeme {q2 ; du/dz}
82 ! (valeur a la fin du pas de temps)
83 ! knq3 : terme en q^3 dans le developpement de kn
84 ! mcstat : valeur de m solution stationnaire du systeme {q2 ; du/dz}
85 ! (valeur a la fin du pas de temps)
86 ! m2cstat : valeur de m2 solution stationnaire du systeme {q2 ; du/dz}
87 ! (valeur a la fin du pas de temps)
88 ! m : valeur a la fin du pas de temps
89 ! mpre : valeur au debut du pas de temps
90 ! m2 : valeur a la fin du pas de temps
91 ! n2 : valeur a la fin du pas de temps
92
93 REAL kmq3
94 REAL kmcstat
95 REAL knq3
96 REAL mcstat
97 REAL m2cstat
98 REAL m(klon, klev+1)
99 REAL mpre(klon, klev+1)
100 REAL m2(klon, klev+1)
101 REAL n2(klon, klev+1)
102
103 ! gn : intermediaire pour les coefficients de stabilite
104 ! gnmin : borne inferieure de gn (-0.23 ou -0.28)
105 ! gnmax : borne superieure de gn (0.0233)
106 ! gninf : vrai si gn est en dessous de sa borne inferieure
107 ! gnsup : vrai si gn est en dessus de sa borne superieure
108 ! ri : nombre de Richardson
109 ! sn : coefficient de stabilite pour n
110 ! snq2 : premier terme du developement limite de sn en q2
111 ! sm : coefficient de stabilite pour m
112 ! smq2 : premier terme du developement limite de sm en q2
113
114 REAL gn
115 REAL gnmin
116 REAL gnmax
117 LOGICAL gninf
118 LOGICAL gnsup
119 REAL sn(klon, klev+1)
120 REAL snq2(klon, klev+1)
121 REAL sm(klon, klev+1)
122 REAL smq2(klon, klev+1)
123
124 ! kappa : consatnte de Von Karman (0.4)
125 ! long00 : longueur de reference pour le calcul de long (160)
126 ! a1, a2, b1, b2, c1 : constantes d'origine pour les coefficients
127 ! de stabilite (0.92/0.74/16.6/10.1/0.08)
128 ! cn1, cn2 : constantes pour sn
129 ! cm1, cm2, cm3, cm4 : constantes pour sm
130
131 REAL kappa
132 REAL long00
133 REAL a1, a2, b1, b2, c1
134 REAL cn1, cn2
135 REAL cm1, cm2, cm3, cm4
136
137 ! termq : termes en $q$ dans l'equation de q2
138 ! termq3 : termes en $q^3$ dans l'equation de q2
139 ! termqm2 : termes en $q*m^2$ dans l'equation de q2
140 ! termq3m2 : termes en $q^3*m^2$ dans l'equation de q2
141
142 REAL termq
143 REAL termq3
144 REAL termqm2
145 REAL termq3m2
146
147 ! q2min : borne inferieure de q2
148 REAL q2min
149
150 ! knmin : borne inferieure de kn
151 ! kmmin : borne inferieure de km
152
153 REAL knmin
154 REAL kmmin
155
156 INTEGER ilay, ilev, igrid
157 REAL tmp1, tmp2
158
159 PARAMETER (kappa=0.4E+0)
160 PARAMETER (long00=160.E+0)
161 ! PARAMETER (gnmin=-10.E+0)
162 PARAMETER (gnmin=-0.28)
163 PARAMETER (gnmax=0.0233E+0)
164 PARAMETER (a1=0.92E+0)
165 PARAMETER (a2=0.74E+0)
166 PARAMETER (b1=16.6E+0)
167 PARAMETER (b2=10.1E+0)
168 PARAMETER (c1=0.08E+0)
169 PARAMETER (knmin=1.E-5)
170 PARAMETER (kmmin=1.E-5)
171 PARAMETER (q2min=1.e-5)
172 PARAMETER (nlay=klev)
173 PARAMETER (nlev=klev+1)
174
175 PARAMETER (cn1=a2*(1.E+0 -6.E+0 *a1/b1))
176 PARAMETER (cn2=-3.E+0 *a2*(6.E+0 *a1+b2))
177 PARAMETER (cm1=a1*(1.E+0 -3.E+0 *c1-6.E+0 *a1/b1))
178 PARAMETER (cm2=a1*(-3.E+0 *a2*((b2-3.E+0 *a2)*(1.E+0 -6.E+0 *a1/b1) &
179 -3.E+0 *c1*(b2+6.E+0 *a1))))
180 PARAMETER (cm3=-3.E+0 *a2*(6.E+0 *a1+b2))
181 PARAMETER (cm4=-9.E+0 *a1*a2)
182
183 logical:: first = .true.
184
185 !------------------------------------------------------------
186
187 ! traitment des valeur de q2 en entree
188
189 ! Initialisation de q2
190
191 call yamada(knon, g, zlev, zlay, u, v, teta, q2diag, km, kn)
192 if (first.and.1.eq.1) then
193 first=.false.
194 q2=q2diag
195 endif
196
197 DO ilev=1, nlev
198 DO igrid=1, knon
199 q2(igrid, ilev)=amax1(q2(igrid, ilev), q2min)
200 q(igrid, ilev)=sqrt(q2(igrid, ilev))
201 ENDDO
202 ENDDO
203
204 DO igrid=1, knon
205 tmp1=cd(igrid)*(u(igrid, 1)**2+v(igrid, 1)**2)
206 q2(igrid, 1)=b1**(2.E+0/3.E+0)*tmp1
207 q2(igrid, 1)=amax1(q2(igrid, 1), q2min)
208 q(igrid, 1)=sqrt(q2(igrid, 1))
209 ENDDO
210
211 ! les increments verticaux
212
213 ! allerte !c
214 ! zlev n'est pas declare a nlev !c
215 DO igrid=1, knon
216 zlev(igrid, nlev)=zlay(igrid, nlay) &
217 +( zlay(igrid, nlay) - zlev(igrid, nlev-1) )
218 ENDDO
219 ! allerte !c
220
221 DO ilay=1, nlay
222 DO igrid=1, knon
223 unsdz(igrid, ilay)=1.E+0/(zlev(igrid, ilay+1)-zlev(igrid, ilay))
224 ENDDO
225 ENDDO
226 DO igrid=1, knon
227 unsdzdec(igrid, 1)=1.E+0/(zlay(igrid, 1)-zlev(igrid, 1))
228 ENDDO
229 DO ilay=2, nlay
230 DO igrid=1, knon
231 unsdzdec(igrid, ilay)=1.E+0/(zlay(igrid, ilay)-zlay(igrid, ilay-1))
232 ENDDO
233 ENDDO
234 DO igrid=1, knon
235 unsdzdec(igrid, nlay+1)=1.E+0/(zlev(igrid, nlay+1)-zlay(igrid, nlay))
236 ENDDO
237
238 ! le cisaillement et le gradient de temperature
239
240 DO igrid=1, knon
241 m2(igrid, 1)=(unsdzdec(igrid, 1) &
242 *u(igrid, 1))**2 &
243 +(unsdzdec(igrid, 1) &
244 *v(igrid, 1))**2
245 m(igrid, 1)=sqrt(m2(igrid, 1))
246 mpre(igrid, 1)=m(igrid, 1)
247 ENDDO
248
249 DO ilev=2, nlev-1
250 DO igrid=1, knon
251
252 n2(igrid, ilev)=g*unsdzdec(igrid, ilev) &
253 *(teta(igrid, ilev)-teta(igrid, ilev-1)) &
254 /(teta(igrid, ilev)+teta(igrid, ilev-1)) *2.E+0
255
256 ! on ne sais traiter que les cas stratifies. et l'ajustement
257 ! convectif est cense faire en sorte que seul des
258 ! configurations stratifiees soient rencontrees en entree de
259 ! cette routine. mais, bon ... on sait jamais (meme on sait
260 ! que n2 prends quelques valeurs negatives ... parfois)
261 ! alors :
262
263 IF (n2(igrid, ilev).lt.0.E+0) THEN
264 n2(igrid, ilev)=0.E+0
265 ENDIF
266
267 m2(igrid, ilev)=(unsdzdec(igrid, ilev) &
268 *(u(igrid, ilev)-u(igrid, ilev-1)))**2 &
269 +(unsdzdec(igrid, ilev) &
270 *(v(igrid, ilev)-v(igrid, ilev-1)))**2
271 m(igrid, ilev)=sqrt(m2(igrid, ilev))
272 mpre(igrid, ilev)=m(igrid, ilev)
273
274 ENDDO
275 ENDDO
276
277 DO igrid=1, knon
278 m2(igrid, nlev)=m2(igrid, nlev-1)
279 m(igrid, nlev)=m(igrid, nlev-1)
280 mpre(igrid, nlev)=m(igrid, nlev)
281 ENDDO
282
283 ! calcul des fonctions de stabilite
284
285 if (l_mix.eq.4) then
286 DO igrid=1, knon
287 sqz(igrid)=1.e-10
288 sq(igrid)=1.e-10
289 ENDDO
290 do ilev=2, nlev-1
291 DO igrid=1, knon
292 zq=sqrt(q2(igrid, ilev))
293 sqz(igrid) &
294 =sqz(igrid)+zq*zlev(igrid, ilev) &
295 *(zlay(igrid, ilev)-zlay(igrid, ilev-1))
296 sq(igrid)=sq(igrid)+zq*(zlay(igrid, ilev)-zlay(igrid, ilev-1))
297 ENDDO
298 enddo
299 DO igrid=1, knon
300 long0(igrid)=0.2*sqz(igrid)/sq(igrid)
301 ENDDO
302 else if (l_mix.eq.3) then
303 long0(igrid)=long00
304 endif
305
306 DO ilev=2, nlev-1
307 DO igrid=1, knon
308 tmp1=kappa*(zlev(igrid, ilev)-zlev(igrid, 1))
309 if (l_mix.ge.10) then
310 long(igrid, ilev)=l_mix
311 else
312 long(igrid, ilev)=tmp1/(1.E+0 + tmp1/long0(igrid))
313 endif
314 long(igrid, ilev)=max(min(long(igrid, ilev) &
315 , 0.5*sqrt(q2(igrid, ilev))/sqrt(max(n2(igrid, ilev), 1.e-10))) &
316 , 5.)
317
318 gn=-long(igrid, ilev)**2 / q2(igrid, ilev) &
319 * n2(igrid, ilev)
320 gninf=.false.
321 gnsup=.false.
322 long(igrid, ilev)=long(igrid, ilev)
323 long(igrid, ilev)=long(igrid, ilev)
324
325 IF (gn.lt.gnmin) THEN
326 gninf=.true.
327 gn=gnmin
328 ENDIF
329
330 IF (gn.gt.gnmax) THEN
331 gnsup=.true.
332 gn=gnmax
333 ENDIF
334
335 sn(igrid, ilev)=cn1/(1.E+0 +cn2*gn)
336 sm(igrid, ilev)= &
337 (cm1+cm2*gn) &
338 /( (1.E+0 +cm3*gn) &
339 *(1.E+0 +cm4*gn) )
340
341 IF ((gninf).or.(gnsup)) THEN
342 snq2(igrid, ilev)=0.E+0
343 smq2(igrid, ilev)=0.E+0
344 ELSE
345 snq2(igrid, ilev)= &
346 -gn &
347 *(-cn1*cn2/(1.E+0 +cn2*gn)**2 )
348 smq2(igrid, ilev)= &
349 -gn &
350 *( cm2*(1.E+0 +cm3*gn) &
351 *(1.E+0 +cm4*gn) &
352 -( cm3*(1.E+0 +cm4*gn) &
353 +cm4*(1.E+0 +cm3*gn) ) &
354 *(cm1+cm2*gn) ) &
355 /( (1.E+0 +cm3*gn) &
356 *(1.E+0 +cm4*gn) )**2
357 ENDIF
358 ! la decomposition de Taylor en q2 n'a de sens que dans les
359 ! cas stratifies ou sn et sm sont quasi proportionnels a
360 ! q2. ailleurs on laisse le meme algorithme car l'ajustement
361 ! convectif fait le travail. mais c'est delirant quand sn
362 ! et snq2 n'ont pas le meme signe : dans ces cas, on ne fait
363 ! pas la decomposition.
364
365 IF (snq2(igrid, ilev)*sn(igrid, ilev).le.0.E+0) &
366 snq2(igrid, ilev)=0.E+0
367 IF (smq2(igrid, ilev)*sm(igrid, ilev).le.0.E+0) &
368 smq2(igrid, ilev)=0.E+0
369
370 ! Correction pour les couches stables.
371 ! Schema repris de JHoltzlag Boville, lui meme venant de...
372
373 snstable=1.-zlev(igrid, ilev) &
374 /(700.*max(ustar(igrid), 0.0001))
375 snstable=1.-zlev(igrid, ilev)/400.
376 snstable=max(snstable, 0.)
377 snstable=snstable*snstable
378
379 if (sn(igrid, ilev).lt.snstable) then
380 sn(igrid, ilev)=snstable
381 snq2(igrid, ilev)=0.
382 endif
383
384 if (sm(igrid, ilev).lt.snstable) then
385 sm(igrid, ilev)=snstable
386 smq2(igrid, ilev)=0.
387 endif
388
389 ! sn : coefficient de stabilite pour n
390 ! snq2 : premier terme du developement limite de sn en q2
391 ENDDO
392 ENDDO
393
394 ! calcul de km et kn au debut du pas de temps
395
396 DO igrid=1, knon
397 kn(igrid, 1)=knmin
398 km(igrid, 1)=kmmin
399 kmpre(igrid, 1)=km(igrid, 1)
400 ENDDO
401
402 DO ilev=2, nlev-1
403 DO igrid=1, knon
404 kn(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev) &
405 *sn(igrid, ilev)
406 km(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev) &
407 *sm(igrid, ilev)
408 kmpre(igrid, ilev)=km(igrid, ilev)
409 ENDDO
410 ENDDO
411
412 DO igrid=1, knon
413 kn(igrid, nlev)=kn(igrid, nlev-1)
414 km(igrid, nlev)=km(igrid, nlev-1)
415 kmpre(igrid, nlev)=km(igrid, nlev)
416 ENDDO
417
418 ! boucle sur les niveaux 2 a nlev-1
419
420 DO ilev=2, nlev-1
421 DO igrid=1, knon
422 ! calcul des termes sources et puits de l'equation de q2
423
424 knq3=kn(igrid, ilev)*snq2(igrid, ilev) &
425 /sn(igrid, ilev)
426 kmq3=km(igrid, ilev)*smq2(igrid, ilev) &
427 /sm(igrid, ilev)
428
429 termq=0.E+0
430 termq3=0.E+0
431 termqm2=0.E+0
432 termq3m2=0.E+0
433
434 tmp1=dt*2.E+0 *km(igrid, ilev)*m2(igrid, ilev)
435 tmp2=dt*2.E+0 *kmq3*m2(igrid, ilev)
436 termqm2=termqm2 &
437 +dt*2.E+0 *km(igrid, ilev)*m2(igrid, ilev) &
438 -dt*2.E+0 *kmq3*m2(igrid, ilev)
439 termq3m2=termq3m2 &
440 +dt*2.E+0 *kmq3*m2(igrid, ilev)
441
442 termq=termq &
443 -dt*2.E+0 *kn(igrid, ilev)*n2(igrid, ilev) &
444 +dt*2.E+0 *knq3*n2(igrid, ilev)
445 termq3=termq3 &
446 -dt*2.E+0 *knq3*n2(igrid, ilev)
447
448 termq3=termq3 &
449 -dt*2.E+0 *q(igrid, ilev)**3 / (b1*long(igrid, ilev))
450
451 ! resolution stationnaire couplee avec le gradient de vitesse local
452
453 ! on cherche le cisaillement qui annule l'equation de q^2
454 ! supposee en q3
455
456 tmp1=termq+termq3
457 tmp2=termqm2+termq3m2
458 m2cstat=m2(igrid, ilev) &
459 -(tmp1+tmp2)/(dt*2.E+0*km(igrid, ilev))
460 mcstat=sqrt(m2cstat)
461
462 ! puis on ecrit la valeur de q qui annule l'equation de m
463 ! supposee en q3
464
465 IF (ilev.eq.2) THEN
466 kmcstat=1.E+0 / mcstat &
467 *( unsdz(igrid, ilev)*kmpre(igrid, ilev+1) &
468 *mpre(igrid, ilev+1) &
469 +unsdz(igrid, ilev-1) &
470 *cd(igrid) &
471 *( sqrt(u(igrid, 3)**2+v(igrid, 3)**2) &
472 -mcstat/unsdzdec(igrid, ilev) &
473 -mpre(igrid, ilev+1)/unsdzdec(igrid, ilev+1) )**2) &
474 /( unsdz(igrid, ilev)+unsdz(igrid, ilev-1) )
475 ELSE
476 kmcstat=1.E+0 / mcstat &
477 *( unsdz(igrid, ilev)*kmpre(igrid, ilev+1) &
478 *mpre(igrid, ilev+1) &
479 +unsdz(igrid, ilev-1)*kmpre(igrid, ilev-1) &
480 *mpre(igrid, ilev-1) ) &
481 /( unsdz(igrid, ilev)+unsdz(igrid, ilev-1) )
482 ENDIF
483 tmp2=kmcstat &
484 /( sm(igrid, ilev)/q2(igrid, ilev) ) &
485 /long(igrid, ilev)
486 qcstat=tmp2**(1.E+0/3.E+0)
487 q2cstat=qcstat**2
488
489 ! choix de la solution finale
490
491 q(igrid, ilev)=qcstat
492 q2(igrid, ilev)=q2cstat
493 m(igrid, ilev)=mcstat
494 m2(igrid, ilev)=m2cstat
495
496 ! pour des raisons simples q2 est minore
497
498 IF (q2(igrid, ilev).lt.q2min) THEN
499 q2(igrid, ilev)=q2min
500 q(igrid, ilev)=sqrt(q2min)
501 ENDIF
502
503 ! calcul final de kn et km
504
505 gn=-long(igrid, ilev)**2 / q2(igrid, ilev) &
506 * n2(igrid, ilev)
507 IF (gn.lt.gnmin) gn=gnmin
508 IF (gn.gt.gnmax) gn=gnmax
509 sn(igrid, ilev)=cn1/(1.E+0 +cn2*gn)
510 sm(igrid, ilev)= &
511 (cm1+cm2*gn) &
512 /( (1.E+0 +cm3*gn)*(1.E+0 +cm4*gn) )
513 kn(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev) &
514 *sn(igrid, ilev)
515 km(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev) &
516 *sm(igrid, ilev)
517 end DO
518 end DO
519
520 DO igrid=1, knon
521 kn(igrid, 1)=knmin
522 km(igrid, 1)=kmmin
523 q2(igrid, nlev)=q2(igrid, nlev-1)
524 q(igrid, nlev)=q(igrid, nlev-1)
525 kn(igrid, nlev)=kn(igrid, nlev-1)
526 km(igrid, nlev)=km(igrid, nlev-1)
527 ENDDO
528
529 ! CALCUL DE LA DIFFUSION VERTICALE DE Q2
530 do ilev=1, nlev
531 do igrid=1, knon
532 q2(igrid, ilev)=max(q2(igrid, ilev), q2min)
533 q(igrid, ilev)=sqrt(q2(igrid, ilev))
534
535 ! calcul final de kn et km
536
537 gn=-long(igrid, ilev)**2 / q2(igrid, ilev) &
538 * n2(igrid, ilev)
539 IF (gn.lt.gnmin) gn=gnmin
540 IF (gn.gt.gnmax) gn=gnmax
541 sn(igrid, ilev)=cn1/(1.E+0 +cn2*gn)
542 sm(igrid, ilev)= &
543 (cm1+cm2*gn) &
544 /( (1.E+0 +cm3*gn)*(1.E+0 +cm4*gn) )
545 ! Correction pour les couches stables.
546 ! Schema repris de JHoltzlag Boville, lui meme venant de...
547
548 snstable=1.-zlev(igrid, ilev) &
549 /(700.*max(ustar(igrid), 0.0001))
550 snstable=1.-zlev(igrid, ilev)/400.
551 snstable=max(snstable, 0.)
552 snstable=snstable*snstable
553
554 if (sn(igrid, ilev).lt.snstable) then
555 sn(igrid, ilev)=snstable
556 snq2(igrid, ilev)=0.
557 endif
558
559 if (sm(igrid, ilev).lt.snstable) then
560 sm(igrid, ilev)=snstable
561 smq2(igrid, ilev)=0.
562 endif
563
564 ! sn : coefficient de stabilite pour n
565 kn(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev) &
566 *sn(igrid, ilev)
567 km(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev)
568 enddo
569 enddo
570
571 END SUBROUTINE vdif_kcay
572
573 end module vdif_kcay_m

  ViewVC Help
Powered by ViewVC 1.1.21