/[lmdze]/trunk/Sources/phylmd/vdif_kcay.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/vdif_kcay.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 226 by guez, Tue Mar 22 16:31:39 2016 UTC revision 227 by guez, Thu Nov 2 15:47:03 2017 UTC
# Line 4  module vdif_kcay_m Line 4  module vdif_kcay_m
4    
5  contains  contains
6    
7    SUBROUTINE vdif_kcay(ngrid, dt, g, plev, zlev, zlay, u, v, teta, cd, q2, &    SUBROUTINE vdif_kcay(knon, dt, g, zlev, zlay, u, v, teta, cd, q2, q2diag, &
8         q2diag, km, kn, ustar, l_mix)         km, kn, ustar, l_mix)
9    
10      ! From LMDZ4/libf/phylmd/vdif_kcay.F, version 1.1 2004/06/22 11:45:36      ! From LMDZ4/libf/phylmd/vdif_kcay.F, version 1.1, 2004/06/22 11:45:36
11    
12      USE dimphy, ONLY: klev, klon      USE dimphy, ONLY: klev, klon
13      use yamada_m, only: yamada      use yamada_m, only: yamada
14    
15      INTEGER ngrid      INTEGER knon
16        ! knon : nombre de points de grille
17    
18        REAL, intent(in):: dt
19      ! dt : pas de temps      ! dt : pas de temps
20        real, intent(in):: g
21      ! g : g      ! g : g
22        REAL zlev(klon, klev+1)
23      ! zlev : altitude a chaque niveau (interface inferieure de la couche      ! zlev : altitude a chaque niveau (interface inferieure de la couche
24      ! de meme indice)      ! de meme indice)
25        REAL zlay(klon, klev)
26      ! zlay : altitude au centre de chaque couche      ! zlay : altitude au centre de chaque couche
27        REAL, intent(in):: u(klon, klev)
28        REAL, intent(in):: v(klon, klev)
29      ! u, v : vitesse au centre de chaque couche      ! u, v : vitesse au centre de chaque couche
30      ! (en entree : la valeur au debut du pas de temps)      ! (en entree : la valeur au debut du pas de temps)
31        REAL teta(klon, klev)
32      ! teta : temperature potentielle au centre de chaque couche      ! teta : temperature potentielle au centre de chaque couche
33      ! (en entree : la valeur au debut du pas de temps)      ! (en entree : la valeur au debut du pas de temps)
34        REAL, intent(in):: cd (:) ! (knon) cdrag, valeur au debut du pas de temps
35        REAL q2(klon, klev+1)
36      ! q2 : $q^2$ au bas de chaque couche      ! q2 : $q^2$ au bas de chaque couche
37      ! (en entree : la valeur au debut du pas de temps)      ! (en entree : la valeur au debut du pas de temps)
38      ! (en sortie : la valeur a la fin du pas de temps)      ! (en sortie : la valeur a la fin du pas de temps)
39        REAL q2diag(klon, klev+1)
40        REAL km(klon, klev+1)
41      ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque      ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
42      ! couche)      ! couche)
43      ! (en sortie : la valeur a la fin du pas de temps)      ! (en sortie : la valeur a la fin du pas de temps)
44        REAL kn(klon, klev+1)
45      ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)      ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
46      ! (en sortie : la valeur a la fin du pas de temps)      ! (en sortie : la valeur a la fin du pas de temps)
47        real, intent(in):: ustar(:) ! (knon)
48        integer l_mix
49    
50      REAL, intent(in):: dt      ! Local:
     real, intent(in):: g  
     real plev(klon, klev+1)  
     real ustar(klon), snstable  
     REAL zlev(klon, klev+1)  
     REAL zlay(klon, klev)  
     REAL u(klon, klev)  
     REAL v(klon, klev)  
     REAL teta(klon, klev)  
     REAL, intent(in):: cd (:) ! (ngrid) cdrag, valeur au debut du pas de temps  
     REAL q2(klon, klev+1)  
     REAL q2diag(klon, klev+1)  
     REAL km(klon, klev+1)  
     REAL kn(klon, klev+1)  
     real sq(klon), sqz(klon), zq, long0(klon)  
51    
52      integer l_mix      real snstable
53        real sq(klon), sqz(klon), zq, long0(klon)
54    
55        INTEGER nlay, nlev
56      ! nlay : nombre de couches      ! nlay : nombre de couches
57      ! nlev : nombre de niveaux      ! nlev : nombre de niveaux
58      ! ngrid : nombre de points de grille      REAL unsdz(klon, klev)
59      ! unsdz : 1 sur l'epaisseur de couche      ! unsdz : 1 sur l'epaisseur de couche
60        REAL unsdzdec(klon, klev+1)
61      ! unsdzdec : 1 sur la distance entre le centre de la couche et le      ! unsdzdec : 1 sur la distance entre le centre de la couche et le
62      ! centre de la couche inferieure      ! centre de la couche inferieure
63        REAL q(klon, klev+1)
64      ! q : echelle de vitesse au bas de chaque couche      ! q : echelle de vitesse au bas de chaque couche
65      ! (valeur a la fin du pas de temps)      ! (valeur a la fin du pas de temps)
66    
67      INTEGER nlay, nlev      REAL kmpre(klon, klev+1)
     REAL unsdz(klon, klev)  
     REAL unsdzdec(klon, klev+1)  
     REAL q(klon, klev+1)  
   
68      ! kmpre : km au debut du pas de temps      ! kmpre : km au debut du pas de temps
69        REAL qcstat
70      ! qcstat : q : solution stationnaire du probleme couple      ! qcstat : q : solution stationnaire du probleme couple
71      ! (valeur a la fin du pas de temps)      ! (valeur a la fin du pas de temps)
72        REAL q2cstat
73      ! q2cstat : q2 : solution stationnaire du probleme couple      ! q2cstat : q2 : solution stationnaire du probleme couple
74      ! (valeur a la fin du pas de temps)      ! (valeur a la fin du pas de temps)
75    
     REAL kmpre(klon, klev+1)  
     REAL qcstat  
     REAL q2cstat  
     real sss, sssq  
   
     ! long : longueur de melange calculee selon Blackadar  
   
76      REAL long(klon, klev+1)      REAL long(klon, klev+1)
77        ! long : longueur de melange calculee selon Blackadar
78    
79      ! kmq3 : terme en q^3 dans le developpement de km      ! kmq3 : terme en q^3 dans le developpement de km
80      ! (valeur au debut du pas de temps)      ! (valeur au debut du pas de temps)
# Line 190  contains Line 188  contains
188    
189      ! Initialisation de q2      ! Initialisation de q2
190    
191      call yamada(ngrid, g, zlev, zlay, u, v, teta, q2diag, km, kn)      call yamada(knon, g, zlev, zlay, u, v, teta, q2diag, km, kn)
192      if (first.and.1.eq.1) then      if (first.and.1.eq.1) then
193         first=.false.         first=.false.
194         q2=q2diag         q2=q2diag
195      endif      endif
196    
197      DO ilev=1, nlev      DO ilev=1, nlev
198         DO igrid=1, ngrid         DO igrid=1, knon
199            q2(igrid, ilev)=amax1(q2(igrid, ilev), q2min)            q2(igrid, ilev)=amax1(q2(igrid, ilev), q2min)
200            q(igrid, ilev)=sqrt(q2(igrid, ilev))            q(igrid, ilev)=sqrt(q2(igrid, ilev))
201         ENDDO         ENDDO
202      ENDDO      ENDDO
203    
204      DO igrid=1, ngrid      DO igrid=1, knon
205         tmp1=cd(igrid)*(u(igrid, 1)**2+v(igrid, 1)**2)         tmp1=cd(igrid)*(u(igrid, 1)**2+v(igrid, 1)**2)
206         q2(igrid, 1)=b1**(2.E+0/3.E+0)*tmp1         q2(igrid, 1)=b1**(2.E+0/3.E+0)*tmp1
207         q2(igrid, 1)=amax1(q2(igrid, 1), q2min)         q2(igrid, 1)=amax1(q2(igrid, 1), q2min)
# Line 214  contains Line 212  contains
212    
213      ! allerte !c      ! allerte !c
214      ! zlev n'est pas declare a nlev !c      ! zlev n'est pas declare a nlev !c
215      DO igrid=1, ngrid      DO igrid=1, knon
216         zlev(igrid, nlev)=zlay(igrid, nlay) &         zlev(igrid, nlev)=zlay(igrid, nlay) &
217              +( zlay(igrid, nlay) - zlev(igrid, nlev-1) )              +( zlay(igrid, nlay) - zlev(igrid, nlev-1) )
218      ENDDO      ENDDO
219      ! allerte !c      ! allerte !c
220    
221      DO ilay=1, nlay      DO ilay=1, nlay
222         DO igrid=1, ngrid         DO igrid=1, knon
223            unsdz(igrid, ilay)=1.E+0/(zlev(igrid, ilay+1)-zlev(igrid, ilay))            unsdz(igrid, ilay)=1.E+0/(zlev(igrid, ilay+1)-zlev(igrid, ilay))
224         ENDDO         ENDDO
225      ENDDO      ENDDO
226      DO igrid=1, ngrid      DO igrid=1, knon
227         unsdzdec(igrid, 1)=1.E+0/(zlay(igrid, 1)-zlev(igrid, 1))         unsdzdec(igrid, 1)=1.E+0/(zlay(igrid, 1)-zlev(igrid, 1))
228      ENDDO      ENDDO
229      DO ilay=2, nlay      DO ilay=2, nlay
230         DO igrid=1, ngrid         DO igrid=1, knon
231            unsdzdec(igrid, ilay)=1.E+0/(zlay(igrid, ilay)-zlay(igrid, ilay-1))            unsdzdec(igrid, ilay)=1.E+0/(zlay(igrid, ilay)-zlay(igrid, ilay-1))
232         ENDDO         ENDDO
233      ENDDO      ENDDO
234      DO igrid=1, ngrid      DO igrid=1, knon
235         unsdzdec(igrid, nlay+1)=1.E+0/(zlev(igrid, nlay+1)-zlay(igrid, nlay))         unsdzdec(igrid, nlay+1)=1.E+0/(zlev(igrid, nlay+1)-zlay(igrid, nlay))
236      ENDDO      ENDDO
237    
238      ! le cisaillement et le gradient de temperature      ! le cisaillement et le gradient de temperature
239    
240      DO igrid=1, ngrid      DO igrid=1, knon
241         m2(igrid, 1)=(unsdzdec(igrid, 1) &         m2(igrid, 1)=(unsdzdec(igrid, 1) &
242              *u(igrid, 1))**2 &              *u(igrid, 1))**2 &
243              +(unsdzdec(igrid, 1) &              +(unsdzdec(igrid, 1) &
# Line 249  contains Line 247  contains
247      ENDDO      ENDDO
248    
249      DO ilev=2, nlev-1      DO ilev=2, nlev-1
250         DO igrid=1, ngrid         DO igrid=1, knon
251    
252            n2(igrid, ilev)=g*unsdzdec(igrid, ilev) &            n2(igrid, ilev)=g*unsdzdec(igrid, ilev) &
253                 *(teta(igrid, ilev)-teta(igrid, ilev-1)) &                 *(teta(igrid, ilev)-teta(igrid, ilev-1)) &
# Line 276  contains Line 274  contains
274         ENDDO         ENDDO
275      ENDDO      ENDDO
276    
277      DO igrid=1, ngrid      DO igrid=1, knon
278         m2(igrid, nlev)=m2(igrid, nlev-1)         m2(igrid, nlev)=m2(igrid, nlev-1)
279         m(igrid, nlev)=m(igrid, nlev-1)         m(igrid, nlev)=m(igrid, nlev-1)
280         mpre(igrid, nlev)=m(igrid, nlev)         mpre(igrid, nlev)=m(igrid, nlev)
# Line 285  contains Line 283  contains
283      ! calcul des fonctions de stabilite      ! calcul des fonctions de stabilite
284    
285      if (l_mix.eq.4) then      if (l_mix.eq.4) then
286         DO igrid=1, ngrid         DO igrid=1, knon
287            sqz(igrid)=1.e-10            sqz(igrid)=1.e-10
288            sq(igrid)=1.e-10            sq(igrid)=1.e-10
289         ENDDO         ENDDO
290         do ilev=2, nlev-1         do ilev=2, nlev-1
291            DO igrid=1, ngrid            DO igrid=1, knon
292               zq=sqrt(q2(igrid, ilev))               zq=sqrt(q2(igrid, ilev))
293               sqz(igrid) &               sqz(igrid) &
294                    =sqz(igrid)+zq*zlev(igrid, ilev) &                    =sqz(igrid)+zq*zlev(igrid, ilev) &
# Line 298  contains Line 296  contains
296               sq(igrid)=sq(igrid)+zq*(zlay(igrid, ilev)-zlay(igrid, ilev-1))               sq(igrid)=sq(igrid)+zq*(zlay(igrid, ilev)-zlay(igrid, ilev-1))
297            ENDDO            ENDDO
298         enddo         enddo
299         DO igrid=1, ngrid         DO igrid=1, knon
300            long0(igrid)=0.2*sqz(igrid)/sq(igrid)            long0(igrid)=0.2*sqz(igrid)/sq(igrid)
301         ENDDO         ENDDO
302      else if (l_mix.eq.3) then      else if (l_mix.eq.3) then
# Line 306  contains Line 304  contains
304      endif      endif
305    
306      DO ilev=2, nlev-1      DO ilev=2, nlev-1
307         DO igrid=1, ngrid         DO igrid=1, knon
308            tmp1=kappa*(zlev(igrid, ilev)-zlev(igrid, 1))            tmp1=kappa*(zlev(igrid, ilev)-zlev(igrid, 1))
309            if (l_mix.ge.10) then            if (l_mix.ge.10) then
310               long(igrid, ilev)=l_mix               long(igrid, ilev)=l_mix
# Line 372  contains Line 370  contains
370            ! Correction pour les couches stables.            ! Correction pour les couches stables.
371            ! Schema repris de JHoltzlag Boville, lui meme venant de...            ! Schema repris de JHoltzlag Boville, lui meme venant de...
372    
373            if (1.eq.1) then            snstable=1.-zlev(igrid, ilev) &
374               snstable=1.-zlev(igrid, ilev) &                 /(700.*max(ustar(igrid), 0.0001))
375                    /(700.*max(ustar(igrid), 0.0001))            snstable=1.-zlev(igrid, ilev)/400.
376               snstable=1.-zlev(igrid, ilev)/400.            snstable=max(snstable, 0.)
377               snstable=max(snstable, 0.)            snstable=snstable*snstable
378               snstable=snstable*snstable  
379              if (sn(igrid, ilev).lt.snstable) then
380               if (sn(igrid, ilev).lt.snstable) then               sn(igrid, ilev)=snstable
381                  sn(igrid, ilev)=snstable               snq2(igrid, ilev)=0.
382                  snq2(igrid, ilev)=0.            endif
383               endif  
384              if (sm(igrid, ilev).lt.snstable) then
385               if (sm(igrid, ilev).lt.snstable) then               sm(igrid, ilev)=snstable
386                  sm(igrid, ilev)=snstable               smq2(igrid, ilev)=0.
                 smq2(igrid, ilev)=0.  
              endif  
387            endif            endif
388    
389            ! sn : coefficient de stabilite pour n            ! sn : coefficient de stabilite pour n
# Line 397  contains Line 393  contains
393    
394      ! calcul de km et kn au debut du pas de temps      ! calcul de km et kn au debut du pas de temps
395    
396      DO igrid=1, ngrid      DO igrid=1, knon
397         kn(igrid, 1)=knmin         kn(igrid, 1)=knmin
398         km(igrid, 1)=kmmin         km(igrid, 1)=kmmin
399         kmpre(igrid, 1)=km(igrid, 1)         kmpre(igrid, 1)=km(igrid, 1)
400      ENDDO      ENDDO
401    
402      DO ilev=2, nlev-1      DO ilev=2, nlev-1
403         DO igrid=1, ngrid         DO igrid=1, knon
404            kn(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev) &            kn(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev) &
405                 *sn(igrid, ilev)                 *sn(igrid, ilev)
406            km(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev) &            km(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev) &
# Line 413  contains Line 409  contains
409         ENDDO         ENDDO
410      ENDDO      ENDDO
411    
412      DO igrid=1, ngrid      DO igrid=1, knon
413         kn(igrid, nlev)=kn(igrid, nlev-1)         kn(igrid, nlev)=kn(igrid, nlev-1)
414         km(igrid, nlev)=km(igrid, nlev-1)         km(igrid, nlev)=km(igrid, nlev-1)
415         kmpre(igrid, nlev)=km(igrid, nlev)         kmpre(igrid, nlev)=km(igrid, nlev)
# Line 422  contains Line 418  contains
418      ! boucle sur les niveaux 2 a nlev-1      ! boucle sur les niveaux 2 a nlev-1
419    
420      DO ilev=2, nlev-1      DO ilev=2, nlev-1
421         DO igrid=1, ngrid         DO igrid=1, knon
422            ! calcul des termes sources et puits de l'equation de q2            ! calcul des termes sources et puits de l'equation de q2
423    
424            knq3=kn(igrid, ilev)*snq2(igrid, ilev) &            knq3=kn(igrid, ilev)*snq2(igrid, ilev) &
# Line 521  contains Line 517  contains
517         end DO         end DO
518      end DO      end DO
519    
520      DO igrid=1, ngrid      DO igrid=1, knon
521         kn(igrid, 1)=knmin         kn(igrid, 1)=knmin
522         km(igrid, 1)=kmmin         km(igrid, 1)=kmmin
523         q2(igrid, nlev)=q2(igrid, nlev-1)         q2(igrid, nlev)=q2(igrid, nlev-1)
# Line 531  contains Line 527  contains
527      ENDDO      ENDDO
528    
529      ! CALCUL DE LA DIFFUSION VERTICALE DE Q2      ! CALCUL DE LA DIFFUSION VERTICALE DE Q2
530      if (1.eq.1) then      do ilev=1, nlev
531         do ilev=2, klev-1         do igrid=1, knon
532            sss=sss+plev(1, ilev-1)-plev(1, ilev+1)            q2(igrid, ilev)=max(q2(igrid, ilev), q2min)
533            sssq=sssq+(plev(1, ilev-1)-plev(1, ilev+1))*q2(1, ilev)            q(igrid, ilev)=sqrt(q2(igrid, ilev))
534         enddo  
535         do ilev=2, klev-1            ! calcul final de kn et km
536            sss=sss+plev(1, ilev-1)-plev(1, ilev+1)  
537            sssq=sssq+(plev(1, ilev-1)-plev(1, ilev+1))*q2(1, ilev)            gn=-long(igrid, ilev)**2 / q2(igrid, ilev) &
538         enddo                 * n2(igrid, ilev)
539         print*, 'Q2moy apres', sssq/sss            IF (gn.lt.gnmin) gn=gnmin
540              IF (gn.gt.gnmax) gn=gnmax
541              sn(igrid, ilev)=cn1/(1.E+0 +cn2*gn)
542              sm(igrid, ilev)= &
543                   (cm1+cm2*gn) &
544                   /( (1.E+0 +cm3*gn)*(1.E+0 +cm4*gn) )
545              ! Correction pour les couches stables.
546              ! Schema repris de JHoltzlag Boville, lui meme venant de...
547    
548         do ilev=1, nlev            snstable=1.-zlev(igrid, ilev) &
549            do igrid=1, ngrid                 /(700.*max(ustar(igrid), 0.0001))
550               q2(igrid, ilev)=max(q2(igrid, ilev), q2min)            snstable=1.-zlev(igrid, ilev)/400.
551               q(igrid, ilev)=sqrt(q2(igrid, ilev))            snstable=max(snstable, 0.)
552              snstable=snstable*snstable
553               ! calcul final de kn et km  
554              if (sn(igrid, ilev).lt.snstable) then
555               gn=-long(igrid, ilev)**2 / q2(igrid, ilev) &               sn(igrid, ilev)=snstable
556                    * n2(igrid, ilev)               snq2(igrid, ilev)=0.
557               IF (gn.lt.gnmin) gn=gnmin            endif
558               IF (gn.gt.gnmax) gn=gnmax  
559               sn(igrid, ilev)=cn1/(1.E+0 +cn2*gn)            if (sm(igrid, ilev).lt.snstable) then
560               sm(igrid, ilev)= &               sm(igrid, ilev)=snstable
561                    (cm1+cm2*gn) &               smq2(igrid, ilev)=0.
562                    /( (1.E+0 +cm3*gn)*(1.E+0 +cm4*gn) )            endif
563               ! Correction pour les couches stables.  
564               ! Schema repris de JHoltzlag Boville, lui meme venant de...            ! sn : coefficient de stabilite pour n
565              kn(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev) &
566               if (1.eq.1) then                 *sn(igrid, ilev)
567                  snstable=1.-zlev(igrid, ilev) &            km(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev)
                      /(700.*max(ustar(igrid), 0.0001))  
                 snstable=1.-zlev(igrid, ilev)/400.  
                 snstable=max(snstable, 0.)  
                 snstable=snstable*snstable  
   
                 if (sn(igrid, ilev).lt.snstable) then  
                    sn(igrid, ilev)=snstable  
                    snq2(igrid, ilev)=0.  
                 endif  
   
                 if (sm(igrid, ilev).lt.snstable) then  
                    sm(igrid, ilev)=snstable  
                    smq2(igrid, ilev)=0.  
                 endif  
              endif  
   
              ! sn : coefficient de stabilite pour n  
              kn(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev) &  
                   *sn(igrid, ilev)  
              km(igrid, ilev)=long(igrid, ilev)*q(igrid, ilev)  
           enddo  
568         enddo         enddo
569      endif      enddo
570    
571    END SUBROUTINE vdif_kcay    END SUBROUTINE vdif_kcay
572    
Legend:
Removed from v.226  
changed lines
  Added in v.227
  ViewVC Help
Powered by ViewVC 1.1.21