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

Contents of /trunk/phylmd/vdif_kcay.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 17 - (show 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 !
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, intent(in):: g
36 real rconst
37 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 call yamada(ngrid,g,rconst,plev,temp
232 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