source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/yamada.F90 @ 2699

Last change on this file since 2699 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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