source: trunk/LMDZ.TITAN/libf/phytitan/yamada.F @ 1243

Last change on this file since 1243 was 102, checked in by slebonnois, 14 years ago

SL : corrections et modifications dans phytitan correspondant a celles
faites apres compilation Venus. Titan pas encore compile.
File size: 5.4 KB
Line 
1!
2! $Header: /home/cvsroot/LMDZ4/libf/phylmd/yamada.F,v 1.1 2004/06/22 11:45:36 lmdzadmin Exp $
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)
7c.......................................................................
8      use dimphy
9      IMPLICIT NONE
10#include "dimensions.h"
11c.......................................................................
12c
13c dt : pas de temps
14c g  : g
15c zlev : altitude a chaque niveau (interface inferieure de la couche
16c        de meme indice)
17c zlay : altitude au centre de chaque couche
18c u,v : vitesse au centre de chaque couche
19c       (en entree : la valeur au debut du pas de temps)
20c teta : temperature potentielle au centre de chaque couche
21c        (en entree : la valeur au debut du pas de temps)
22c cd : cdrag
23c      (en entree : la valeur au debut du pas de temps)
24c q2 : $q^2$ au bas de chaque couche
25c      (en entree : la valeur au debut du pas de temps)
26c      (en sortie : la valeur a la fin du pas de temps)
27c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
28c      couche)
29c      (en sortie : la valeur a la fin du pas de temps)
30c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
31c      (en sortie : la valeur a la fin du pas de temps)
32c
33c.......................................................................
34      REAL dt,g,rconst
35      real plev(klon,klev+1),temp(klon,klev)
36      real ustar(klon),snstable
37      REAL zlev(klon,klev+1)
38      REAL zlay(klon,klev)
39      REAL u(klon,klev)
40      REAL v(klon,klev)
41      REAL teta(klon,klev)
42      REAL cd(klon)
43      REAL q2(klon,klev+1)
44      REAL km(klon,klev+1)
45      REAL kn(klon,klev+1)
46      integer l_mix,ngrid
47
48      logical first
49      save first
50      data first/.true./
51
52      integer ig,k
53
54      real ri,zrif,zalpha,zsm
55      real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
56
57      real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
58      real l(klon,klev+1),l0(klon)
59
60      real sq(klon),sqz(klon),zz(klon,klev+1)
61      integer iter
62
63      real ric,rifc,b1,kap
64      save ric,rifc,b1,kap
65      data ric,rifc,b1,kap/0.195,0.191,16.6,0.3/
66
67      real frif,falpha,fsm
68
69      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
70      falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
71      fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
72
73      if (0.eq.1.and.first) then
74      do ig=1,1000
75         ri=(ig-800.)/500.
76         if (ri.lt.ric) then
77            zrif=frif(ri)
78         else
79            zrif=rifc
80         endif
81         if(zrif.lt.0.16) then
82            zalpha=falpha(zrif)
83            zsm=fsm(zrif)
84         else
85            zalpha=1.12
86            zsm=0.085
87         endif
88         print*,ri,rif,zalpha,zsm
89      enddo
90      first=.false.
91      endif
92
93c  Correction d'un bug sauvage a verifier.
94c      do k=2,klevp1
95      do k=2,klev
96                                                          do ig=1,ngrid
97         dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
98         m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)
99     s             /(dz(ig,k)*dz(ig,k))
100         n2(ig,k)=g*2.*(teta(ig,k)-teta(ig,k-1))
101     s            /(teta(ig,k-1)+teta(ig,k))  /dz(ig,k)
102         ri=n2(ig,k)/max(m2(ig,k),1.e-10)
103         if (ri.lt.ric) then
104            rif(ig,k)=frif(ri)
105         else
106            rif(ig,k)=rifc
107         endif
108         if(rif(ig,k).lt.0.16) then
109            alpha(ig,k)=falpha(rif(ig,k))
110            sm(ig,k)=fsm(rif(ig,k))
111         else
112            alpha(ig,k)=1.12
113            sm(ig,k)=0.085
114         endif
115         zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
116                                                          enddo
117      enddo
118
119c iterration pour determiner la longueur de melange
120
121                                                          do ig=1,ngrid
122      l0(ig)=100.
123                                                          enddo
124      do k=2,klev-1
125                                                          do ig=1,ngrid
126        l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
127                                                          enddo
128      enddo
129
130      do iter=1,10
131                                                          do ig=1,ngrid
132         sq(ig)=1.e-10
133         sqz(ig)=1.e-10
134                                                          enddo
135         do k=2,klev-1
136                                                          do ig=1,ngrid
137           q2(ig,k)=l(ig,k)**2*zz(ig,k)
138           l(ig,k)=min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
139     s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10)))
140           zq=sqrt(q2(ig,k))
141           sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
142           sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
143                                                          enddo
144         enddo
145                                                          do ig=1,ngrid
146         l0(ig)=0.2*sqz(ig)/sq(ig)
147                                                          enddo
148c(abd 3 5 2)         print*,'ITER=',iter,'  L0=',l0
149
150      enddo
151
152      do k=2,klev
153                                                          do ig=1,ngrid
154         l(ig,k)=min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
155     s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10)))
156         q2(ig,k)=l(ig,k)**2*zz(ig,k)
157         km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
158         kn(ig,k)=km(ig,k)*alpha(ig,k)
159                                                          enddo
160      enddo
161
162      return
163      end
Note: See TracBrowser for help on using the repository browser.