source: LMDZ6/branches/contrails/libf/phylmd/vdif_kcay.f90 @ 5456

Last change on this file since 5456 was 5268, checked in by abarral, 2 months ago

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