My Project
 All Classes Files Functions Variables Macros
vdif_kcay.F
Go to the documentation of this file.
1 !
2 ! $Header$
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 dimphy
8  IMPLICIT NONE
9 c.......................................................................
10 cym#include "dimensions.h"
11 cym#include "dimphy.h"
12 c.......................................................................
13 c
14 c dt : pas de temps
15 c g : g
16 c zlev : altitude a chaque niveau (interface inferieure de la couche
17 c de meme indice)
18 c zlay : altitude au centre de chaque couche
19 c u,v : vitesse au centre de chaque couche
20 c (en entree : la valeur au debut du pas de temps)
21 c teta : temperature potentielle au centre de chaque couche
22 c (en entree : la valeur au debut du pas de temps)
23 c cd : cdrag
24 c (en entree : la valeur au debut du pas de temps)
25 c q2 : $q^2$ au bas de chaque couche
26 c (en entree : la valeur au debut du pas de temps)
27 c (en sortie : la valeur a la fin du pas de temps)
28 c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
29 c couche)
30 c (en sortie : la valeur a la fin du pas de temps)
31 c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
32 c (en sortie : la valeur a la fin du pas de temps)
33 c
34 c.......................................................................
35  REAL dt,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 cym PARAMETER (nlay=klev)
199 cym 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$OMP THREADPRIVATE(first)
225 c.......................................................................
226 c traitment des valeur de q2 en entree
227 c.......................................................................
228 c
229 c Initialisation de q2
230  nlay=klev
231  nlev=klev+1
232 
233  call yamada(ngrid,dt,g,rconst,plev,temp
234  s ,zlev,zlay,u,v,teta,cd,q2diag,km,kn,ustar
235  s ,l_mix)
236  if (first.and.1.eq.1) then
237  first=.false.
238  q2=q2diag
239  endif
240 
241  DO ilev=1,nlev
242  DO igrid=1,ngrid
243  q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
244  q(igrid,ilev)=sqrt(q2(igrid,ilev))
245  ENDDO
246  ENDDO
247 c
248  DO igrid=1,ngrid
249  tmp1=cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2)
250  q2(igrid,1)=b1**(2.e+0/3.e+0)*tmp1
251  q2(igrid,1)=amax1(q2(igrid,1),q2min)
252  q(igrid,1)=sqrt(q2(igrid,1))
253  ENDDO
254 c
255 c.......................................................................
256 c les increments verticaux
257 c.......................................................................
258 c
259 c!!!!! allerte !!!!!c
260 c!!!!! zlev n'est pas declare a nlev !!!!!c
261 c!!!!! ---->
262  DO igrid=1,ngrid
263  zlev(igrid,nlev)=zlay(igrid,nlay)
264  & +( zlay(igrid,nlay) - zlev(igrid,nlev-1) )
265  ENDDO
266 c!!!!! <----
267 c!!!!! allerte !!!!!c
268 c
269  DO ilay=1,nlay
270  DO igrid=1,ngrid
271  unsdz(igrid,ilay)=1.e+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
272  ENDDO
273  ENDDO
274  DO igrid=1,ngrid
275  unsdzdec(igrid,1)=1.e+0/(zlay(igrid,1)-zlev(igrid,1))
276  ENDDO
277  DO ilay=2,nlay
278  DO igrid=1,ngrid
279  unsdzdec(igrid,ilay)=1.e+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1))
280  ENDDO
281  ENDDO
282  DO igrid=1,ngrid
283  unsdzdec(igrid,nlay+1)=1.e+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay))
284  ENDDO
285 c
286 c.......................................................................
287 c le cisaillement et le gradient de temperature
288 c.......................................................................
289 c
290  DO igrid=1,ngrid
291  m2(igrid,1)=(unsdzdec(igrid,1)
292  & *u(igrid,1))**2
293  & +(unsdzdec(igrid,1)
294  & *v(igrid,1))**2
295  m(igrid,1)=sqrt(m2(igrid,1))
296  mpre(igrid,1)=m(igrid,1)
297  ENDDO
298 c
299 c-----------------------------------------------------------------------
300  DO ilev=2,nlev-1
301  DO igrid=1,ngrid
302 c-----------------------------------------------------------------------
303 c
304  n2(igrid,ilev)=g*unsdzdec(igrid,ilev)
305  & *(teta(igrid,ilev)-teta(igrid,ilev-1))
306  & /(teta(igrid,ilev)+teta(igrid,ilev-1)) *2.e+0
307 c n2(igrid,ilev)=0.
308 c
309 c --->
310 c on ne sais traiter que les cas stratifies. et l'ajustement
311 c convectif est cense faire en sorte que seul des configurations
312 c stratifiees soient rencontrees en entree de cette routine.
313 c mais, bon ... on sait jamais (meme on sait que n2 prends
314 c quelques valeurs negatives ... parfois) alors :
315 c<---
316 c
317  IF (n2(igrid,ilev).lt.0.e+0) THEN
318  n2(igrid,ilev)=0.e+0
319  ENDIF
320 c
321  m2(igrid,ilev)=(unsdzdec(igrid,ilev)
322  & *(u(igrid,ilev)-u(igrid,ilev-1)))**2
323  & +(unsdzdec(igrid,ilev)
324  & *(v(igrid,ilev)-v(igrid,ilev-1)))**2
325  m(igrid,ilev)=sqrt(m2(igrid,ilev))
326  mpre(igrid,ilev)=m(igrid,ilev)
327 c
328 c-----------------------------------------------------------------------
329  ENDDO
330  ENDDO
331 c-----------------------------------------------------------------------
332 c
333  DO igrid=1,ngrid
334  m2(igrid,nlev)=m2(igrid,nlev-1)
335  m(igrid,nlev)=m(igrid,nlev-1)
336  mpre(igrid,nlev)=m(igrid,nlev)
337  ENDDO
338 c
339 c.......................................................................
340 c calcul des fonctions de stabilite
341 c.......................................................................
342 c
343  if (l_mix.eq.4) then
344  DO igrid=1,ngrid
345  sqz(igrid)=1.e-10
346  sq(igrid)=1.e-10
347  ENDDO
348  do ilev=2,nlev-1
349  DO igrid=1,ngrid
350  zq=sqrt(q2(igrid,ilev))
351  sqz(igrid)
352  . =sqz(igrid)+zq*zlev(igrid,ilev)
353  . *(zlay(igrid,ilev)-zlay(igrid,ilev-1))
354  sq(igrid)=sq(igrid)+zq*(zlay(igrid,ilev)-zlay(igrid,ilev-1))
355  ENDDO
356  enddo
357  DO igrid=1,ngrid
358  long0(igrid)=0.2*sqz(igrid)/sq(igrid)
359  ENDDO
360  else if (l_mix.eq.3) then
361  long0(igrid)=long00
362  endif
363 
364 c (abd 5 2) print*,'LONG0=',long0
365 
366 c-----------------------------------------------------------------------
367  DO ilev=2,nlev-1
368  DO igrid=1,ngrid
369 c-----------------------------------------------------------------------
370 c
371  tmp1=kappa*(zlev(igrid,ilev)-zlev(igrid,1))
372  if (l_mix.ge.10) then
373  long(igrid,ilev)=l_mix
374  else
375  long(igrid,ilev)=tmp1/(1.e+0 + tmp1/long0(igrid))
376  endif
377  long(igrid,ilev)=max(min(long(igrid,ilev)
378  s ,0.5*sqrt(q2(igrid,ilev))/sqrt(max(n2(igrid,ilev),1.e-10)))
379  s ,5.)
380 
381  gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
382  & * n2(igrid,ilev)
383  gm=long(igrid,ilev)**2 / q2(igrid,ilev)
384  & * m2(igrid,ilev)
385 c
386  gninf=.false.
387  gnsup=.false.
388  long(igrid,ilev)=long(igrid,ilev)
389  long(igrid,ilev)=long(igrid,ilev)
390 c
391  IF (gn.lt.gnmin) THEN
392  gninf=.true.
393  gn=gnmin
394  ENDIF
395 c
396  IF (gn.gt.gnmax) THEN
397  gnsup=.true.
398  gn=gnmax
399  ENDIF
400 c
401  sn(igrid,ilev)=cn1/(1.e+0 +cn2*gn)
402  sm(igrid,ilev)=
403  & (cm1+cm2*gn)
404  & /( (1.e+0 +cm3*gn)
405  & *(1.e+0 +cm4*gn) )
406 c
407  IF ((gninf).or.(gnsup)) THEN
408  snq2(igrid,ilev)=0.e+0
409  smq2(igrid,ilev)=0.e+0
410  ELSE
411  snq2(igrid,ilev)=
412  & -gn
413  & *(-cn1*cn2/(1.e+0 +cn2*gn)**2 )
414  smq2(igrid,ilev)=
415  & -gn
416  & *( cm2*(1.e+0 +cm3*gn)
417  & *(1.e+0 +cm4*gn)
418  & -( cm3*(1.e+0 +cm4*gn)
419  & +cm4*(1.e+0 +cm3*gn) )
420  & *(cm1+cm2*gn) )
421  & /( (1.e+0 +cm3*gn)
422  & *(1.e+0 +cm4*gn) )**2
423  ENDIF
424 c
425 c abd
426 c if(ilev.le.57.and.ilev.ge.37) then
427 c print*,'L=',ilev,' GN=',gn,' SM=',sm(igrid,ilev)
428 c endif
429 c --->
430 c la decomposition de Taylor en q2 n'a de sens que
431 c dans les cas stratifies ou sn et sm sont quasi
432 c proportionnels a q2. ailleurs on laisse le meme
433 c algorithme car l'ajustement convectif fait le travail.
434 c mais c'est delirant quand sn et snq2 n'ont pas le meme
435 c signe : dans ces cas, on ne fait pas la decomposition.
436 c<---
437 c
438  IF (snq2(igrid,ilev)*sn(igrid,ilev).le.0.e+0)
439  & snq2(igrid,ilev)=0.e+0
440  IF (smq2(igrid,ilev)*sm(igrid,ilev).le.0.e+0)
441  & smq2(igrid,ilev)=0.e+0
442 c
443 C Correction pour les couches stables.
444 C Schema repris de JHoltzlag Boville, lui meme venant de...
445 
446  if (1.eq.1) then
447  snstable=1.-zlev(igrid,ilev)
448  s /(700.*max(ustar(igrid),0.0001))
449  snstable=1.-zlev(igrid,ilev)/400.
450  snstable=max(snstable,0.)
451  snstable=snstable*snstable
452 
453 c abde print*,'SN ',ilev,sn(1,ilev),snstable
454  if (sn(igrid,ilev).lt.snstable) then
455  sn(igrid,ilev)=snstable
456  snq2(igrid,ilev)=0.
457  endif
458 
459  if (sm(igrid,ilev).lt.snstable) then
460  sm(igrid,ilev)=snstable
461  smq2(igrid,ilev)=0.
462  endif
463 
464  endif
465 
466 c sn : coefficient de stabilite pour n
467 c snq2 : premier terme du developement limite de sn en q2
468 c-----------------------------------------------------------------------
469  ENDDO
470  ENDDO
471 c-----------------------------------------------------------------------
472 c
473 c.......................................................................
474 c calcul de km et kn au debut du pas de temps
475 c.......................................................................
476 c
477  DO igrid=1,ngrid
478  kn(igrid,1)=knmin
479  km(igrid,1)=kmmin
480  kmpre(igrid,1)=km(igrid,1)
481  ENDDO
482 c
483 c-----------------------------------------------------------------------
484  DO ilev=2,nlev-1
485  DO igrid=1,ngrid
486 c-----------------------------------------------------------------------
487 c
488  kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
489  & *sn(igrid,ilev)
490  km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
491  & *sm(igrid,ilev)
492  kmpre(igrid,ilev)=km(igrid,ilev)
493 c
494 c-----------------------------------------------------------------------
495  ENDDO
496  ENDDO
497 c-----------------------------------------------------------------------
498 c
499  DO igrid=1,ngrid
500  kn(igrid,nlev)=kn(igrid,nlev-1)
501  km(igrid,nlev)=km(igrid,nlev-1)
502  kmpre(igrid,nlev)=km(igrid,nlev)
503  ENDDO
504 c
505 c.......................................................................
506 c boucle sur les niveaux 2 a nlev-1
507 c.......................................................................
508 c
509 c---->
510  DO 10001 ilev=2,nlev-1
511 c---->
512  DO 10002 igrid=1,ngrid
513 c
514 c.......................................................................
515 c
516 c calcul des termes sources et puits de l'equation de q2
517 c ------------------------------------------------------
518 c
519  knq3=kn(igrid,ilev)*snq2(igrid,ilev)
520  & /sn(igrid,ilev)
521  kmq3=km(igrid,ilev)*smq2(igrid,ilev)
522  & /sm(igrid,ilev)
523 c
524  termq=0.e+0
525  termq3=0.e+0
526  termqm2=0.e+0
527  termq3m2=0.e+0
528 c
529  tmp1=dt*2.e+0 *km(igrid,ilev)*m2(igrid,ilev)
530  tmp2=dt*2.e+0 *kmq3*m2(igrid,ilev)
531  termqm2=termqm2
532  & +dt*2.e+0 *km(igrid,ilev)*m2(igrid,ilev)
533  & -dt*2.e+0 *kmq3*m2(igrid,ilev)
534  termq3m2=termq3m2
535  & +dt*2.e+0 *kmq3*m2(igrid,ilev)
536 c
537  termq=termq
538  & -dt*2.e+0 *kn(igrid,ilev)*n2(igrid,ilev)
539  & +dt*2.e+0 *knq3*n2(igrid,ilev)
540  termq3=termq3
541  & -dt*2.e+0 *knq3*n2(igrid,ilev)
542 c
543  termq3=termq3
544  & -dt*2.e+0 *q(igrid,ilev)**3 / (b1*long(igrid,ilev))
545 c
546 c.......................................................................
547 c
548 c resolution stationnaire couplee avec le gradient de vitesse local
549 c -----------------------------------------------------------------
550 c
551 c -----{on cherche le cisaillement qui annule l'equation de q^2
552 c supposee en q3}
553 c
554  tmp1=termq+termq3
555  tmp2=termqm2+termq3m2
556  m2cstat=m2(igrid,ilev)
557  & -(tmp1+tmp2)/(dt*2.e+0*km(igrid,ilev))
558  mcstat=sqrt(m2cstat)
559 
560 c abde print*,'M2 L=',ilev,mpre(igrid,ilev),mcstat
561 c
562 c -----{puis on ecrit la valeur de q qui annule l'equation de m
563 c supposee en q3}
564 c
565  IF (ilev.eq.2) THEN
566  kmcstat=1.e+0 / mcstat
567  & *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
568  & *mpre(igrid,ilev+1)
569  & +unsdz(igrid,ilev-1)
570  & *cd(igrid)
571  & *( sqrt(u(igrid,3)**2+v(igrid,3)**2)
572  & -mcstat/unsdzdec(igrid,ilev)
573  & -mpre(igrid,ilev+1)/unsdzdec(igrid,ilev+1) )**2)
574  & /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
575  ELSE
576  kmcstat=1.e+0 / mcstat
577  & *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
578  & *mpre(igrid,ilev+1)
579  & +unsdz(igrid,ilev-1)*kmpre(igrid,ilev-1)
580  & *mpre(igrid,ilev-1) )
581  & /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
582  ENDIF
583  tmp2=kmcstat
584  & /( sm(igrid,ilev)/q2(igrid,ilev) )
585  & /long(igrid,ilev)
586  qcstat=tmp2**(1.e+0/3.e+0)
587  q2cstat=qcstat**2
588 c
589 c.......................................................................
590 c
591 c choix de la solution finale
592 c ---------------------------
593 c
594  q(igrid,ilev)=qcstat
595  q2(igrid,ilev)=q2cstat
596  m(igrid,ilev)=mcstat
597 c abd if(ilev.le.57.and.ilev.ge.37) then
598 c print*,'L=',ilev,' M2=',m2(igrid,ilev),m2cstat,
599 c s 'N2=',n2(igrid,ilev)
600 c abd endif
601  m2(igrid,ilev)=m2cstat
602 c
603 c --->
604 c pour des raisons simples q2 est minore
605 c<---
606 c
607  IF (q2(igrid,ilev).lt.q2min) THEN
608  q2(igrid,ilev)=q2min
609  q(igrid,ilev)=sqrt(q2min)
610  ENDIF
611 c
612 c.......................................................................
613 c
614 c calcul final de kn et km
615 c ------------------------
616 c
617  gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
618  & * n2(igrid,ilev)
619  IF (gn.lt.gnmin) gn=gnmin
620  IF (gn.gt.gnmax) gn=gnmax
621  sn(igrid,ilev)=cn1/(1.e+0 +cn2*gn)
622  sm(igrid,ilev)=
623  & (cm1+cm2*gn)
624  & /( (1.e+0 +cm3*gn)*(1.e+0 +cm4*gn) )
625  kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
626  & *sn(igrid,ilev)
627  km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
628  & *sm(igrid,ilev)
629 c abd
630 c if(ilev.le.57.and.ilev.ge.37) then
631 c print*,'L=',ilev,' GN=',gn,' SM=',sm(igrid,ilev)
632 c endif
633 c
634 c.......................................................................
635 c
636 10002 CONTINUE
637 c
638 10001 CONTINUE
639 c
640 c.......................................................................
641 c
642 c
643  DO igrid=1,ngrid
644  kn(igrid,1)=knmin
645  km(igrid,1)=kmmin
646 c kn(igrid,1)=cd(igrid)
647 c km(igrid,1)=cd(igrid)
648  q2(igrid,nlev)=q2(igrid,nlev-1)
649  q(igrid,nlev)=q(igrid,nlev-1)
650  kn(igrid,nlev)=kn(igrid,nlev-1)
651  km(igrid,nlev)=km(igrid,nlev-1)
652  ENDDO
653 c
654 c CALCUL DE LA DIFFUSION VERTICALE DE Q2
655  if (1.eq.1) then
656 
657  do ilev=2,klev-1
658  sss=sss+plev(1,ilev-1)-plev(1,ilev+1)
659  sssq=sssq+(plev(1,ilev-1)-plev(1,ilev+1))*q2(1,ilev)
660  enddo
661 c print*,'Q2moy avant',sssq/sss
662 c print*,'Q2q20 ',(q2(1,ilev),ilev=1,10)
663 c print*,'Q2km0 ',(km(1,ilev),ilev=1,10)
664 c ! C'est quoi ca qu'etait dans l'original???
665 c do igrid=1,ngrid
666 c q2(igrid,1)=10.
667 c enddo
668 c q2s=q2
669 c do iii=1,10
670 c call vdif_q2(dt,g,rconst,plev,temp,km,q2)
671 c do ilev=1,klev+1
672 c write(iii+49,*) q2(1,ilev),zlev(1,ilev)
673 c enddo
674 c enddo
675 c stop
676 c do ilev=1,klev
677 c print*,zlev(1,ilev),q2s(1,ilev),q2(1,ilev)
678 c enddo
679 c q2s=q2-q2s
680 c do ilev=1,klev
681 c print*,q2s(1,ilev),zlev(1,ilev)
682 c enddo
683  do ilev=2,klev-1
684  sss=sss+plev(1,ilev-1)-plev(1,ilev+1)
685  sssq=sssq+(plev(1,ilev-1)-plev(1,ilev+1))*q2(1,ilev)
686  enddo
687  print*,'Q2moy apres',sssq/sss
688 c
689 c
690  do ilev=1,nlev
691  do igrid=1,ngrid
692  q2(igrid,ilev)=max(q2(igrid,ilev),q2min)
693  q(igrid,ilev)=sqrt(q2(igrid,ilev))
694 
695 c.......................................................................
696 c
697 c calcul final de kn et km
698 c ------------------------
699 c
700  gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
701  & * n2(igrid,ilev)
702  IF (gn.lt.gnmin) gn=gnmin
703  IF (gn.gt.gnmax) gn=gnmax
704  sn(igrid,ilev)=cn1/(1.e+0 +cn2*gn)
705  sm(igrid,ilev)=
706  & (cm1+cm2*gn)
707  & /( (1.e+0 +cm3*gn)*(1.e+0 +cm4*gn) )
708 C Correction pour les couches stables.
709 C Schema repris de JHoltzlag Boville, lui meme venant de...
710 
711  if (1.eq.1) then
712  snstable=1.-zlev(igrid,ilev)
713  s /(700.*max(ustar(igrid),0.0001))
714  snstable=1.-zlev(igrid,ilev)/400.
715  snstable=max(snstable,0.)
716  snstable=snstable*snstable
717 
718 c abde print*,'SN ',ilev,sn(1,ilev),snstable
719  if (sn(igrid,ilev).lt.snstable) then
720  sn(igrid,ilev)=snstable
721  snq2(igrid,ilev)=0.
722  endif
723 
724  if (sm(igrid,ilev).lt.snstable) then
725  sm(igrid,ilev)=snstable
726  smq2(igrid,ilev)=0.
727  endif
728 
729  endif
730 
731 c sn : coefficient de stabilite pour n
732  kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
733  & *sn(igrid,ilev)
734  km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
735 c
736  enddo
737  enddo
738 c print*,'Q2km1 ',(km(1,ilev),ilev=1,10)
739 
740  endif
741 
742  RETURN
743  END