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

Contents of /trunk/libf/phylmd/vdif_kcay.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 12 - (show annotations)
Mon Jul 21 16:05:07 2008 UTC (15 years, 10 months ago) by guez
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 !
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 REAL, intent(in):: dt
35 real g,rconst
36 real plev(klon,klev+1),temp(klon,klev)
37 real ustar(klon),snstable
38 REAL zlev(klon,klev+1)
39 REAL zlay(klon,klev)
40 REAL u(klon,klev)
41 REAL v(klon,klev)
42 REAL teta(klon,klev)
43 REAL 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 call yamada(ngrid,g,rconst,plev,temp
231 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