source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/vdif_kcay.F90 @ 5085

Last change on this file since 5085 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.
  • 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.