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

Annotation of /trunk/phylmd/vdif_kcay.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 62 - (hide annotations)
Thu Jul 26 14:37:37 2012 UTC (11 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/vdif_kcay.f90
File size: 18963 byte(s) 
Changed handling of compiler in compilation system.

Removed the prefix letters "y", "p", "t" or "z" in some names of variables.

Replaced calls to NetCDF by calls to NetCDF95.

Extracted "ioget_calendar" procedures from "calendar.f90" into a
separate file.

Extracted to a separate file, "mathop2.f90", procedures that were not
part of the generic interface "mathop" in "mathop.f90".

Removed computation of "dq" in "bilan_dyn", which was not used.

In "iniadvtrac", removed schemes 20 Slopes and 30 Prather. Was not
compatible with declarations of array sizes.

In "clcdrag", "ustarhb", "vdif_kcay", "yamada4" and "coefkz", changed
the size of some arrays from "klon" to "knon".

Removed possible call to "conema3" in "physiq".

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