/[lmdze]/trunk/Sources/phylmd/coefkz.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/coefkz.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 221 - (show annotations)
Thu Apr 20 14:44:47 2017 UTC (6 years, 11 months ago) by guez
File size: 8007 byte(s)
clcdrag is no longer used in LMDZ. Replaced by cdrag in LMDZ. In cdrag
in LMDZ, zxli is a symbolic constant, false. So removed case zxli true
in LMDZE.

read_sst is called zero (if no ocean point on the whole planet) time or
once per call of physiq. If mod(itap - 1, lmt_pas) == 0 then we have
advanced in time of lmt_pas and deja_lu is necessarily false.

qsat[sl] and dqsat[sl] were never called.

Added output of qsurf in histins, following LMDZ.

Last dummy argument dtime of phystokenc is always the same as first
dummy argument pdtphys, removed dtime.

Removed make rules for nag_xref95, since it does not exist any longer.

1 module coefkz_m
2
3 IMPLICIT none
4
5 contains
6
7 SUBROUTINE coefkz(nsrf, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, t, &
8 q, qsurf, coefm, coefh)
9
10 ! Authors: F. Hourdin, M. Forichon, Z. X. Li (LMD/CNRS)
11 ! Date: September 22nd, 1993
12 ! Objet : calculer le coefficient de frottement du sol ("Cdrag") et les
13 ! coefficients d'échange turbulent dans l'atmosphère.
14
15 use clcdrag_m, only: clcdrag
16 USE conf_phys_m, ONLY: iflag_pbl
17 USE dimphy, ONLY: klev, klon
18 USE fcttre, ONLY: foede, foeew
19 USE indicesol, ONLY: is_oce
20 USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
21 USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
22
23 integer, intent(in):: nsrf ! indicateur de la nature du sol
24
25 REAL, intent(in):: paprs(:, :) ! (klon, klev+1)
26 ! pression a chaque intercouche (en Pa)
27
28 real, intent(in):: pplay(:, :) ! (klon, klev)
29 ! pression au milieu de chaque couche (en Pa)
30
31 REAL, intent(in):: ksta, ksta_ter
32 REAL, intent(in):: ts(:) ! (knon) temperature du sol (en Kelvin)
33 REAL, intent(in):: rugos(:) ! (klon) longeur de rugosite (en m)
34 REAL, intent(in):: u(:, :), v(:, :) ! (klon, klev) wind
35 REAL, intent(in):: t(:, :) ! (klon, klev) temperature (K)
36 real, intent(in):: q(:, :) ! (klon, klev) vapeur d'eau (kg/kg)
37 real, intent(in):: qsurf(:) ! (knon)
38 REAL, intent(out):: coefm(:, :) ! (knon, klev) coefficient, vitesse
39
40 real, intent(out):: coefh(:, :) ! (knon, klev)
41 ! coefficient, chaleur et humidité
42
43 ! Local:
44
45 INTEGER knon ! nombre de points a traiter
46
47 INTEGER itop(size(coefm, 1))
48 ! (knon) numero de couche du sommet de la couche limite
49
50 ! Quelques constantes et options:
51
52 REAL, PARAMETER:: cepdu2 =0.1**2
53 REAL, PARAMETER:: CKAP = 0.4
54 REAL, PARAMETER:: cb = 5.
55 REAL, PARAMETER:: cc = 5.
56 REAL, PARAMETER:: cd = 5.
57 REAL, PARAMETER:: clam = 160.
58 REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
59
60 LOGICAL, PARAMETER:: richum = .TRUE.
61 ! utilise le nombre de Richardson humide
62
63 REAL, PARAMETER:: ric = 0.4 ! nombre de Richardson critique
64 REAL, PARAMETER:: prandtl = 0.4
65
66 REAL kstable ! diffusion minimale (situation stable)
67 REAL, PARAMETER:: mixlen = 35. ! constante contrôlant longueur de mélange
68 INTEGER, PARAMETER:: isommet = klev ! sommet de la couche limite
69
70 LOGICAL, PARAMETER:: tvirtu = .TRUE.
71 ! calculer Ri d'une maniere plus performante
72
73 LOGICAL, PARAMETER:: opt_ec = .FALSE.
74 ! formule du Centre Europeen dans l'atmosphere
75
76 INTEGER i, k
77 REAL zgeop(klon, klev)
78 REAL zmgeom(klon)
79 REAL ri(klon)
80 REAL l2(klon)
81
82 REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
83
84 REAL zdphi, zdu2, ztvd, ztvu, cdn
85 REAL scf
86 REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
87 logical zdelta
88 REAL z2geomf, zalh2, alm2, zscfh, scfm
89 REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
90
91 !--------------------------------------------------------------------
92
93 knon = size(coefm, 1)
94
95 ! Prescrire la valeur de contre-gradient
96 if (iflag_pbl.eq.1) then
97 DO k = 3, klev
98 gamt(k) = -1.0E-03
99 ENDDO
100 gamt(2) = -2.5E-03
101 else
102 DO k = 2, klev
103 gamt(k) = 0.0
104 ENDDO
105 ENDIF
106
107 IF ( nsrf .NE. is_oce ) THEN
108 kstable = ksta_ter
109 ELSE
110 kstable = ksta
111 ENDIF
112
113 ! Calculer les géopotentiels de chaque couche
114 DO i = 1, knon
115 zgeop(i, 1) = RD * t(i, 1) / (0.5 * (paprs(i, 1) + pplay(i, 1))) &
116 * (paprs(i, 1) - pplay(i, 1))
117 ENDDO
118 DO k = 2, klev
119 DO i = 1, knon
120 zgeop(i, k) = zgeop(i, k-1) &
121 + RD * 0.5*(t(i, k-1)+t(i, k)) / paprs(i, k) &
122 * (pplay(i, k-1)-pplay(i, k))
123 ENDDO
124 ENDDO
125
126 ! Calculer le frottement au sol (Cdrag)
127
128 DO i = 1, knon
129 u1(i) = u(i, 1)
130 v1(i) = v(i, 1)
131 t1(i) = t(i, 1)
132 q1(i) = q(i, 1)
133 z1(i) = zgeop(i, 1)
134 ENDDO
135
136 CALL clcdrag(nsrf, u1, v1, t1, q1, z1, ts, qsurf, rugos, coefm(:, 1), &
137 coefh(:, 1))
138
139 ! Calculer les coefficients turbulents dans l'atmosphere
140
141 itop = isommet
142
143 loop_vertical: DO k = 2, isommet
144 loop_horiz: DO i = 1, knon
145 zdu2 = MAX(cepdu2, (u(i, k)-u(i, k-1))**2 &
146 +(v(i, k)-v(i, k-1))**2)
147 zmgeom(i) = zgeop(i, k)-zgeop(i, k-1)
148 zdphi =zmgeom(i) / 2.0
149 zt = (t(i, k)+t(i, k-1)) * 0.5
150 zq = (q(i, k)+q(i, k-1)) * 0.5
151
152 ! calculer Qs et dQs/dT:
153
154 zdelta = RTT >=zt
155 zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
156 / (1. + RVTMP2*zq)
157 zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
158 zqs = MIN(0.5, zqs)
159 zcor = 1./(1.-RETV*zqs)
160 zqs = zqs*zcor
161 zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
162
163 ! calculer la fraction nuageuse (processus humide):
164
165 zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
166 zfr = MAX(0.0, MIN(1.0, zfr))
167 IF (.NOT.richum) zfr = 0.0
168
169 ! calculer le nombre de Richardson:
170
171 IF (tvirtu) THEN
172 ztvd =( t(i, k) &
173 + zdphi/RCPD/(1.+RVTMP2*zq) &
174 *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
175 )*(1.+RETV*q(i, k))
176 ztvu =( t(i, k-1) &
177 - zdphi/RCPD/(1.+RVTMP2*zq) &
178 *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) &
179 )*(1.+RETV*q(i, k-1))
180 ri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
181 ri(i) = ri(i) &
182 + zmgeom(i)*zmgeom(i)/RG*gamt(k) &
183 *(paprs(i, k)/101325.0)**RKAPPA &
184 /(zdu2*0.5*(ztvd+ztvu))
185 ELSE
186 ! calcul de Ridchardson compatible LMD5
187 ri(i) =(RCPD*(t(i, k)-t(i, k-1)) &
188 -RD*0.5*(t(i, k)+t(i, k-1))/paprs(i, k) &
189 *(pplay(i, k)-pplay(i, k-1)) &
190 )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i, k-1)+t(i, k)))
191 ri(i) = ri(i) + &
192 zmgeom(i)*zmgeom(i)*gamt(k)/RG &
193 *(paprs(i, k)/101325.0)**RKAPPA &
194 /(zdu2*0.5*(t(i, k-1)+t(i, k)))
195 ENDIF
196
197 ! finalement, les coefficients d'echange sont obtenus:
198
199 cdn = SQRT(zdu2) / zmgeom(i) * RG
200
201 IF (opt_ec) THEN
202 z2geomf = zgeop(i, k-1)+zgeop(i, k)
203 alm2 = (0.5*ckap/RG*z2geomf &
204 /(1.+0.5*ckap/rg/clam*z2geomf))**2
205 zalh2 = (0.5*ckap/rg*z2geomf &
206 /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
207 IF (ri(i) < 0.) THEN
208 ! situation instable
209 scf = ((zgeop(i, k)/zgeop(i, k-1))**(1./3.)-1.)**3 &
210 / (zmgeom(i)/RG)**3 / (zgeop(i, k-1)/RG)
211 scf = SQRT(-ri(i)*scf)
212 scfm = 1.0 / (1.0+3.0*cb*cc*alm2*scf)
213 zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*scf)
214 coefm(i, k) = cdn * alm2 * (1. - 2. * cb * ri(i) * scfm)
215 coefh(i, k) = cdn*zalh2*(1.-3.0*cb*ri(i)*zscfh)
216 ELSE
217 ! situation stable
218 scf = SQRT(1.+cd*ri(i))
219 coefm(i, k) = cdn * alm2 / (1. + 2. * cb * ri(i) / scf)
220 coefh(i, k) = cdn*zalh2/(1.+3.0*cb*ri(i)*scf)
221 ENDIF
222 ELSE
223 l2(i) = (mixlen*MAX(0.0, (paprs(i, k)-paprs(i, itop(i)+1)) &
224 /(paprs(i, 2)-paprs(i, itop(i)+1)) ))**2
225 coefm(i, k) = sqrt(max(cdn**2 * (ric - ri(i)) / ric, kstable))
226 coefm(i, k)= l2(i) * coefm(i, k)
227 coefh(i, k) = coefm(i, k) / prandtl ! h et m different
228 ENDIF
229 ENDDO loop_horiz
230 ENDDO loop_vertical
231
232 ! Au-delà du sommet, pas de diffusion turbulente :
233 forall (i = 1: knon)
234 coefh(i, itop(i) + 1:) = 0.
235 coefm(i, itop(i) + 1:) = 0.
236 END forall
237
238 END SUBROUTINE coefkz
239
240 end module coefkz_m

  ViewVC Help
Powered by ViewVC 1.1.21