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

Last change on this file since 3557 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

  • Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
  • got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

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