source: LMDZ5/trunk/libf/phylmd/vdif_kcay.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

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