source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/vdif_kcay.F90 @ 3839

Last change on this file since 3839 was 1999, checked in by Laurent Fairhead, 10 years ago

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