source: trunk/LMDZ.VENUS/libf/phyvenus/vdif_kcay.F @ 1530

Last change on this file since 1530 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

  • Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
  • got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

EM

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