source: LMDZ6/branches/cirrus/libf/phylmd/yamada.F90 @ 5409

Last change on this file since 5409 was 4649, checked in by fhourdin, 16 months ago

Supression de depassements de tableaux

Routines non utilisees dans les configurations standard

  • 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
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  ! 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)
28
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
43
44
45  INTEGER nlay, nlev
46  ! ym      PARAMETER (nlay=klev)
47  ! ym      PARAMETER (nlev=klev+1)
48
49  LOGICAL first
50  SAVE first
51  DATA first/.TRUE./
52  !$OMP THREADPRIVATE(first)
53
54  INTEGER ig, k
55
56  REAL ri, zrif, zalpha, zsm
57  REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
58
59  REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
60  REAL l(klon, klev+1), l0(klon)
61
62  REAL sq(klon), sqz(klon), zz(klon, klev+1)
63  INTEGER iter
64
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)
69
70  REAL frif, falpha, fsm
71
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))
75
76  nlay = klev
77  nlev = klev + 1
78
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
98
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
124
125  ! iterration pour determiner la longueur de melange
126
127  DO ig = 1, ngrid
128    l0(ig) = 100.
129  END DO
130  DO k = 2, klev
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 k = 2, klev
143      DO ig = 1, ngrid
144        q2(ig, k) = l(ig, k)**2*zz(ig, k)
145        l(ig, k) = min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig, &
146          k)+l0(ig)), 0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.E-10)))
147        zq = sqrt(q2(ig,k))
148        sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
149        sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
150      END DO
151    END DO
152    DO ig = 1, ngrid
153      l0(ig) = 0.2*sqz(ig)/sq(ig)
154    END DO
155    ! (abd 3 5 2)         print*,'ITER=',iter,'  L0=',l0
156
157  END DO
158
159  DO k = 2, klev
160    DO ig = 1, ngrid
161      l(ig, k) = min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig, &
162        k)+l0(ig)), 0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.E-10)))
163      q2(ig, k) = l(ig, k)**2*zz(ig, k)
164      km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
165      kn(ig, k) = km(ig, k)*alpha(ig, k)
166    END DO
167  END DO
168
169  RETURN
170END SUBROUTINE yamada
Note: See TracBrowser for help on using the repository browser.