My Project
 All Classes Files Functions Variables Macros
yamada.F
Go to the documentation of this file.
1 !
2 ! $Header$
3 !
4  SUBROUTINE yamada(ngrid,dt,g,rconst,plev,temp
5  s ,zlev,zlay,u,v,teta,cd,q2,km,kn,ustar
6  s ,l_mix)
7  use dimphy
8  IMPLICIT NONE
9 c.......................................................................
10 cym#include "dimensions.h"
11 cym#include "dimphy.h"
12 c.......................................................................
13 c
14 c dt : pas de temps
15 c g : g
16 c zlev : altitude a chaque niveau (interface inferieure de la couche
17 c de meme indice)
18 c zlay : altitude au centre de chaque couche
19 c u,v : vitesse au centre de chaque couche
20 c (en entree : la valeur au debut du pas de temps)
21 c teta : temperature potentielle au centre de chaque couche
22 c (en entree : la valeur au debut du pas de temps)
23 c cd : cdrag
24 c (en entree : la valeur au debut du pas de temps)
25 c q2 : $q^2$ au bas de chaque couche
26 c (en entree : la valeur au debut du pas de temps)
27 c (en sortie : la valeur a la fin du pas de temps)
28 c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
29 c couche)
30 c (en sortie : la valeur a la fin du pas de temps)
31 c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
32 c (en sortie : la valeur a la fin du pas de temps)
33 c
34 c.......................................................................
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)
45  REAL km(klon,klev+1)
46  REAL kn(klon,klev+1)
47  integer l_mix,ngrid
48 
49 
50  integer nlay,nlev
51 cym PARAMETER (nlay=klev)
52 cym PARAMETER (nlev=klev+1)
53 
54  logical first
55  save first
56  data first/.true./
57 c$OMP THREADPRIVATE(first)
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 c$OMP THREADPRIVATE(ric,rifc,b1,kap)
74 
75  real frif,falpha,fsm
76 
77  frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
78  falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
79  fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
80 
81  nlay=klev
82  nlev=klev+1
83 
84  if (0.eq.1.and.first) then
85  do ig=1,1000
86  ri=(ig-800.)/500.
87  if (ri.lt.ric) then
88  zrif=frif(ri)
89  else
90  zrif=rifc
91  endif
92  if(zrif.lt.0.16) then
93  zalpha=falpha(zrif)
94  zsm=fsm(zrif)
95  else
96  zalpha=1.12
97  zsm=0.085
98  endif
99  print*,ri,rif,zalpha,zsm
100  enddo
101  first=.false.
102  endif
103 
104 c Correction d'un bug sauvage a verifier.
105 c do k=2,nlev
106  do k=2,nlay
107  do ig=1,ngrid
108  dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
109  m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)
110  s /(dz(ig,k)*dz(ig,k))
111  n2(ig,k)=g*2.*(teta(ig,k)-teta(ig,k-1))
112  s /(teta(ig,k-1)+teta(ig,k)) /dz(ig,k)
113  ri=n2(ig,k)/max(m2(ig,k),1.e-10)
114  if (ri.lt.ric) then
115  rif(ig,k)=frif(ri)
116  else
117  rif(ig,k)=rifc
118  endif
119  if(rif(ig,k).lt.0.16) then
120  alpha(ig,k)=falpha(rif(ig,k))
121  sm(ig,k)=fsm(rif(ig,k))
122  else
123  alpha(ig,k)=1.12
124  sm(ig,k)=0.085
125  endif
126  zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
127  enddo
128  enddo
129 
130 c iterration pour determiner la longueur de melange
131 
132  do ig=1,ngrid
133  l0(ig)=100.
134  enddo
135  do k=2,klev-1
136  do ig=1,ngrid
137  l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
138  enddo
139  enddo
140 
141  do iter=1,10
142  do ig=1,ngrid
143  sq(ig)=1.e-10
144  sqz(ig)=1.e-10
145  enddo
146  do k=2,klev-1
147  do ig=1,ngrid
148  q2(ig,k)=l(ig,k)**2*zz(ig,k)
149  l(ig,k)=min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
150  s ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10)))
151  zq=sqrt(q2(ig,k))
152  sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
153  sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
154  enddo
155  enddo
156  do ig=1,ngrid
157  l0(ig)=0.2*sqz(ig)/sq(ig)
158  enddo
159 c(abd 3 5 2) print*,'ITER=',iter,' L0=',l0
160 
161  enddo
162 
163  do k=2,klev
164  do ig=1,ngrid
165  l(ig,k)=min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
166  s ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10)))
167  q2(ig,k)=l(ig,k)**2*zz(ig,k)
168  km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
169  kn(ig,k)=km(ig,k)*alpha(ig,k)
170  enddo
171  enddo
172 
173  return
174  end