source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/yamada.F @ 1373

Last change on this file since 1373 was 987, checked in by Laurent Fairhead, 16 years ago

Du nettoyage sur le parallelisme, inclusion de nouvelles interfaces pour OPA9
et ORCHIDEE YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.5 KB
Line 
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
9c.......................................................................
10cym#include "dimensions.h"
11cym#include "dimphy.h"
12c.......................................................................
13c
14c dt : pas de temps
15c g  : g
16c zlev : altitude a chaque niveau (interface inferieure de la couche
17c        de meme indice)
18c zlay : altitude au centre de chaque couche
19c u,v : vitesse au centre de chaque couche
20c       (en entree : la valeur au debut du pas de temps)
21c teta : temperature potentielle au centre de chaque couche
22c        (en entree : la valeur au debut du pas de temps)
23c cd : cdrag
24c      (en entree : la valeur au debut du pas de temps)
25c q2 : $q^2$ au bas de chaque couche
26c      (en entree : la valeur au debut du pas de temps)
27c      (en sortie : la valeur a la fin du pas de temps)
28c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
29c      couche)
30c      (en sortie : la valeur a la fin du pas de temps)
31c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
32c      (en sortie : la valeur a la fin du pas de temps)
33c
34c.......................................................................
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
51cym      PARAMETER (nlay=klev)
52cym      PARAMETER (nlev=klev+1)
53
54      logical first
55      save first
56      data first/.true./
57c$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/
73c$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
104c  Correction d'un bug sauvage a verifier.
105c      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
130c 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
159c(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
Note: See TracBrowser for help on using the repository browser.