source: LMDZ5/branches/testing/libf/phylmd/yamada.F90 @ 3046

Last change on this file since 3046 was 2408, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2298:2396 into testing branch

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