/[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 178 - (hide annotations)
Fri Mar 11 18:47:26 2016 UTC (8 years, 2 months ago) by guez
File size: 18706 byte(s) 
Moved variables date0, deltat, datasz_max, ncvar_ids, point, buff_pos,
buffer, regular from module histcom_var to modules where they are
defined.

Removed procedure ioipslmpp, useless for a sequential program.

Added argument datasz_max to histwrite_real (to avoid circular
dependency with histwrite).

Removed useless variables and computations everywhere.

Changed real litteral constants from default kind to double precision
in lwb, lwu, lwvn, sw1s, swtt, swtt1, swu.

Removed unused arguments: paer of sw, sw1s, sw2s, swclr; pcldsw of
sw1s, sw2s; pdsig, prayl of swr; co2_ppm of clmain, clqh; tsol of
transp_lay; nsrf of screenp; kcrit and kknu of gwstress; pstd of
orosetup.

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