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

Annotation of /trunk/phylmd/vdif_kcay.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 17 - (hide annotations)
Tue Aug 5 13:31:32 2008 UTC (15 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/vdif_kcay.f
File size: 25146 byte(s) 
Created rule for "compare_sampl_*" files in
"Documentation/Manuel_LMDZE.texfol/Graphiques/GNUmakefile".

Extracted "qcheck", "radiornpb", "minmaxqfi" into separate files.

Read pressure coordinate of ozone coefficients once per run instead of
every day.

Added some "intent" attributes.

Added argument "nq" to "ini_histday". Replaced calls to "gr_fi_ecrit"
by calls to "gr_phy_write_2d". "Sigma_O3_Royer" is written to
"histday.nc" only if "nq >= 4". Moved "ini_histrac" to module
"ini_hist".

Compute "zmasse" in "physiq", pass it to "phytrac".

Removed computations of "pftsol*" and "ppsrf*" in "phytrac".

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