source: LMDZ5/trunk/libf/phylmd/vdif_kcay.F90 @ 2920

Last change on this file since 2920 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

EM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.9 KB
Line 
1
2! $Header$
3
4SUBROUTINE 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
674END SUBROUTINE vdif_kcay
Note: See TracBrowser for help on using the repository browser.