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

Annotation of /trunk/Sources/phylmd/yamada.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 76 - (hide annotations)
Fri Nov 15 18:45:49 2013 UTC (10 years, 6 months ago) by guez
Original Path: trunk/phylmd/yamada.f
File size: 5499 byte(s)
Moved everything out of libf.
1 guez 3 !
2     ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/yamada.F,v 1.1 2004/06/22 11:45:36 lmdzadmin Exp $
3     !
4 guez 12 SUBROUTINE yamada(ngrid,g,rconst,plev,temp
5 guez 62 s ,zlev,zlay,u,v,teta,q2,km,kn,ustar
6 guez 3 s ,l_mix)
7     use dimens_m
8     use dimphy
9     IMPLICIT NONE
10     c.......................................................................
11     c.......................................................................
12     c
13     c g : g
14     c zlev : altitude a chaque niveau (interface inferieure de la couche
15     c de meme indice)
16     c zlay : altitude au centre de chaque couche
17     c u,v : vitesse au centre de chaque couche
18     c (en entree : la valeur au debut du pas de temps)
19     c teta : temperature potentielle au centre de chaque couche
20     c (en entree : la valeur au debut du pas de temps)
21     c q2 : $q^2$ au bas de chaque couche
22     c (en entree : la valeur au debut du pas de temps)
23     c (en sortie : la valeur a la fin du pas de temps)
24     c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
25     c couche)
26     c (en sortie : la valeur a la fin du pas de temps)
27     c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
28     c (en sortie : la valeur a la fin du pas de temps)
29     c
30     c.......................................................................
31 guez 17 REAL, intent(in):: g
32     real rconst
33 guez 3 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.eq.1.and.first) then
77     do ig=1,1000
78     ri=(ig-800.)/500.
79     if (ri.lt.ric) then
80     zrif=frif(ri)
81     else
82     zrif=rifc
83     endif
84     if(zrif.lt.0.16) then
85     zalpha=falpha(zrif)
86     zsm=fsm(zrif)
87     else
88     zalpha=1.12
89     zsm=0.085
90     endif
91     print*,ri,rif,zalpha,zsm
92     enddo
93     first=.false.
94     endif
95    
96     c Correction d'un bug sauvage a verifier.
97     c 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,k-1))**2)
102     s /(dz(ig,k)*dz(ig,k))
103     n2(ig,k)=g*2.*(teta(ig,k)-teta(ig,k-1))
104     s /(teta(ig,k-1)+teta(ig,k)) /dz(ig,k)
105     ri=n2(ig,k)/max(m2(ig,k),1.e-10)
106     if (ri.lt.ric) then
107     rif(ig,k)=frif(ri)
108     else
109     rif(ig,k)=rifc
110     endif
111     if(rif(ig,k).lt.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     endif
118     zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
119     enddo
120     enddo
121    
122     c iterration pour determiner la longueur de melange
123    
124     do ig=1,ngrid
125     l0(ig)=100.
126     enddo
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     enddo
131     enddo
132    
133     do iter=1,10
134     do ig=1,ngrid
135     sq(ig)=1.e-10
136     sqz(ig)=1.e-10
137     enddo
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,k)+l0(ig))
142     s ,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     enddo
147     enddo
148     do ig=1,ngrid
149     l0(ig)=0.2*sqz(ig)/sq(ig)
150     enddo
151     c(abd 3 5 2) print*,'ITER=',iter,' L0=',l0
152    
153     enddo
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,k)+l0(ig))
158     s ,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     enddo
163     enddo
164    
165     return
166     end

  ViewVC Help
Powered by ViewVC 1.1.21