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

Annotation of /trunk/phylmd/vdif_kcay.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 118 - (hide annotations)
Thu Dec 18 17:30:24 2014 UTC (9 years, 5 months ago) by guez
File size: 18963 byte(s) 
In file grilles_gcm.nc, renamed variable phis to orog, deleted
variable presnivs.

Removed variable bug_ozone from module clesphys.

In procedure ozonecm, moved computation of sint and cost out of the
loops on horizontal position and vertical level. Inverted the order of
the two loops. We can then move all computations from slat to aprim
out of the loop on vertical levels. Created variable slat2, following
LMDZ. Moved the limitation of column-density of ozone in cell at 1e-12
from radlwsw to ozonecm, following LMDZ.

Removed unused arguments u, albsol, rh, cldfra, rneb, diafra, cldliq,
pmflxr, pmflxs, prfl, psfl of phytrac.

In procedure yamada4, for all the arrays, replaced the dimension klon
by ngrid. At the end of the procedure, for the computation of kmn,kn,
kq and q2, changed the upper limit of the loop index from klon to ngrid.

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