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

Diff of /trunk/phylmd/yamada.f

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

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