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

Annotation of /trunk/phylmd/vdif_kcay.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 12 - (hide annotations)
Mon Jul 21 16:05:07 2008 UTC (15 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/vdif_kcay.f
File size: 25121 byte(s) 
-- Minor modification of input/output:

Created procedure "read_logic". Variables of module "logic" are read
by "read_logic" instead of "conf_gcm". Variable "offline" of module
"conf_gcm" is read from namelist instead of "*.def".

Deleted arguments "dtime", "co2_ppm_etat0", "solaire_etat0",
"tabcntr0" and local variables "radpas", "tab_cntrl" of
"phyetat0". "phyetat0" does not read "controle" in "startphy.nc" any
longer. "phyetat0" now reads global attribute "itau_phy" from
"startphy.nc". "phyredem" does not create variable "controle" in
"startphy.nc" any longer. "phyredem" now writes global attribute
"itau_phy" of "startphy.nc". Deleted argument "tabcntr0" of
"printflag". Removed diagnostic messages written by "printflag" for
comparison of the variable "controle" of "startphy.nc" and the
variables read from "*.def" or namelist input.

-- Removing unwanted functionality:

Removed variable "lunout" from module "iniprint", replaced everywhere
by standard output.

Removed case "ocean == 'couple'" in "clmain", "interfsurf_hq" and
"physiq". Removed procedure "interfoce_cpl".

-- Should not change anything at run time:

Automated creation of graphs in documentation. More documentation on
input files.

Converted Fortran files to free format: "phyredem.f90", "printflag.f90".

Split module "clesphy" into "clesphys" and "clesphys2".

Removed variables "conser", "leapf", "forward", "apphys", "apdiss" and
"statcl" from module "logic". Added arguments "conser" to "advect",
"leapf" to "integrd". Added local variables "forward", "leapf",
"apphys", "conser", "apdiss" in "leapfrog".

Added intent attributes.

Deleted arguments "dtime" of "phyredem", "pdtime" of "flxdtdq", "sh"
of "phytrac", "dt" of "yamada".

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