source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/vdif_kcay.F90 @ 3820

Last change on this file since 3820 was 3818, checked in by millour, 10 years ago

Some partial cleanup on uses of "dimensions.h" in physics.
At this point 3D gcm compiles and bench seems to run fine :-)
EM

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.