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

Annotation of /trunk/Sources/phylmd/vdif_kcay.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 227 - (hide annotations)
Thu Nov 2 15:47:03 2017 UTC (6 years, 7 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 guez 62 module vdif_kcay_m
2 guez 3
3 guez 62 IMPLICIT NONE
4 guez 3
5 guez 62 contains
6 guez 3
7 guez 227 SUBROUTINE vdif_kcay(knon, dt, g, zlev, zlay, u, v, teta, cd, q2, q2diag, &
8     km, kn, ustar, l_mix)
9 guez 3
10 guez 227 ! From LMDZ4/libf/phylmd/vdif_kcay.F, version 1.1, 2004/06/22 11:45:36
11 guez 3
12 guez 62 USE dimphy, ONLY: klev, klon
13 guez 108 use yamada_m, only: yamada
14 guez 3
15 guez 227 INTEGER knon
16     ! knon : nombre de points de grille
17    
18     REAL, intent(in):: dt
19 guez 62 ! dt : pas de temps
20 guez 227 real, intent(in):: g
21 guez 62 ! g : g
22 guez 227 REAL zlev(klon, klev+1)
23 guez 62 ! zlev : altitude a chaque niveau (interface inferieure de la couche
24     ! de meme indice)
25 guez 227 REAL zlay(klon, klev)
26 guez 62 ! zlay : altitude au centre de chaque couche
27 guez 227 REAL, intent(in):: u(klon, klev)
28     REAL, intent(in):: v(klon, klev)
29 guez 62 ! u, v : vitesse au centre de chaque couche
30     ! (en entree : la valeur au debut du pas de temps)
31 guez 227 REAL teta(klon, klev)
32 guez 62 ! teta : temperature potentielle au centre de chaque couche
33     ! (en entree : la valeur au debut du pas de temps)
34 guez 227 REAL, intent(in):: cd (:) ! (knon) cdrag, valeur au debut du pas de temps
35     REAL q2(klon, klev+1)
36 guez 62 ! 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 guez 227 REAL q2diag(klon, klev+1)
40     REAL km(klon, klev+1)
41 guez 62 ! 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 guez 227 REAL kn(klon, klev+1)
45 guez 62 ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
46     ! (en sortie : la valeur a la fin du pas de temps)
47 guez 227 real, intent(in):: ustar(:) ! (knon)
48     integer l_mix
49 guez 3
50 guez 227 ! Local:
51    
52     real snstable
53 guez 105 real sq(klon), sqz(klon), zq, long0(klon)
54 guez 3
55 guez 227 INTEGER nlay, nlev
56 guez 62 ! nlay : nombre de couches
57     ! nlev : nombre de niveaux
58 guez 227 REAL unsdz(klon, klev)
59 guez 62 ! unsdz : 1 sur l'epaisseur de couche
60 guez 227 REAL unsdzdec(klon, klev+1)
61 guez 62 ! unsdzdec : 1 sur la distance entre le centre de la couche et le
62     ! centre de la couche inferieure
63 guez 227 REAL q(klon, klev+1)
64 guez 62 ! q : echelle de vitesse au bas de chaque couche
65     ! (valeur a la fin du pas de temps)
66 guez 3
67 guez 227 REAL kmpre(klon, klev+1)
68 guez 62 ! kmpre : km au debut du pas de temps
69 guez 227 REAL qcstat
70 guez 62 ! qcstat : q : solution stationnaire du probleme couple
71     ! (valeur a la fin du pas de temps)
72 guez 227 REAL q2cstat
73 guez 62 ! q2cstat : q2 : solution stationnaire du probleme couple
74     ! (valeur a la fin du pas de temps)
75 guez 3
76 guez 227 REAL long(klon, klev+1)
77 guez 62 ! long : longueur de melange calculee selon Blackadar
78 guez 3
79 guez 62 ! 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 guez 3
93 guez 62 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 guez 3
103 guez 62 ! 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 guez 3
114 guez 62 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 guez 3
124 guez 62 ! 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 guez 3
131 guez 62 REAL kappa
132     REAL long00
133     REAL a1, a2, b1, b2, c1
134     REAL cn1, cn2
135     REAL cm1, cm2, cm3, cm4
136 guez 3
137 guez 62 ! 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 guez 3
142 guez 62 REAL termq
143     REAL termq3
144     REAL termqm2
145     REAL termq3m2
146 guez 3
147 guez 62 ! 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 guez 227 call yamada(knon, g, zlev, zlay, u, v, teta, q2diag, km, kn)
192 guez 62 if (first.and.1.eq.1) then
193     first=.false.
194     q2=q2diag
195     endif
196    
197     DO ilev=1, nlev
198 guez 227 DO igrid=1, knon
199 guez 62 q2(igrid, ilev)=amax1(q2(igrid, ilev), q2min)
200     q(igrid, ilev)=sqrt(q2(igrid, ilev))
201     ENDDO
202     ENDDO
203    
204 guez 227 DO igrid=1, knon
205 guez 62 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 guez 227 DO igrid=1, knon
216 guez 62 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 guez 227 DO igrid=1, knon
223 guez 62 unsdz(igrid, ilay)=1.E+0/(zlev(igrid, ilay+1)-zlev(igrid, ilay))
224     ENDDO
225     ENDDO
226 guez 227 DO igrid=1, knon
227 guez 62 unsdzdec(igrid, 1)=1.E+0/(zlay(igrid, 1)-zlev(igrid, 1))
228     ENDDO
229     DO ilay=2, nlay
230 guez 227 DO igrid=1, knon
231 guez 62 unsdzdec(igrid, ilay)=1.E+0/(zlay(igrid, ilay)-zlay(igrid, ilay-1))
232     ENDDO
233     ENDDO
234 guez 227 DO igrid=1, knon
235 guez 62 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 guez 227 DO igrid=1, knon
241 guez 62 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 guez 227 DO igrid=1, knon
251 guez 62
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 guez 227 DO igrid=1, knon
278 guez 62 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 guez 227 DO igrid=1, knon
287 guez 62 sqz(igrid)=1.e-10
288     sq(igrid)=1.e-10
289     ENDDO
290     do ilev=2, nlev-1
291 guez 227 DO igrid=1, knon
292 guez 62 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 guez 227 DO igrid=1, knon
300 guez 62 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 guez 227 DO igrid=1, knon
308 guez 62 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 guez 227 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 guez 62
379 guez 227 if (sn(igrid, ilev).lt.snstable) then
380     sn(igrid, ilev)=snstable
381     snq2(igrid, ilev)=0.
382     endif
383 guez 62
384 guez 227 if (sm(igrid, ilev).lt.snstable) then
385     sm(igrid, ilev)=snstable
386     smq2(igrid, ilev)=0.
387 guez 62 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 guez 227 DO igrid=1, knon
397 guez 62 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 guez 227 DO igrid=1, knon
404 guez 62 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 guez 227 DO igrid=1, knon
413 guez 62 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 guez 227 DO igrid=1, knon
422 guez 62 ! 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 guez 227 DO igrid=1, knon
521 guez 62 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 guez 227 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 guez 62
535 guez 227 ! calcul final de kn et km
536 guez 62
537 guez 227 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 guez 62
548 guez 227 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 guez 62
554 guez 227 if (sn(igrid, ilev).lt.snstable) then
555     sn(igrid, ilev)=snstable
556     snq2(igrid, ilev)=0.
557     endif
558 guez 62
559 guez 227 if (sm(igrid, ilev).lt.snstable) then
560     sm(igrid, ilev)=snstable
561     smq2(igrid, ilev)=0.
562     endif
563 guez 62
564 guez 227 ! 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 guez 62 enddo
569 guez 227 enddo
570 guez 62
571     END SUBROUTINE vdif_kcay
572    
573     end module vdif_kcay_m

  ViewVC Help
Powered by ViewVC 1.1.21