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

Annotation of /trunk/phylmd/vdif_kcay.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 108 - (hide annotations)
Tue Sep 16 14:00:41 2014 UTC (9 years, 8 months ago) by guez
File size: 18952 byte(s) 
Imported writefield from LMDZ. Close at the end of gcm the files which
were created by writefiled (not done in LMDZ).

Removed procedures for the output of Grads files. Removed calls to
dump2d. In guide, replaced calls to wrgrads by calls to writefield.

In vlspltqs, removed redundant programming of saturation
pressure. Call foeew from module FCTTRE instead.

Bug fix in interpre: size of w exceeding size of correponding actual
argument wg in advtrac.

In leapfrog, call guide until the end of the run, instead of six hours
before the end.

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