/[lmdze]/trunk/phylmd/yamada.f
ViewVC logotype

Diff of /trunk/phylmd/yamada.f

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

trunk/phylmd/yamada.f revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/phylmd/yamada.f90 revision 81 by guez, Wed Mar 5 14:38:41 2014 UTC
# Line 1  Line 1 
 !  
 ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/yamada.F,v 1.1 2004/06/22 11:45:36 lmdzadmin Exp $  
 !  
       SUBROUTINE yamada(ngrid,g,rconst,plev,temp  
      s   ,zlev,zlay,u,v,teta,q2,km,kn,ustar  
      s   ,l_mix)  
       use dimens_m  
       use dimphy  
       IMPLICIT NONE  
 c.......................................................................  
 c.......................................................................  
 c  
 c g  : g  
 c zlev : altitude a chaque niveau (interface inferieure de la couche  
 c        de meme indice)  
 c zlay : altitude au centre de chaque couche  
 c u,v : vitesse au centre de chaque couche  
 c       (en entree : la valeur au debut du pas de temps)  
 c teta : temperature potentielle au centre de chaque couche  
 c        (en entree : la valeur au debut du pas de temps)  
 c q2 : $q^2$ au bas de chaque couche  
 c      (en entree : la valeur au debut du pas de temps)  
 c      (en sortie : la valeur a la fin du pas de temps)  
 c km : diffusivite turbulente de quantite de mouvement (au bas de chaque  
 c      couche)  
 c      (en sortie : la valeur a la fin du pas de temps)  
 c kn : diffusivite turbulente des scalaires (au bas de chaque couche)  
 c      (en sortie : la valeur a la fin du pas de temps)  
 c  
 c.......................................................................  
       REAL, intent(in):: g  
       real rconst  
       real plev(klon,klev+1),temp(klon,klev)  
       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 q2(klon,klev+1)  
       REAL km(klon,klev+1)  
       REAL kn(klon,klev+1)  
       integer l_mix,ngrid  
   
   
       integer nlay,nlev  
       PARAMETER (nlay=klev)  
       PARAMETER (nlev=klev+1)  
   
       logical first  
       save first  
       data first/.true./  
   
   
       integer ig,k  
   
       real ri,zrif,zalpha,zsm  
       real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)  
   
       real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)  
       real l(klon,klev+1),l0(klon)  
   
       real sq(klon),sqz(klon),zz(klon,klev+1)  
       integer iter  
   
       real ric,rifc,b1,kap  
       save ric,rifc,b1,kap  
       data ric,rifc,b1,kap/0.195,0.191,16.6,0.3/  
   
       real frif,falpha,fsm  
   
       frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))  
       falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)  
       fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))  
   
       if (0.eq.1.and.first) then  
       do ig=1,1000  
          ri=(ig-800.)/500.  
          if (ri.lt.ric) then  
             zrif=frif(ri)  
          else  
             zrif=rifc  
          endif  
          if(zrif.lt.0.16) then  
             zalpha=falpha(zrif)  
             zsm=fsm(zrif)  
          else  
             zalpha=1.12  
             zsm=0.085  
          endif  
          print*,ri,rif,zalpha,zsm  
       enddo  
       first=.false.  
       endif  
   
 c  Correction d'un bug sauvage a verifier.  
 c      do k=2,nlev  
       do k=2,nlay  
                                                           do ig=1,ngrid  
          dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)  
          m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)  
      s             /(dz(ig,k)*dz(ig,k))  
          n2(ig,k)=g*2.*(teta(ig,k)-teta(ig,k-1))  
      s            /(teta(ig,k-1)+teta(ig,k))  /dz(ig,k)  
          ri=n2(ig,k)/max(m2(ig,k),1.e-10)  
          if (ri.lt.ric) then  
             rif(ig,k)=frif(ri)  
          else  
             rif(ig,k)=rifc  
          endif  
          if(rif(ig,k).lt.0.16) then  
             alpha(ig,k)=falpha(rif(ig,k))  
             sm(ig,k)=fsm(rif(ig,k))  
          else  
             alpha(ig,k)=1.12  
             sm(ig,k)=0.085  
          endif  
          zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)  
                                                           enddo  
       enddo  
   
 c iterration pour determiner la longueur de melange  
   
                                                           do ig=1,ngrid  
       l0(ig)=100.  
                                                           enddo  
       do k=2,klev-1  
                                                           do ig=1,ngrid  
         l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))  
                                                           enddo  
       enddo  
   
       do iter=1,10  
                                                           do ig=1,ngrid  
          sq(ig)=1.e-10  
          sqz(ig)=1.e-10  
                                                           enddo  
          do k=2,klev-1  
                                                           do ig=1,ngrid  
            q2(ig,k)=l(ig,k)**2*zz(ig,k)  
            l(ig,k)=min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))  
      s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10)))  
            zq=sqrt(q2(ig,k))  
            sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))  
            sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))  
                                                           enddo  
          enddo  
                                                           do ig=1,ngrid  
          l0(ig)=0.2*sqz(ig)/sq(ig)  
                                                           enddo  
 c(abd 3 5 2)         print*,'ITER=',iter,'  L0=',l0  
   
       enddo  
   
       do k=2,klev  
                                                           do ig=1,ngrid  
          l(ig,k)=min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))  
      s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10)))  
          q2(ig,k)=l(ig,k)**2*zz(ig,k)  
          km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)  
          kn(ig,k)=km(ig,k)*alpha(ig,k)  
                                                           enddo  
       enddo  
1    
2        return  ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/yamada.F,v 1.1 2004/06/22 11:45:36
3        end  ! lmdzadmin Exp $
4    
5    SUBROUTINE yamada(ngrid, g, rconst, plev, temp, zlev, zlay, u, v, teta, q2, &
6        km, kn, ustar, l_mix)
7      USE dimens_m
8      USE dimphy
9      IMPLICIT NONE
10      ! .......................................................................
11      ! .......................................................................
12    
13      ! g  : g
14      ! zlev : altitude a chaque niveau (interface inferieure de la couche
15      ! de meme indice)
16      ! zlay : altitude au centre de chaque couche
17      ! u,v : vitesse au centre de chaque couche
18      ! (en entree : la valeur au debut du pas de temps)
19      ! teta : temperature potentielle au centre de chaque couche
20      ! (en entree : la valeur au debut du pas de temps)
21      ! q2 : $q^2$ au bas de chaque couche
22      ! (en entree : la valeur au debut du pas de temps)
23      ! (en sortie : la valeur a la fin du pas de temps)
24      ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
25      ! couche)
26      ! (en sortie : la valeur a la fin du pas de temps)
27      ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
28      ! (en sortie : la valeur a la fin du pas de temps)
29    
30      ! .......................................................................
31      REAL, INTENT (IN) :: g
32      REAL rconst
33      REAL plev(klon, klev+1), temp(klon, klev)
34      REAL ustar(klon), snstable
35      REAL zlev(klon, klev+1)
36      REAL zlay(klon, klev)
37      REAL u(klon, klev)
38      REAL v(klon, klev)
39      REAL teta(klon, klev)
40      REAL q2(klon, klev+1)
41      REAL km(klon, klev+1)
42      REAL kn(klon, klev+1)
43      INTEGER l_mix, ngrid
44    
45    
46      INTEGER nlay, nlev
47      PARAMETER (nlay=klev)
48      PARAMETER (nlev=klev+1)
49    
50      LOGICAL first
51      SAVE first
52      DATA first/.TRUE./
53    
54    
55      INTEGER ig, k
56    
57      REAL ri, zrif, zalpha, zsm
58      REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
59    
60      REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
61      REAL l(klon, klev+1), l0(klon)
62    
63      REAL sq(klon), sqz(klon), zz(klon, klev+1)
64      INTEGER iter
65    
66      REAL ric, rifc, b1, kap
67      SAVE ric, rifc, b1, kap
68      DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.3/
69    
70      REAL frif, falpha, fsm
71    
72      frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
73      falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
74      fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
75    
76      IF (0==1 .AND. first) THEN
77        DO ig = 1, 1000
78          ri = (ig-800.)/500.
79          IF (ri<ric) THEN
80            zrif = frif(ri)
81          ELSE
82            zrif = rifc
83          END IF
84          IF (zrif<0.16) THEN
85            zalpha = falpha(zrif)
86            zsm = fsm(zrif)
87          ELSE
88            zalpha = 1.12
89            zsm = 0.085
90          END IF
91          PRINT *, ri, rif, zalpha, zsm
92        END DO
93        first = .FALSE.
94      END IF
95    
96      ! Correction d'un bug sauvage a verifier.
97      ! do k=2,nlev
98      DO k = 2, nlay
99        DO ig = 1, ngrid
100          dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
101          m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &
102            k-1))**2)/(dz(ig,k)*dz(ig,k))
103          n2(ig, k) = g*2.*(teta(ig,k)-teta(ig,k-1))/(teta(ig,k-1)+teta(ig,k))/ &
104            dz(ig, k)
105          ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
106          IF (ri<ric) THEN
107            rif(ig, k) = frif(ri)
108          ELSE
109            rif(ig, k) = rifc
110          END IF
111          IF (rif(ig,k)<0.16) THEN
112            alpha(ig, k) = falpha(rif(ig,k))
113            sm(ig, k) = fsm(rif(ig,k))
114          ELSE
115            alpha(ig, k) = 1.12
116            sm(ig, k) = 0.085
117          END IF
118          zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
119        END DO
120      END DO
121    
122      ! iterration pour determiner la longueur de melange
123    
124      DO ig = 1, ngrid
125        l0(ig) = 100.
126      END DO
127      DO k = 2, klev - 1
128        DO ig = 1, ngrid
129          l(ig, k) = l0(ig)*kap*zlev(ig, k)/(kap*zlev(ig,k)+l0(ig))
130        END DO
131      END DO
132    
133      DO iter = 1, 10
134        DO ig = 1, ngrid
135          sq(ig) = 1.E-10
136          sqz(ig) = 1.E-10
137        END DO
138        DO k = 2, klev - 1
139          DO ig = 1, ngrid
140            q2(ig, k) = l(ig, k)**2*zz(ig, k)
141            l(ig, k) = min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig, &
142              k)+l0(ig)), 0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.E-10)))
143            zq = sqrt(q2(ig,k))
144            sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
145            sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
146          END DO
147        END DO
148        DO ig = 1, ngrid
149          l0(ig) = 0.2*sqz(ig)/sq(ig)
150        END DO
151        ! (abd 3 5 2)         print*,'ITER=',iter,'  L0=',l0
152    
153      END DO
154    
155      DO k = 2, klev
156        DO ig = 1, ngrid
157          l(ig, k) = min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig, &
158            k)+l0(ig)), 0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.E-10)))
159          q2(ig, k) = l(ig, k)**2*zz(ig, k)
160          km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
161          kn(ig, k) = km(ig, k)*alpha(ig, k)
162        END DO
163      END DO
164    
165      RETURN
166    END SUBROUTINE yamada
Legend:
Removed from v.76  
changed lines
  Added in v.81
  ViewVC Help
Powered by ViewVC 1.1.21