source: LMDZ5/trunk/libf/phylmd/yamada.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • 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.