1 |
|
|
! |
2 |
|
|
MODULE coef_diff_turb_mod |
3 |
|
|
! |
4 |
|
|
! This module contains some procedures for calculation of the coefficients of the |
5 |
|
|
! turbulent diffusion in the atmosphere and coefficients for turbulent diffusion |
6 |
|
|
! at surface(cdrag) |
7 |
|
|
! |
8 |
|
|
IMPLICIT NONE |
9 |
|
|
|
10 |
|
|
CONTAINS |
11 |
|
|
! |
12 |
|
|
!**************************************************************************************** |
13 |
|
|
! |
14 |
|
1152 |
SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, & |
15 |
|
1152 |
ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, & |
16 |
|
|
ycoefm, ycoefh ,yq2, ydrgpro) |
17 |
|
|
|
18 |
|
|
USE dimphy |
19 |
|
|
USE indice_sol_mod |
20 |
|
|
USE print_control_mod, ONLY: prt_level, lunout |
21 |
|
|
! |
22 |
|
|
! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the |
23 |
|
|
! atmosphere |
24 |
|
|
! NB! No values are calculated between surface and the first model layer. |
25 |
|
|
! ycoefm(:,1) and ycoefh(:,1) are not valid !!! |
26 |
|
|
! |
27 |
|
|
! |
28 |
|
|
! Input arguments |
29 |
|
|
!**************************************************************************************** |
30 |
|
|
REAL, INTENT(IN) :: dtime |
31 |
|
|
INTEGER, INTENT(IN) :: nsrf, knon |
32 |
|
|
INTEGER, DIMENSION(klon), INTENT(IN) :: ni |
33 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: ypaprs |
34 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: ypplay |
35 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: yu, yv |
36 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: yq, yt |
37 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: yts, yqsurf |
38 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: ycdragm |
39 |
|
|
!FC |
40 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: ydrgpro |
41 |
|
|
|
42 |
|
|
|
43 |
|
|
! InOutput arguments |
44 |
|
|
!**************************************************************************************** |
45 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(INOUT):: yq2 |
46 |
|
|
|
47 |
|
|
! Output arguments |
48 |
|
|
!**************************************************************************************** |
49 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefh |
50 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefm |
51 |
|
|
|
52 |
|
|
! Other local variables |
53 |
|
|
!**************************************************************************************** |
54 |
|
|
INTEGER :: k, i, j |
55 |
|
2304 |
REAL, DIMENSION(klon,klev) :: ycoefm0, ycoefh0, yzlay, yteta |
56 |
|
2304 |
REAL, DIMENSION(klon,klev+1) :: yzlev, q2diag, ykmm, ykmn, ykmq |
57 |
|
1152 |
REAL, DIMENSION(klon) :: yustar |
58 |
|
|
|
59 |
|
|
! Include |
60 |
|
|
!**************************************************************************************** |
61 |
|
|
INCLUDE "clesphys.h" |
62 |
|
|
INCLUDE "compbl.h" |
63 |
|
|
INCLUDE "YOETHF.h" |
64 |
|
|
INCLUDE "YOMCST.h" |
65 |
|
|
|
66 |
|
|
|
67 |
✓✓✓✓
|
45850752 |
ykmm = 0 !ym missing init |
68 |
✓✓✓✓
|
45850752 |
ykmn = 0 !ym missing init |
69 |
✓✓✓✓
|
45850752 |
ykmq = 0 !ym missing init |
70 |
|
|
|
71 |
|
|
|
72 |
|
|
!**************************************************************************************** |
73 |
|
|
! Calcul de coefficients de diffusion turbulent de l'atmosphere : |
74 |
|
|
! ycoefm(:,2:klev), ycoefh(:,2:klev) |
75 |
|
|
! |
76 |
|
|
!**************************************************************************************** |
77 |
|
|
|
78 |
|
|
CALL coefkz(nsrf, knon, ypaprs, ypplay, & |
79 |
|
|
ksta, ksta_ter, & |
80 |
|
|
yts, yu, yv, yt, yq, & |
81 |
|
|
yqsurf, & |
82 |
|
1152 |
ycoefm, ycoefh) |
83 |
|
|
|
84 |
|
|
!**************************************************************************************** |
85 |
|
|
! Eventuelle recalcule des coeffeicients de diffusion turbulent de l'atmosphere : |
86 |
|
|
! ycoefm(:,2:klev), ycoefh(:,2:klev) |
87 |
|
|
! |
88 |
|
|
!**************************************************************************************** |
89 |
|
|
|
90 |
|
1152 |
IF (iflag_pbl.EQ.1) THEN |
91 |
|
|
CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, & |
92 |
|
|
ycoefm0, ycoefh0) |
93 |
|
|
|
94 |
|
|
DO k = 2, klev |
95 |
|
|
DO i = 1, knon |
96 |
|
|
ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k)) |
97 |
|
|
ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k)) |
98 |
|
|
ENDDO |
99 |
|
|
ENDDO |
100 |
|
|
ENDIF |
101 |
|
|
|
102 |
|
|
|
103 |
|
|
!**************************************************************************************** |
104 |
|
|
! Calcul d'une diffusion minimale pour les conditions tres stables |
105 |
|
|
! |
106 |
|
|
!**************************************************************************************** |
107 |
|
1152 |
IF (ok_kzmin) THEN |
108 |
|
|
CALL coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycdragm, & |
109 |
|
|
ycoefm0,ycoefh0) |
110 |
|
|
|
111 |
|
|
DO k = 2, klev |
112 |
|
|
DO i = 1, knon |
113 |
|
|
ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k)) |
114 |
|
|
ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k)) |
115 |
|
|
ENDDO |
116 |
|
|
ENDDO |
117 |
|
|
|
118 |
|
|
ENDIF |
119 |
|
|
|
120 |
|
|
|
121 |
|
|
!**************************************************************************************** |
122 |
|
|
! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin |
123 |
|
|
! |
124 |
|
|
!**************************************************************************************** |
125 |
|
|
|
126 |
|
1152 |
IF (iflag_pbl.GE.3) THEN |
127 |
|
|
|
128 |
|
|
yzlay(1:knon,1)= & |
129 |
|
|
RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1))) & |
130 |
✓✓ |
473954 |
*(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG |
131 |
✓✓ |
44928 |
DO k=2,klev |
132 |
✓✓ |
18011404 |
DO i = 1, knon |
133 |
|
|
yzlay(i,k)= & |
134 |
|
|
yzlay(i,k-1)+RD*0.5*(yt(i,k-1)+yt(i,k)) & |
135 |
|
18010252 |
/ypaprs(i,k)*(ypplay(i,k-1)-ypplay(i,k))/RG |
136 |
|
|
END DO |
137 |
|
|
END DO |
138 |
|
|
|
139 |
✓✓ |
46080 |
DO k=1,klev |
140 |
✓✓ |
18485358 |
DO i = 1, knon |
141 |
|
|
yteta(i,k)= & |
142 |
|
|
yt(i,k)*(ypaprs(i,1)/ypplay(i,k))**RKAPPA & |
143 |
|
44928 |
*(1.+0.61*yq(i,k)) |
144 |
|
|
END DO |
145 |
|
|
END DO |
146 |
|
|
|
147 |
✓✓ |
473954 |
yzlev(1:knon,1)=0. |
148 |
✓✓ |
473954 |
yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1) |
149 |
✓✓ |
44928 |
DO k=2,klev |
150 |
✓✓ |
18011404 |
DO i = 1, knon |
151 |
|
18010252 |
yzlev(i,k)=0.5*(yzlay(i,k)+yzlay(i,k-1)) |
152 |
|
|
END DO |
153 |
|
|
END DO |
154 |
|
|
|
155 |
|
|
!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
156 |
|
|
!!$! Pour memoire, le papier Hourdin et al. 2002 a ete obtenur avec un |
157 |
|
|
!!$! bug sur les coefficients de surface : |
158 |
|
|
!!$! ycdragh(1:knon) = ycoefm(1:knon,1) |
159 |
|
|
!!$! ycdragm(1:knon) = ycoefh(1:knon,1) |
160 |
|
|
!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
161 |
|
1152 |
CALL ustarhb(knon,yu,yv,ycdragm, yustar) |
162 |
|
|
|
163 |
✗✓ |
1152 |
IF (prt_level > 9) THEN |
164 |
|
|
WRITE(lunout,*) 'USTAR = ',(yustar(i),i=1,knon) |
165 |
|
|
ENDIF |
166 |
|
|
|
167 |
|
|
! iflag_pbl peut etre utilise comme longuer de melange |
168 |
|
1152 |
IF (iflag_pbl.GE.31) THEN |
169 |
|
|
CALL vdif_kcay(knon,dtime,RG,RD,ypaprs,yt, & |
170 |
|
|
yzlev,yzlay,yu,yv,yteta, & |
171 |
|
|
ycdragm,yq2,q2diag,ykmm,ykmn,yustar, & |
172 |
|
|
iflag_pbl) |
173 |
✓✗ |
1152 |
ELSE IF (iflag_pbl<20) THEN |
174 |
|
|
CALL yamada4(ni,nsrf,knon,dtime,RG,RD,ypaprs,yt, & |
175 |
|
|
yzlev,yzlay,yu,yv,yteta, & |
176 |
|
|
ycdragm,yq2,ykmm,ykmn,ykmq,yustar, & |
177 |
|
1152 |
iflag_pbl,ydrgpro) |
178 |
|
|
!FC |
179 |
|
|
ENDIF |
180 |
|
|
|
181 |
✓✓✓✓
|
18011404 |
ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev) |
182 |
✓✓✓✓
|
18011404 |
ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev) |
183 |
|
|
|
184 |
|
|
ELSE |
185 |
|
|
! No TKE for Standard Physics |
186 |
|
|
yq2=0. |
187 |
|
|
ENDIF !(iflag_pbl.ge.3) |
188 |
|
|
|
189 |
|
1152 |
END SUBROUTINE coef_diff_turb |
190 |
|
|
! |
191 |
|
|
!**************************************************************************************** |
192 |
|
|
! |
193 |
|
1152 |
SUBROUTINE coefkz(nsrf, knon, paprs, pplay, & |
194 |
|
|
ksta, ksta_ter, & |
195 |
|
|
ts, & |
196 |
|
|
u,v,t,q, & |
197 |
|
|
qsurf, & |
198 |
|
1152 |
pcfm, pcfh) |
199 |
|
|
|
200 |
✗✓✗✓
✓✗✗✓
|
18441582 |
USE dimphy |
201 |
|
|
USE indice_sol_mod |
202 |
|
|
USE print_control_mod, ONLY: prt_level, lunout |
203 |
|
|
|
204 |
|
|
!====================================================================== |
205 |
|
|
! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922 |
206 |
|
|
! (une version strictement identique a l'ancien modele) |
207 |
|
|
! Objet: calculer le coefficient du frottement du sol (Cdrag) et les |
208 |
|
|
! coefficients d'echange turbulent dans l'atmosphere. |
209 |
|
|
! Arguments: |
210 |
|
|
! nsrf-----input-I- indicateur de la nature du sol |
211 |
|
|
! knon-----input-I- nombre de points a traiter |
212 |
|
|
! paprs----input-R- pregssion a chaque intercouche (en Pa) |
213 |
|
|
! pplay----input-R- pression au milieu de chaque couche (en Pa) |
214 |
|
|
! ts-------input-R- temperature du sol (en Kelvin) |
215 |
|
|
! u--------input-R- vitesse u |
216 |
|
|
! v--------input-R- vitesse v |
217 |
|
|
! t--------input-R- temperature (K) |
218 |
|
|
! q--------input-R- vapeur d'eau (kg/kg) |
219 |
|
|
! |
220 |
|
|
! pcfm-----output-R- coefficients a calculer (vitesse) |
221 |
|
|
! pcfh-----output-R- coefficients a calculer (chaleur et humidite) |
222 |
|
|
!====================================================================== |
223 |
|
|
INCLUDE "YOETHF.h" |
224 |
|
|
INCLUDE "YOMCST.h" |
225 |
|
|
INCLUDE "FCTTRE.h" |
226 |
|
|
INCLUDE "compbl.h" |
227 |
|
|
! |
228 |
|
|
! Arguments: |
229 |
|
|
! |
230 |
|
|
INTEGER, INTENT(IN) :: knon, nsrf |
231 |
|
|
REAL, INTENT(IN) :: ksta, ksta_ter |
232 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: ts |
233 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs |
234 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay |
235 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: u, v, t, q |
236 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: qsurf |
237 |
|
|
|
238 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: pcfm, pcfh |
239 |
|
|
|
240 |
|
|
! |
241 |
|
|
! Local variables: |
242 |
|
|
! |
243 |
|
2304 |
INTEGER, DIMENSION(klon) :: itop ! numero de couche du sommet de la couche limite |
244 |
|
|
! |
245 |
|
|
! Quelques constantes et options: |
246 |
|
|
! |
247 |
|
|
REAL, PARAMETER :: cepdu2=0.1**2 |
248 |
|
|
REAL, PARAMETER :: CKAP=0.4 |
249 |
|
|
REAL, PARAMETER :: cb=5.0 |
250 |
|
|
REAL, PARAMETER :: cc=5.0 |
251 |
|
|
REAL, PARAMETER :: cd=5.0 |
252 |
|
|
REAL, PARAMETER :: clam=160.0 |
253 |
|
|
REAL, PARAMETER :: ratqs=0.05 ! largeur de distribution de vapeur d'eau |
254 |
|
|
LOGICAL, PARAMETER :: richum=.TRUE. ! utilise le nombre de Richardson humide |
255 |
|
|
REAL, PARAMETER :: ric=0.4 ! nombre de Richardson critique |
256 |
|
|
REAL, PARAMETER :: prandtl=0.4 |
257 |
|
|
REAL kstable ! diffusion minimale (situation stable) |
258 |
|
|
! GKtest |
259 |
|
|
! PARAMETER (kstable=1.0e-10) |
260 |
|
|
!IM: 261103 REAL kstable_ter, kstable_sinon |
261 |
|
|
!IM: 211003 cf GK PARAMETER (kstable_ter = 1.0e-6) |
262 |
|
|
!IM: 261103 PARAMETER (kstable_ter = 1.0e-8) |
263 |
|
|
!IM: 261103 PARAMETER (kstable_ter = 1.0e-10) |
264 |
|
|
!IM: 261103 PARAMETER (kstable_sinon = 1.0e-10) |
265 |
|
|
! fin GKtest |
266 |
|
|
REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange |
267 |
|
|
INTEGER isommet ! le sommet de la couche limite |
268 |
|
|
LOGICAL, PARAMETER :: tvirtu=.TRUE. ! calculer Ri d'une maniere plus performante |
269 |
|
|
LOGICAL, PARAMETER :: opt_ec=.FALSE.! formule du Centre Europeen dans l'atmosphere |
270 |
|
|
|
271 |
|
|
! |
272 |
|
|
! Variables locales: |
273 |
|
|
INTEGER i, k !IM 120704 |
274 |
|
2304 |
REAL zgeop(klon,klev) |
275 |
|
2304 |
REAL zmgeom(klon) |
276 |
|
2304 |
REAL zri(klon) |
277 |
|
2304 |
REAL zl2(klon) |
278 |
|
|
REAL zdphi, zdu2, ztvd, ztvu, zcdn |
279 |
|
|
REAL zscf |
280 |
|
|
REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs |
281 |
|
|
REAL z2geomf, zalh2, zalm2, zscfh, zscfm |
282 |
|
|
REAL, PARAMETER :: t_coup=273.15 |
283 |
|
|
LOGICAL, PARAMETER :: check=.FALSE. |
284 |
|
|
! |
285 |
|
|
! contre-gradient pour la chaleur sensible: Kelvin/metre |
286 |
|
2304 |
REAL gamt(2:klev) |
287 |
|
|
|
288 |
|
|
LOGICAL, SAVE :: appel1er=.TRUE. |
289 |
|
|
!$OMP THREADPRIVATE(appel1er) |
290 |
|
|
! |
291 |
|
|
! Fonctions thermodynamiques et fonctions d'instabilite |
292 |
|
|
REAL fsta, fins, x |
293 |
|
|
|
294 |
|
|
fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x)) |
295 |
|
|
fins(x) = SQRT(1.0-18.0*x) |
296 |
|
|
|
297 |
|
1152 |
isommet=klev |
298 |
|
|
|
299 |
✓✗ |
1152 |
IF (appel1er) THEN |
300 |
✗✓ |
1152 |
IF (prt_level > 9) THEN |
301 |
|
|
WRITE(lunout,*)'coefkz, opt_ec:', opt_ec |
302 |
|
|
WRITE(lunout,*)'coefkz, richum:', richum |
303 |
|
|
IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs |
304 |
|
|
WRITE(lunout,*)'coefkz, isommet:', isommet |
305 |
|
|
WRITE(lunout,*)'coefkz, tvirtu:', tvirtu |
306 |
|
|
appel1er = .FALSE. |
307 |
|
|
ENDIF |
308 |
|
|
ENDIF |
309 |
|
|
! |
310 |
|
|
! Initialiser les sorties |
311 |
|
|
! |
312 |
✓✓ |
46080 |
DO k = 1, klev |
313 |
✓✓ |
18485358 |
DO i = 1, knon |
314 |
|
18439278 |
pcfm(i,k) = 0.0 |
315 |
|
18484206 |
pcfh(i,k) = 0.0 |
316 |
|
|
ENDDO |
317 |
|
|
ENDDO |
318 |
✓✓ |
473954 |
DO i = 1, knon |
319 |
|
473954 |
itop(i) = 0 |
320 |
|
|
ENDDO |
321 |
|
|
|
322 |
|
|
! |
323 |
|
|
! Prescrire la valeur de contre-gradient |
324 |
|
|
! |
325 |
|
1152 |
IF (iflag_pbl.EQ.1) THEN |
326 |
|
|
DO k = 3, klev |
327 |
|
|
gamt(k) = -1.0E-03 |
328 |
|
|
ENDDO |
329 |
|
|
gamt(2) = -2.5E-03 |
330 |
|
|
ELSE |
331 |
✓✓ |
44928 |
DO k = 2, klev |
332 |
|
44928 |
gamt(k) = 0.0 |
333 |
|
|
ENDDO |
334 |
|
|
ENDIF |
335 |
|
|
!IM cf JLD/ GKtest |
336 |
✓✓ |
1152 |
IF ( nsrf .NE. is_oce ) THEN |
337 |
|
|
!IM 261103 kstable = kstable_ter |
338 |
|
864 |
kstable = ksta_ter |
339 |
|
|
ELSE |
340 |
|
|
!IM 261103 kstable = kstable_sinon |
341 |
|
288 |
kstable = ksta |
342 |
|
|
ENDIF |
343 |
|
|
!IM cf JLD/ GKtest fin |
344 |
|
|
|
345 |
|
|
! |
346 |
|
|
! Calculer les geopotentiels de chaque couche |
347 |
|
|
! |
348 |
✓✓ |
473954 |
DO i = 1, knon |
349 |
|
|
zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) & |
350 |
|
1152 |
* (paprs(i,1)-pplay(i,1)) |
351 |
|
|
ENDDO |
352 |
✓✓ |
44928 |
DO k = 2, klev |
353 |
✓✓ |
18011404 |
DO i = 1, knon |
354 |
|
|
zgeop(i,k) = zgeop(i,k-1) & |
355 |
|
|
+ RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) & |
356 |
|
43776 |
* (pplay(i,k-1)-pplay(i,k)) |
357 |
|
|
ENDDO |
358 |
|
|
ENDDO |
359 |
|
|
|
360 |
|
|
! |
361 |
|
|
! Calculer les coefficients turbulents dans l'atmosphere |
362 |
|
|
! |
363 |
✓✓ |
473954 |
DO i = 1, knon |
364 |
|
473954 |
itop(i) = isommet |
365 |
|
|
ENDDO |
366 |
|
|
|
367 |
|
|
|
368 |
✓✓ |
44928 |
DO k = 2, isommet |
369 |
✓✓ |
18011404 |
DO i = 1, knon |
370 |
|
|
zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 & |
371 |
|
17966476 |
+(v(i,k)-v(i,k-1))**2) |
372 |
|
17966476 |
zmgeom(i)=zgeop(i,k)-zgeop(i,k-1) |
373 |
|
17966476 |
zdphi =zmgeom(i) / 2.0 |
374 |
|
17966476 |
zt = (t(i,k)+t(i,k-1)) * 0.5 |
375 |
|
17966476 |
zq = (q(i,k)+q(i,k-1)) * 0.5 |
376 |
|
|
|
377 |
|
|
! |
378 |
|
|
! Calculer Qs et dQs/dT: |
379 |
|
|
! |
380 |
|
|
IF (thermcep) THEN |
381 |
|
17966476 |
zdelta = MAX(0.,SIGN(1.,RTT-zt)) |
382 |
|
|
zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta) & |
383 |
|
17966476 |
+ R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta |
384 |
|
17966476 |
zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k) |
385 |
|
17966476 |
zqs = MIN(0.5,zqs) |
386 |
|
17966476 |
zcor = 1./(1.-RETV*zqs) |
387 |
|
17966476 |
zqs = zqs*zcor |
388 |
|
17966476 |
zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor) |
389 |
|
|
ELSE |
390 |
|
|
IF (zt .LT. t_coup) THEN |
391 |
|
|
zqs = qsats(zt) / pplay(i,k) |
392 |
|
|
zdqs = dqsats(zt,zqs) |
393 |
|
|
ELSE |
394 |
|
|
zqs = qsatl(zt) / pplay(i,k) |
395 |
|
|
zdqs = dqsatl(zt,zqs) |
396 |
|
|
ENDIF |
397 |
|
|
ENDIF |
398 |
|
|
! |
399 |
|
|
! calculer la fraction nuageuse (processus humide): |
400 |
|
|
! |
401 |
|
17966476 |
if (zq /= 0.) then |
402 |
|
17966476 |
zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq) |
403 |
|
|
else |
404 |
|
|
zfr = 0. |
405 |
|
|
end if |
406 |
|
17966476 |
zfr = MAX(0.0,MIN(1.0,zfr)) |
407 |
|
|
IF (.NOT.richum) zfr = 0.0 |
408 |
|
|
! |
409 |
|
|
! calculer le nombre de Richardson: |
410 |
|
|
! |
411 |
|
|
IF (tvirtu) THEN |
412 |
|
|
ztvd =( t(i,k) & |
413 |
|
|
+ zdphi/RCPD/(1.+RVTMP2*zq) & |
414 |
|
|
*( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) & |
415 |
|
17966476 |
)*(1.+RETV*q(i,k)) |
416 |
|
|
ztvu =( t(i,k-1) & |
417 |
|
|
- zdphi/RCPD/(1.+RVTMP2*zq) & |
418 |
|
|
*( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) & |
419 |
|
17966476 |
)*(1.+RETV*q(i,k-1)) |
420 |
|
17966476 |
zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu)) |
421 |
|
|
zri(i) = zri(i) & |
422 |
|
|
+ zmgeom(i)*zmgeom(i)/RG*gamt(k) & |
423 |
|
|
*(paprs(i,k)/101325.0)**RKAPPA & |
424 |
|
17966476 |
/(zdu2*0.5*(ztvd+ztvu)) |
425 |
|
|
|
426 |
|
|
ELSE ! calcul de Ridchardson compatible LMD5 |
427 |
|
|
|
428 |
|
|
zri(i) =(RCPD*(t(i,k)-t(i,k-1)) & |
429 |
|
|
-RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k) & |
430 |
|
|
*(pplay(i,k)-pplay(i,k-1)) & |
431 |
|
|
)*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k))) |
432 |
|
|
zri(i) = zri(i) + & |
433 |
|
|
zmgeom(i)*zmgeom(i)*gamt(k)/RG & |
434 |
|
|
*(paprs(i,k)/101325.0)**RKAPPA & |
435 |
|
|
/(zdu2*0.5*(t(i,k-1)+t(i,k))) |
436 |
|
|
ENDIF |
437 |
|
|
! |
438 |
|
|
! finalement, les coefficients d'echange sont obtenus: |
439 |
|
|
! |
440 |
|
17966476 |
zcdn=SQRT(zdu2) / zmgeom(i) * RG |
441 |
|
|
|
442 |
|
43776 |
IF (opt_ec) THEN |
443 |
|
|
z2geomf=zgeop(i,k-1)+zgeop(i,k) |
444 |
|
|
zalm2=(0.5*ckap/RG*z2geomf & |
445 |
|
|
/(1.+0.5*ckap/rg/clam*z2geomf))**2 |
446 |
|
|
zalh2=(0.5*ckap/rg*z2geomf & |
447 |
|
|
/(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2 |
448 |
|
|
IF (zri(i).LT.0.0) THEN ! situation instable |
449 |
|
|
zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3 & |
450 |
|
|
/ (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG) |
451 |
|
|
zscf = SQRT(-zri(i)*zscf) |
452 |
|
|
zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf) |
453 |
|
|
zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf) |
454 |
|
|
pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm) |
455 |
|
|
pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh) |
456 |
|
|
ELSE ! situation stable |
457 |
|
|
zscf=SQRT(1.+cd*zri(i)) |
458 |
|
|
pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf) |
459 |
|
|
pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf) |
460 |
|
|
ENDIF |
461 |
|
|
ELSE |
462 |
|
|
zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) & |
463 |
|
17966476 |
/(paprs(i,2)-paprs(i,itop(i)+1)) ))**2 |
464 |
|
17966476 |
pcfm(i,k)=SQRT(MAX(zcdn*zcdn*(ric-zri(i))/ric, kstable)) |
465 |
|
17966476 |
pcfm(i,k)= zl2(i)* pcfm(i,k) |
466 |
|
17966476 |
pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different |
467 |
|
|
ENDIF |
468 |
|
|
ENDDO |
469 |
|
|
ENDDO |
470 |
|
|
|
471 |
|
|
! |
472 |
|
|
! Au-dela du sommet, pas de diffusion turbulente: |
473 |
|
|
! |
474 |
✓✓ |
473954 |
DO i = 1, knon |
475 |
✗✓ |
473954 |
IF (itop(i)+1 .LE. klev) THEN |
476 |
|
|
DO k = itop(i)+1, klev |
477 |
|
|
pcfh(i,k) = 0.0 |
478 |
|
|
pcfm(i,k) = 0.0 |
479 |
|
|
ENDDO |
480 |
|
|
ENDIF |
481 |
|
|
ENDDO |
482 |
|
|
|
483 |
|
1152 |
END SUBROUTINE coefkz |
484 |
|
|
! |
485 |
|
|
!**************************************************************************************** |
486 |
|
|
! |
487 |
|
|
SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, & |
488 |
|
|
pcfm, pcfh) |
489 |
|
|
|
490 |
✗✓✓✗
|
54373382 |
USE dimphy |
491 |
|
|
USE indice_sol_mod |
492 |
|
|
|
493 |
|
|
!====================================================================== |
494 |
|
|
! J'introduit un peu de diffusion sauf dans les endroits |
495 |
|
|
! ou une forte inversion est presente |
496 |
|
|
! On peut dire qu'il represente la convection peu profonde |
497 |
|
|
! |
498 |
|
|
! Arguments: |
499 |
|
|
! nsrf-----input-I- indicateur de la nature du sol |
500 |
|
|
! knon-----input-I- nombre de points a traiter |
501 |
|
|
! paprs----input-R- pression a chaque intercouche (en Pa) |
502 |
|
|
! pplay----input-R- pression au milieu de chaque couche (en Pa) |
503 |
|
|
! t--------input-R- temperature (K) |
504 |
|
|
! |
505 |
|
|
! pcfm-----output-R- coefficients a calculer (vitesse) |
506 |
|
|
! pcfh-----output-R- coefficients a calculer (chaleur et humidite) |
507 |
|
|
!====================================================================== |
508 |
|
|
! |
509 |
|
|
! Arguments: |
510 |
|
|
! |
511 |
|
|
INTEGER, INTENT(IN) :: knon, nsrf |
512 |
|
|
REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs |
513 |
|
|
REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay |
514 |
|
|
REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon,klev) |
515 |
|
|
|
516 |
|
|
REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh |
517 |
|
|
! |
518 |
|
|
! Quelques constantes et options: |
519 |
|
|
! |
520 |
|
|
REAL, PARAMETER :: prandtl=0.4 |
521 |
|
|
REAL, PARAMETER :: kstable=0.002 |
522 |
|
|
! REAL, PARAMETER :: kstable=0.001 |
523 |
|
|
REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange |
524 |
|
|
REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible |
525 |
|
|
! PARAMETER (seuil=-0.04) |
526 |
|
|
! PARAMETER (seuil=-0.06) |
527 |
|
|
! PARAMETER (seuil=-0.09) |
528 |
|
|
|
529 |
|
|
! |
530 |
|
|
! Variables locales: |
531 |
|
|
! |
532 |
|
|
INTEGER i, k, invb(knon) |
533 |
|
|
REAL zl2(knon) |
534 |
|
|
REAL zdthmin(knon), zdthdp |
535 |
|
|
|
536 |
|
|
INCLUDE "YOMCST.h" |
537 |
|
|
! |
538 |
|
|
! Initialiser les sorties |
539 |
|
|
! |
540 |
|
|
DO k = 1, klev |
541 |
|
|
DO i = 1, knon |
542 |
|
|
pcfm(i,k) = 0.0 |
543 |
|
|
pcfh(i,k) = 0.0 |
544 |
|
|
ENDDO |
545 |
|
|
ENDDO |
546 |
|
|
|
547 |
|
|
! |
548 |
|
|
! Chercher la zone d'inversion forte |
549 |
|
|
! |
550 |
|
|
DO i = 1, knon |
551 |
|
|
invb(i) = klev |
552 |
|
|
zdthmin(i)=0.0 |
553 |
|
|
ENDDO |
554 |
|
|
DO k = 2, klev/2-1 |
555 |
|
|
DO i = 1, knon |
556 |
|
|
zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) & |
557 |
|
|
- RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1) |
558 |
|
|
zdthdp = zdthdp * 100.0 |
559 |
|
|
IF (pplay(i,k).GT.0.8*paprs(i,1) .AND. & |
560 |
|
|
zdthdp.LT.zdthmin(i) ) THEN |
561 |
|
|
zdthmin(i) = zdthdp |
562 |
|
|
invb(i) = k |
563 |
|
|
ENDIF |
564 |
|
|
ENDDO |
565 |
|
|
ENDDO |
566 |
|
|
|
567 |
|
|
! |
568 |
|
|
! Introduire une diffusion: |
569 |
|
|
! |
570 |
|
|
IF ( nsrf.EQ.is_oce ) THEN |
571 |
|
|
DO k = 2, klev |
572 |
|
|
DO i = 1, knon |
573 |
|
|
!IM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean |
574 |
|
|
!IM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion |
575 |
|
|
!IM cf JLD/ GKtest TERkz2 |
576 |
|
|
! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre |
577 |
|
|
! fin GKtest |
578 |
|
|
|
579 |
|
|
|
580 |
|
|
! s'il n'y a pas d'inversion ou si l'inversion est trop faible |
581 |
|
|
! IF ( (nsrf.EQ.is_oce) .AND. & |
582 |
|
|
IF ( (invb(i).EQ.klev) .OR. (zdthmin(i).GT.seuil) ) THEN |
583 |
|
|
zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1)) & |
584 |
|
|
/(paprs(i,2)-paprs(i,klev+1)) ))**2 |
585 |
|
|
pcfm(i,k)= zl2(i)* kstable |
586 |
|
|
pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different |
587 |
|
|
ENDIF |
588 |
|
|
ENDDO |
589 |
|
|
ENDDO |
590 |
|
|
ENDIF |
591 |
|
|
|
592 |
|
|
END SUBROUTINE coefkz2 |
593 |
|
|
! |
594 |
|
|
!**************************************************************************************** |
595 |
|
|
! |
596 |
|
|
END MODULE coef_diff_turb_mod |