GCC Code Coverage Report

Directory: ./
File: phys/coef_diff_turb_mod.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 94 147 63.9%
Branches: 73 124 58.9%
Line Branch Exec Source
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 1920 SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, &
15 1920 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 3840 REAL, DIMENSION(klon,klev) :: ycoefm0, ycoefh0, yzlay, yteta
56 3840 REAL, DIMENSION(klon,klev+1) :: yzlev, q2diag, ykmm, ykmn, ykmq
57 1920 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
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 76339200 times.
✓ Branch 3 taken 76800 times.
 76417920 ykmm = 0 !ym missing init
68
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 76339200 times.
✓ Branch 3 taken 76800 times.
 76417920 ykmn = 0 !ym missing init
69
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 76339200 times.
✓ Branch 3 taken 76800 times.
 76417920 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 1920 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 1920 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 1920 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 1920 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
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
131
2/2
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
 74880 DO k=2,klev
132
2/2
✓ Branch 0 taken 29961176 times.
✓ Branch 1 taken 72960 times.
 30036056 DO i = 1, knon
133 yzlay(i,k)= &
134 yzlay(i,k-1)+RD*0.5*(yt(i,k-1)+yt(i,k)) &
135 30034136 /ypaprs(i,k)*(ypplay(i,k-1)-ypplay(i,k))/RG
136 END DO
137 END DO
138
139
2/2
✓ Branch 0 taken 74880 times.
✓ Branch 1 taken 1920 times.
 76800 DO k=1,klev
140
2/2
✓ Branch 0 taken 30749628 times.
✓ Branch 1 taken 74880 times.
 30826428 DO i = 1, knon
141 yteta(i,k)= &
142 yt(i,k)*(ypaprs(i,1)/ypplay(i,k))**RKAPPA &
143 74880 *(1.+0.61*yq(i,k))
144 END DO
145 END DO
146
147
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 yzlev(1:knon,1)=0.
148
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
149
2/2
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
 74880 DO k=2,klev
150
2/2
✓ Branch 0 taken 29961176 times.
✓ Branch 1 taken 72960 times.
 30036056 DO i = 1, knon
151 30034136 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 1920 CALL ustarhb(knon,yu,yv,ycdragm, yustar)
162
163
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 IF (prt_level > 9) THEN
164 ✗ WRITE(lunout,*) 'USTAR = ',yustar
165 ENDIF
166
167 ! iflag_pbl peut etre utilise comme longuer de melange
168 1920 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
1/2
✓ Branch 0 taken 1920 times.
✗ Branch 1 not taken.
 1920 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 1920 iflag_pbl,ydrgpro)
178 !FC
179 ENDIF
180
181
4/4
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 29961176 times.
✓ Branch 3 taken 72960 times.
 30036056 ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
182
4/4
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 29961176 times.
✓ Branch 3 taken 72960 times.
 30036056 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 1920 END SUBROUTINE coef_diff_turb
190 !
191 !****************************************************************************************
192 !
193 1920 SUBROUTINE coefkz(nsrf, knon, paprs, pplay, &
194 ksta, ksta_ter, &
195 ts, &
196 u,v,t,q, &
197 qsurf, &
198 1920 pcfm, pcfh)
199
200
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1920 times.
✓ Branch 4 taken 1920 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1920 times.
 30753468 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 3840 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 3840 REAL zgeop(klon,klev)
275 3840 REAL zmgeom(klon)
276 3840 REAL zri(klon)
277 3840 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 3840 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 1920 isommet=klev
298
299
1/2
✓ Branch 0 taken 1920 times.
✗ Branch 1 not taken.
 1920 IF (appel1er) THEN
300
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 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
2/2
✓ Branch 0 taken 74880 times.
✓ Branch 1 taken 1920 times.
 76800 DO k = 1, klev
313
2/2
✓ Branch 0 taken 30749628 times.
✓ Branch 1 taken 74880 times.
 30826428 DO i = 1, knon
314 30749628 pcfm(i,k) = 0.0
315 30824508 pcfh(i,k) = 0.0
316 ENDDO
317 ENDDO
318
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 DO i = 1, knon
319 790372 itop(i) = 0
320 ENDDO
321
322 !
323 ! Prescrire la valeur de contre-gradient
324 !
325 1920 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
2/2
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
 74880 DO k = 2, klev
332 74880 gamt(k) = 0.0
333 ENDDO
334 ENDIF
335 !IM cf JLD/ GKtest
336
2/2
✓ Branch 0 taken 1440 times.
✓ Branch 1 taken 480 times.
 1920 IF ( nsrf .NE. is_oce ) THEN
337 !IM 261103 kstable = kstable_ter
338 1440 kstable = ksta_ter
339 ELSE
340 !IM 261103 kstable = kstable_sinon
341 480 kstable = ksta
342 ENDIF
343 !IM cf JLD/ GKtest fin
344
345 !
346 ! Calculer les geopotentiels de chaque couche
347 !
348
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 DO i = 1, knon
349 zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) &
350 1920 * (paprs(i,1)-pplay(i,1))
351 ENDDO
352
2/2
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
 74880 DO k = 2, klev
353
2/2
✓ Branch 0 taken 29961176 times.
✓ Branch 1 taken 72960 times.
 30036056 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 72960 * (pplay(i,k-1)-pplay(i,k))
357 ENDDO
358 ENDDO
359
360 !
361 ! Calculer les coefficients turbulents dans l'atmosphere
362 !
363
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 DO i = 1, knon
364 790372 itop(i) = isommet
365 ENDDO
366
367
368
2/2
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
 74880 DO k = 2, isommet
369
2/2
✓ Branch 0 taken 29961176 times.
✓ Branch 1 taken 72960 times.
 30036056 DO i = 1, knon
370 zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 &
371 29961176 +(v(i,k)-v(i,k-1))**2)
372 29961176 zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
373 29961176 zdphi =zmgeom(i) / 2.0
374 29961176 zt = (t(i,k)+t(i,k-1)) * 0.5
375 29961176 zq = (q(i,k)+q(i,k-1)) * 0.5
376
377 !
378 ! Calculer Qs et dQs/dT:
379 !
380 IF (thermcep) THEN
381 29961176 zdelta = MAX(0.,SIGN(1.,RTT-zt))
382 zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta) &
383 29961176 + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
384 29961176 zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
385 29961176 zqs = MIN(0.5,zqs)
386 29961176 zcor = 1./(1.-RETV*zqs)
387 29961176 zqs = zqs*zcor
388 29961176 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 29961176 if (zq /= 0.) then
402 29961176 zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
403 else
404 zfr = 0.
405 end if
406 29961176 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 29961176 )*(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 29961176 )*(1.+RETV*q(i,k-1))
420 29961176 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 29961176 /(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 29961176 zcdn=SQRT(zdu2) / zmgeom(i) * RG
441
442 72960 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 29961176 /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
464 29961176 pcfm(i,k)=SQRT(MAX(zcdn*zcdn*(ric-zri(i))/ric, kstable))
465 29961176 pcfm(i,k)= zl2(i)* pcfm(i,k)
466 29961176 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
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 DO i = 1, knon
475
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 788452 times.
 790372 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 1920 END SUBROUTINE coefkz
484 !
485 !****************************************************************************************
486 !
487 ✗ SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, &
488 ✗ pcfm, pcfh)
489
490
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 29961176 times.
✗ Branch 3 not taken.
 90673900 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
597