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