GCC Code Coverage Report

Directory: ./
File: phys/yamada4.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 204 362 56.4%
Branches: 185 328 56.4%
Line Branch Exec Source
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2
3 3841 SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
4 1920 cd, tke, km, kn, kq, ustar, iflag_pbl, drgpro)
5
6 USE dimphy
7 USE ioipsl_getin_p_mod, ONLY : getin_p
8 USE phys_local_var_mod, only: tke_dissip
9
10 IMPLICIT NONE
11 include "iniprint.h"
12 ! .......................................................................
13 ! ym#include "dimensions.h"
14 ! ym#include "dimphy.h"
15 ! ************************************************************************************************
16 !
17 ! yamada4: subroutine qui calcule le transfert turbulent avec une fermeture d'ordre 2 ou 2.5
18 !
19 ! Reference: Simulation of nocturnal drainage flows by a q2l Turbulence Closure Model
20 ! Yamada T.
21 ! J. Atmos. Sci, 40, 91-106, 1983
22 !
23 !************************************************************************************************
24 ! Input :
25 !'======
26 ! ni: indice horizontal sur la grille de base, non restreinte
27 ! nsrf: type de surface
28 ! ngrid: nombre de mailles concern??es sur l'horizontal
29 ! dt : pas de temps
30 ! g : g
31 ! rconst: constante de l'air sec
32 ! zlev : altitude a chaque niveau (interface inferieure de la couche
33 ! de meme indice)
34 ! zlay : altitude au centre de chaque couche
35 ! u,v : vitesse au centre de chaque couche
36 ! (en entree : la valeur au debut du pas de temps)
37 ! teta : temperature potentielle virtuelle au centre de chaque couche
38 ! (en entree : la valeur au debut du pas de temps)
39 ! cd : cdrag pour la quantit?? de mouvement
40 ! (en entree : la valeur au debut du pas de temps)
41 ! ustar: vitesse de friction calcul??e par une formule diagnostique
42 ! iflag_pbl: flag pour choisir des options du sch??ma de turbulence
43 !
44 ! iflag_pbl doit valoir entre 6 et 9
45 ! l=6, on prend systematiquement une longueur d'equilibre
46 ! iflag_pbl=6 : MY 2.0
47 ! iflag_pbl=7 : MY 2.0.Fournier
48 ! iflag_pbl=8/9 : MY 2.5
49 ! iflag_pbl=8 with special obsolete treatments for convergence
50 ! with Cmpi5 NPv3.1 simulations
51 ! iflag_pbl=10/11 : New scheme M2 and N2 explicit and dissiptation exact
52 ! iflag_pbl=12 = 11 with vertical diffusion off q2
53 !
54 ! 2013/04/01 (FH hourdin@lmd.jussieu.fr)
55 ! Correction for very stable PBLs (iflag_pbl=10 and 11)
56 ! iflag_pbl=8 converges numerically with NPv3.1
57 ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l
58 ! -> the model can run with longer time-steps.
59 ! 2016/11/30 (EV etienne.vignon@lmd.ipsl.fr)
60 ! On met tke (=q2/2) en entr??e plut??t que q2
61 ! On corrige l'update de la tke
62 ! 2020/10/01 (EV)
63 ! On ajoute la dissipation de la TKE en diagnostique de sortie
64 !
65 ! Inpout/Output :
66 !==============
67 ! tke : tke au bas de chaque couche
68 ! (en entree : la valeur au debut du pas de temps)
69 ! (en sortie : la valeur a la fin du pas de temps)
70
71 ! Outputs:
72 !==========
73 ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
74 ! couche)
75 ! (en sortie : la valeur a la fin du pas de temps)
76 ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
77 ! (en sortie : la valeur a la fin du pas de temps)
78 !
79 !.......................................................................
80
81 !=======================================================================
82 ! Declarations:
83 !=======================================================================
84
85
86 ! Inputs/Outputs
87 !----------------
88 REAL dt, g, rconst
89 REAL plev(klon, klev+1), temp(klon, klev)
90 REAL ustar(klon)
91 3840 REAL kmin, qmin, pblhmin(klon), coriol(klon)
92 REAL zlev(klon, klev+1)
93 REAL zlay(klon, klev)
94 REAL u(klon, klev)
95 REAL v(klon, klev)
96 REAL teta(klon, klev)
97 REAL cd(klon)
98 REAL tke(klon, klev+1)
99 3840 REAL unsdz(klon, klev)
100 3840 REAL unsdzdec(klon, klev+1)
101 REAL kn(klon, klev+1)
102 REAL km(klon, klev+1)
103 INTEGER iflag_pbl, ngrid, nsrf
104 INTEGER ni(klon)
105
106 !FC
107 REAL drgpro(klon,klev)
108 3840 REAL winds(klon,klev)
109
110 ! Local
111 !-------
112
113 INCLUDE "clesphys.h"
114
115 3840 REAL q2(klon, klev+1)
116 3840 REAL kmpre(klon, klev+1), tmp2, qpre
117 3840 REAL mpre(klon, klev+1)
118 REAL kq(klon, klev+1)
119 3840 REAL ff(klon, klev+1), delta(klon, klev+1)
120 3840 REAL aa(klon, klev+1), aa0, aa1
121 INTEGER nlay, nlev
122
123 LOGICAL,SAVE :: hboville=.TRUE.
124 REAL,SAVE :: viscom,viscoh
125 !$OMP THREADPRIVATE( hboville,viscom,viscoh)
126 INTEGER ig, jg, k
127 REAL ri, zrif, zalpha, zsm, zsn
128 3840 REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
129 3840 REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
130 3840 REAL dtetadz(klon, klev+1)
131 REAL m2cstat, mcstat, kmcstat
132 3840 REAL l(klon, klev+1)
133 3840 REAL zz(klon, klev+1)
134 INTEGER iter
135 3840 REAL dissip(klon,klev), tkeprov,tkeexp, shear(klon,klev), buoy(klon,klev)
136 REAL :: disseff
137
138 REAL,SAVE :: ric0,ric,rifc, b1, kap
139 !$OMP THREADPRIVATE(ric0,ric,rifc,b1,kap)
140 DATA b1, kap/16.6, 0.4/
141 REAL,SAVE :: seuilsm, seuilalpha
142 !$OMP THREADPRIVATE(seuilsm, seuilalpha)
143 REAL,SAVE :: lmixmin
144 !$OMP THREADPRIVATE(lmixmin)
145 LOGICAL, SAVE :: new_yamada4
146 INTEGER, SAVE :: yamada4_num
147 !$OMP THREADPRIVATE(new_yamada4,yamada4_num)
148 REAL, SAVE :: yun,ydeux
149 !$OMP THREADPRIVATE(yun,ydeux)
150
151 REAL frif, falpha, fsm
152 REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &
153 lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev)
154 LOGICAL, SAVE :: firstcall = .TRUE.
155 !$OMP THREADPRIVATE(firstcall)
156
157 CHARACTER (len = 20) :: modname = 'yamada4'
158 CHARACTER (len = 80) :: abort_message
159
160
161
162 ! Fonctions utilis??es
163 !--------------------
164
165 frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
166 falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
167 fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
168
169
170
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1919 times.
 1920 IF (firstcall) THEN
171 ! Seuil dans le code de turbulence
172 1 new_yamada4=.false.
173 1 CALL getin_p('new_yamada4',new_yamada4)
174
175
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 IF (new_yamada4) THEN
176 ! Corrections et reglages issus du travail de these d'Etienne Vignon.
177 1 ric=0.143 ! qui donne des valeurs proches des seuils proposes
178 ! dans YAMADA 1983 : sm=0.0845 (0.085 dans Y83)
179 ! sm=1.1213 (1.12 dans Y83)
180 1 CALL getin_p('yamada4_ric',ric)
181 1 ric0=0.19489 ! ric=0.195 originalement, mais produisait sm<0
182 1 ric=min(ric,ric0) ! Au dela de ric0, sm devient n??gatif
183 1 rifc=frif(ric)
184 1 seuilsm=fsm(frif(ric))
185 1 seuilalpha=falpha(frif(ric))
186 1 yun=1.
187 1 ydeux=2.
188 1 hboville=.FALSE.
189 1 viscom=1.46E-5
190 1 viscoh=2.06E-5
191 !lmixmin=1.0E-3
192 1 lmixmin=0.
193 1 yamada4_num=5
194 ELSE
195 ✗ ric=0.195
196 ✗ rifc=0.191
197 ✗ seuilalpha=1.12
198 ✗ seuilsm=0.085
199 ✗ yun=2.
200 ✗ ydeux=1.
201 ✗ hboville=.TRUE.
202 ✗ viscom=0.
203 ✗ viscoh=0.
204 ✗ lmixmin=1.
205 ✗ yamada4_num=0
206 ENDIF
207
208 1 WRITE(lunout,*)'YAMADA4 RIc, RIfc, Sm_min, Alpha_min',ric,rifc,seuilsm,seuilalpha
209 1 firstcall = .FALSE.
210 1 CALL getin_p('lmixmin',lmixmin)
211 1 CALL getin_p('yamada4_hboville',hboville)
212 1 CALL getin_p('yamada4_num',yamada4_num)
213 END IF
214
215
216
217 !===============================================================================
218 ! Flags, tests et ??valuations de constantes
219 !===============================================================================
220
221 ! On utilise ou non la routine de Holstalg Boville pour les cas tres stables
222
223
224
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 IF (.NOT. (iflag_pbl>=6 .AND. iflag_pbl<=12)) THEN
225 ✗ abort_message='probleme de coherence dans appel a MY'
226 ✗ CALL abort_physic(modname,abort_message,1)
227 END IF
228
229
230 1920 nlay = klev
231 1920 nlev = klev + 1
232
233
234 !========================================================================
235 ! Calcul des increments verticaux
236 !=========================================================================
237
238
239 ! Attention: zlev n'est pas declare a nlev
240
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 DO ig = 1, ngrid
241 790372 zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1))
242 END DO
243
244
245
2/2
✓ Branch 0 taken 74880 times.
✓ Branch 1 taken 1920 times.
 76800 DO k = 1, nlay
246
2/2
✓ Branch 0 taken 30749628 times.
✓ Branch 1 taken 74880 times.
 30826428 DO ig = 1, ngrid
247 30824508 unsdz(ig, k) = 1.E+0/(zlev(ig,k+1)-zlev(ig,k))
248 END DO
249 END DO
250
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 DO ig = 1, ngrid
251 790372 unsdzdec(ig, 1) = 1.E+0/(zlay(ig,1)-zlev(ig,1))
252 END DO
253
2/2
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
 74880 DO k = 2, nlay
254
2/2
✓ Branch 0 taken 29961176 times.
✓ Branch 1 taken 72960 times.
 30036056 DO ig = 1, ngrid
255 30034136 unsdzdec(ig, k) = 1.E+0/(zlay(ig,k)-zlay(ig,k-1))
256 END DO
257 END DO
258
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 DO ig = 1, ngrid
259 790372 unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
260 END DO
261
262 !=========================================================================
263 ! Richardson number and stability functions
264 !=========================================================================
265
266 ! initialize arrays:
267
268
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 31538080 times.
✓ Branch 3 taken 76800 times.
 31616800 m2(1:ngrid, :) = 0.0
269
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 31538080 times.
✓ Branch 3 taken 76800 times.
 31616800 sm(1:ngrid, :) = 0.0
270
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 31538080 times.
✓ Branch 3 taken 76800 times.
 31616800 rif(1:ngrid, :) = 0.0
271
272 !------------------------------------------------------------
273
2/2
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
 74880 DO k = 2, klev
274
275
2/2
✓ Branch 0 taken 29961176 times.
✓ Branch 1 taken 72960 times.
 30036056 DO ig = 1, ngrid
276 29961176 dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
277 m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &
278 29961176 k-1))**2)/(dz(ig,k)*dz(ig,k))
279 29961176 dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k)
280 29961176 n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k))
281 29961176 ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
282
2/2
✓ Branch 0 taken 960908 times.
✓ Branch 1 taken 29000268 times.
 29961176 IF (ri<ric) THEN
283 960908 rif(ig, k) = frif(ri)
284 ELSE
285 29000268 rif(ig, k) = rifc
286 END IF
287
1/2
✓ Branch 0 taken 29961176 times.
✗ Branch 1 not taken.
 29961176 if (new_yamada4) then
288 29961176 alpha(ig, k) = max(falpha(rif(ig,k)),seuilalpha)
289 29961176 sm(ig, k) = max(fsm(rif(ig,k)),seuilsm)
290 else
291 ✗ IF (rif(ig,k)<0.16) THEN
292 ✗ alpha(ig, k) = falpha(rif(ig,k))
293 ✗ sm(ig, k) = fsm(rif(ig,k))
294 ELSE
295 ✗ alpha(ig, k) = seuilalpha
296 ✗ sm(ig, k) = seuilsm
297 END IF
298
299 end if
300 30034136 zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
301 END DO
302 END DO
303
304
305
306
307
308 !=======================================================================
309 ! DIFFERENT TYPES DE SCHEMA de YAMADA
310 !=======================================================================
311
312 ! On commence par calculer q2 a partir de la tke
313
314
1/2
✓ Branch 0 taken 1920 times.
✗ Branch 1 not taken.
 1920 IF (new_yamada4) THEN
315
2/2
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
 78720 DO k=1,klev+1
316
2/2
✓ Branch 0 taken 31538080 times.
✓ Branch 1 taken 76800 times.
 31616800 q2(1:ngrid,k)=tke(1:ngrid,k)*ydeux
317 ENDDO
318 ELSE
319 ✗ DO k=1,klev+1
320 ✗ q2(1:ngrid,k)=tke(1:ngrid,k)
321 ENDDO
322 ENDIF
323
324 ! ====================================================================
325 ! Computing the mixing length
326 ! ====================================================================
327
328
329 1920 CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l)
330
331
332 !--------------
333 ! Yamada 2.0
334 !--------------
335
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 IF (iflag_pbl==6) THEN
336
337 ✗ DO k = 2, klev
338 ✗ q2(1:ngrid, k) = l(1:ngrid, k)**2*zz(1:ngrid, k)
339 END DO
340
341
342 !------------------
343 ! Yamada 2.Fournier
344 !------------------
345
346
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 ELSE IF (iflag_pbl==7) THEN
347
348
349 ! Calcul de l, km, au pas precedent
350 !....................................
351 ✗ DO k = 2, klev
352 ✗ DO ig = 1, ngrid
353 ✗ delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
354 ✗ kmpre(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
355 ✗ mpre(ig, k) = sqrt(m2(ig,k))
356 END DO
357 END DO
358
359 ✗ DO k = 2, klev - 1
360 ✗ DO ig = 1, ngrid
361 ✗ m2cstat = max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1, 1.E-12)
362 ✗ mcstat = sqrt(m2cstat)
363
364 ! Puis on ecrit la valeur de q qui annule l'equation de m supposee en q3
365 !.........................................................................
366
367 ✗ IF (k==2) THEN
368 kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ &
369 unsdz(ig,k-1)*cd(ig)*(sqrt(u(ig,3)**2+v(ig,3)**2)-mcstat/unsdzdec &
370 (ig,k)-mpre(ig,k+1)/unsdzdec(ig,k+1))**2)/(unsdz(ig,k)+unsdz(ig,k &
371 ✗ -1))
372 ELSE
373 kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ &
374 unsdz(ig,k-1)*kmpre(ig,k-1)*mpre(ig,k-1))/ &
375 ✗ (unsdz(ig,k)+unsdz(ig,k-1))
376 END IF
377
378 ✗ tmp2 = kmcstat/(sm(ig,k)/q2(ig,k))/l(ig, k)
379 ✗ q2(ig, k) = max(tmp2, 1.E-12)**(2./3.)
380
381 END DO
382 END DO
383
384
385 ! ------------------------
386 ! Yamada 2.5 a la Didi
387 !-------------------------
388
389
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 ELSE IF (iflag_pbl==8 .OR. iflag_pbl==9) THEN
390
391 ! Calcul de l, km, au pas precedent
392 !....................................
393 ✗ DO k = 2, klev
394 ✗ DO ig = 1, ngrid
395 ✗ delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
396 ✗ IF (delta(ig,k)<1.E-20) THEN
397 ✗ delta(ig, k) = 1.E-20
398 END IF
399 ✗ km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
400 ✗ aa0 = (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
401 ✗ aa1 = (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
402 ✗ aa(ig, k) = aa1*dt/(delta(ig,k)*l(ig,k))
403 qpre = sqrt(q2(ig,k))
404 ✗ IF (aa(ig,k)>0.) THEN
405 ✗ q2(ig, k) = (qpre+aa(ig,k)*qpre*qpre)**2
406 ELSE
407 ✗ q2(ig, k) = (qpre/(1.-aa(ig,k)*qpre))**2
408 END IF
409 ! else ! iflag_pbl=9
410 ! if (aa(ig,k)*qpre.gt.0.9) then
411 ! q2(ig,k)=(qpre*10.)**2
412 ! else
413 ! q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
414 ! endif
415 ! endif
416 ✗ q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
417 END DO
418 END DO
419
420
1/2
✓ Branch 0 taken 1920 times.
✗ Branch 1 not taken.
 1920 ELSE IF (iflag_pbl>=10) THEN
421
422
4/4
✓ Branch 0 taken 74880 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 74430720 times.
✓ Branch 3 taken 74880 times.
 74507520 shear(:,:)=0.
423
4/4
✓ Branch 0 taken 74880 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 74430720 times.
✓ Branch 3 taken 74880 times.
 74507520 buoy(:,:)=0.
424
4/4
✓ Branch 0 taken 74880 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 74430720 times.
✓ Branch 3 taken 74880 times.
 74507520 dissip(:,:)=0.
425
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 76339200 times.
✓ Branch 3 taken 76800 times.
 76417920 km(:,:)=0.
426
427
1/2
✓ Branch 0 taken 1920 times.
✗ Branch 1 not taken.
 1920 IF (yamada4_num>=1) THEN
428
429
2/2
✓ Branch 0 taken 71040 times.
✓ Branch 1 taken 1920 times.
 72960 DO k = 2, klev - 1
430
2/2
✓ Branch 0 taken 29172724 times.
✓ Branch 1 taken 71040 times.
 29245684 DO ig=1,ngrid
431 29172724 q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
432 29172724 km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
433 29172724 shear(ig,k)=km(ig, k)*m2(ig, k)
434 29172724 buoy(ig,k)=km(ig, k)*m2(ig, k)*(-1.*rif(ig,k))
435 ! dissip(ig,k)=min(max(((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)),1.E-12),1.E4)
436 29243764 dissip(ig,k)=((sqrt(q2(ig,k)))**3)/(b1*l(ig,k))
437 ENDDO
438 ENDDO
439
440
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 IF (yamada4_num==1) THEN ! Schema du MAR tel quel
441 ✗ DO k = 2, klev - 1
442 ✗ DO ig=1,ngrid
443 ✗ tkeprov=q2(ig,k)/ydeux
444 tkeprov= tkeprov* &
445 & (tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))))/ &
446 ✗ & (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k)))
447 ✗ q2(ig,k)=tkeprov*ydeux
448 ENDDO
449 ENDDO
450
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 ELSE IF (yamada4_num==2) THEN ! version modifiee avec integration exacte pour la dissipation
451 ✗ DO k = 2, klev - 1
452 ✗ DO ig=1,ngrid
453 ✗ tkeprov=q2(ig,k)/ydeux
454 ✗ disseff=dissip(ig,k)-min(0.,buoy(ig,k))
455 ✗ tkeprov = tkeprov/(1.+dt*disseff/(2.*tkeprov))**2
456 ✗ tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
457 ✗ q2(ig,k)=tkeprov*ydeux
458 ! En cas stable, on traite la flotabilite comme la
459 ! dissipation, en supposant que buoy/q2^3 est constant.
460 ! Puis on prend la solution exacte
461 ENDDO
462 ENDDO
463
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 ELSE IF (yamada4_num==3) THEN ! version modifiee avec integration exacte pour la dissipation
464 ✗ DO k = 2, klev - 1
465 ✗ DO ig=1,ngrid
466 ✗ tkeprov=q2(ig,k)/ydeux
467 ✗ disseff=dissip(ig,k)-min(0.,buoy(ig,k))
468 ✗ tkeprov=tkeprov*exp(-dt*disseff/tkeprov)
469 ✗ tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
470 ✗ q2(ig,k)=tkeprov*ydeux
471 ! En cas stable, on traite la flotabilite comme la
472 ! dissipation, en supposant que buoy/q2^3 est constant.
473 ! Puis on prend la solution exacte
474 ENDDO
475 ENDDO
476
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 ELSE IF (yamada4_num==4) THEN ! version modifiee avec integration exacte pour la dissipation
477 ✗ DO k = 2, klev - 1
478 ✗ DO ig=1,ngrid
479 ✗ tkeprov=q2(ig,k)/ydeux
480 ✗ tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k)))
481 tkeprov= tkeprov* &
482 & tkeprov/ &
483 ✗ & (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k)))
484 ✗ q2(ig,k)=tkeprov*ydeux
485 ! En cas stable, on traite la flotabilite comme la
486 ! dissipation, en supposant que buoy/q2^3 est constant.
487 ! Puis on prend la solution exacte
488 ENDDO
489 ENDDO
490
491 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
492 !! Attention, yamada4_num=5 est inexacte car néglige les termes de flottabilité
493 !! en conditions instables
494 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
495
1/2
✓ Branch 0 taken 1920 times.
✗ Branch 1 not taken.
 1920 ELSE IF (yamada4_num==5) THEN ! version modifiee avec integration exacte pour la dissipation
496
2/2
✓ Branch 0 taken 1920 times.
✓ Branch 1 taken 71040 times.
 72960 DO k = 2, klev - 1
497
2/2
✓ Branch 0 taken 29172724 times.
✓ Branch 1 taken 71040 times.
 29245684 DO ig=1,ngrid
498 29172724 tkeprov=q2(ig,k)/ydeux
499
500 ! if(ifl_pbltree .eq. 0) then
501 ! disseff=dissip(ig,k)-min(0.,buoy(ig,k))
502 ! tkeexp=exp(-dt*disseff/tkeprov)
503 ! tkeprov= shear(ig,k)*tkeprov/disseff*(1.-tkeexp)+tkeprov*tkeexp
504 ! else
505 !FC on ajoute la dissipation due aux arbres
506 29172724 disseff=dissip(ig,k)-min(0.,buoy(ig,k)) + drgpro(ig,k)*tkeprov
507 29172724 tkeexp=exp(-dt*disseff/tkeprov)
508 ! on prend en compte la tke cree par les arbres
509 29172724 winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2)
510 tkeprov= (shear(ig,k)+ &
511 29172724 & drgpro(ig,k)*(winds(ig,k))**3)*tkeprov/disseff*(1.-tkeexp)+tkeprov*tkeexp
512 ! endif
513
514 29243764 q2(ig,k)=tkeprov*ydeux
515
516 ! En cas stable, on traite la flotabilite comme la
517 ! dissipation, en supposant que buoy/q2^3 est constant.
518 ! Puis on prend la solution exacte
519 ENDDO
520 ENDDO
521 ✗ ELSE IF (yamada4_num==6) THEN ! version modifiee avec integration exacte pour la dissipation
522 ✗ DO k = 2, klev - 1
523 ✗ DO ig=1,ngrid
524 !tkeprov=q2(ig,k)/ydeux
525 !tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k),0.)*dt
526 !disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k))
527 !tkeexp=exp(-dt*disseff/tkeprov)
528 !tkeprov= tkeprov*tkeexp
529 !q2(ig,k)=tkeprov*ydeux
530 ! En cas stable, on traite la flotabilite comme la
531 ! dissipation, en supposant que dissipeff/TKE est constant.
532 ! Puis on prend la solution exacte
533 !
534 ! With drag and dissipation from high vegetation (EV & FC, 05/10/2020)
535 ✗ winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2)
536 ✗ tkeprov=q2(ig,k)/ydeux
537 ✗ tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3,0.)*dt
538 ✗ disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3) + drgpro(ig,k)*tkeprov
539 ✗ tkeexp=exp(-dt*disseff/tkeprov)
540 ✗ tkeprov= tkeprov*tkeexp
541 ✗ q2(ig,k)=tkeprov*ydeux
542
543 ENDDO
544 ENDDO
545 ENDIF
546
547
2/2
✓ Branch 0 taken 1920 times.
✓ Branch 1 taken 71040 times.
 72960 DO k = 2, klev - 1
548
2/2
✓ Branch 0 taken 29172724 times.
✓ Branch 1 taken 71040 times.
 29245684 DO ig=1,ngrid
549 29243764 q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
550 ENDDO
551 ENDDO
552
553 ELSE
554
555 ✗ DO k = 2, klev - 1
556 ✗ km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k)
557 ✗ q2(1:ngrid, k) = q2(1:ngrid, k) + ydeux*dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
558 ! q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
559 ✗ q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4)
560 ✗ q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1))
561 ! q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1))
562 ✗ q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k)
563 END DO
564
565 ENDIF
566
567 ELSE
568 ✗ abort_message='Cas nom prevu dans yamada4'
569 ✗ CALL abort_physic(modname,abort_message,1)
570
571 END IF ! Fin du cas 8
572
573
574 ! ====================================================================
575 ! Calcul des coefficients de melange
576 ! ====================================================================
577
578
2/2
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
 74880 DO k = 2, klev
579
2/2
✓ Branch 0 taken 29961176 times.
✓ Branch 1 taken 72960 times.
 30036056 DO ig = 1, ngrid
580 29961176 zq = sqrt(q2(ig,k))
581 29961176 km(ig, k) = l(ig, k)*zq*sm(ig, k) ! For momentum
582 29961176 kn(ig, k) = km(ig, k)*alpha(ig, k) ! For scalars
583 30034136 kq(ig, k) = l(ig, k)*zq*0.2 ! For TKE
584 END DO
585 END DO
586
587
588 !====================================================================
589 ! Transport diffusif vertical de la TKE par la TKE
590 !====================================================================
591
592
593 ! initialize near-surface and top-layer mixing coefficients
594 !...........................................................
595
596
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface
597
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 kq(1:ngrid, klev+1) = 0 ! zero at the top
598
599 ! Transport diffusif vertical de la TKE.
600 !.......................................
601
602
1/2
✓ Branch 0 taken 1920 times.
✗ Branch 1 not taken.
 1920 IF (iflag_pbl>=12) THEN
603
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 q2(1:ngrid, 1) = q2(1:ngrid, 2)
604 1920 CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2)
605 END IF
606
607
608 !====================================================================
609 ! Traitement particulier pour les cas tres stables, introduction d'une
610 ! longueur de m??lange minimale
611 !====================================================================
612 !
613 ! Reference: Local versus Nonlocal boundary-layer diffusion in a global climate model
614 ! Holtslag A.A.M. and Boville B.A.
615 ! J. Clim., 6, 1825-1842, 1993
616
617
618
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 IF (hboville) THEN
619
620
621 ✗ IF (prt_level>1) THEN
622 ✗ WRITE(lunout,*) 'YAMADA4 0'
623 END IF
624
625 ✗ DO ig = 1, ngrid
626 ✗ coriol(ig) = 1.E-4
627 ✗ pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5)
628 END DO
629
630 IF (1==1) THEN
631 ✗ IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
632
633 ✗ DO k = 2, klev
634 ✗ DO ig = 1, ngrid
635 ✗ IF (teta(ig,2)>teta(ig,1)) THEN
636 ✗ qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
637 ✗ kmin = kap*zlev(ig, k)*qmin
638 ELSE
639 kmin = -1. ! kmin n'est utilise que pour les SL stables.
640 END IF
641 ✗ IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
642
643 ✗ kn(ig, k) = kmin
644 ✗ km(ig, k) = kmin
645 ✗ kq(ig, k) = kmin
646
647 ! la longueur de melange est suposee etre l= kap z
648 ! K=l q Sm d'ou q2=(K/l Sm)**2
649
650 ✗ q2(ig, k) = (qmin/sm(ig,k))**2
651 END IF
652 END DO
653 END DO
654
655 ELSE
656 ✗ DO k = 2, klev
657 ✗ DO ig = 1, ngrid
658 ✗ IF (teta(ig,2)>teta(ig,1)) THEN
659 ✗ qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
660 ✗ kmin = kap*zlev(ig, k)*qmin
661 ELSE
662 kmin = -1. ! kmin n'est utilise que pour les SL stables.
663 END IF
664 ✗ IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
665 kn(ig, k) = kmin
666 km(ig, k) = kmin
667 kq(ig, k) = kmin
668 ! la longueur de melange est suposee etre l= kap z
669 ! K=l q Sm d'ou q2=(K/l Sm)**2
670 ✗ sm(ig, k) = 1.
671 ✗ alpha(ig, k) = 1.
672 ✗ q2(ig, k) = min((qmin/sm(ig,k))**2, 10.)
673 ✗ zq = sqrt(q2(ig,k))
674 ✗ km(ig, k) = l(ig, k)*zq*sm(ig, k)
675 ✗ kn(ig, k) = km(ig, k)*alpha(ig, k)
676 ✗ kq(ig, k) = l(ig, k)*zq*0.2
677 END IF
678 END DO
679 END DO
680 END IF
681
682 END IF
683
684 END IF ! hboville
685
686 ! Ajout d'une viscosite moleculaire
687
4/4
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 29961176 times.
✓ Branch 3 taken 72960 times.
 30036056 km(1:ngrid,2:klev)=km(1:ngrid,2:klev)+viscom
688
4/4
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 29961176 times.
✓ Branch 3 taken 72960 times.
 30036056 kn(1:ngrid,2:klev)=kn(1:ngrid,2:klev)+viscoh
689
4/4
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 29961176 times.
✓ Branch 3 taken 72960 times.
 30036056 kq(1:ngrid,2:klev)=kq(1:ngrid,2:klev)+viscoh
690
691
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 IF (prt_level>1) THEN
692 ✗ WRITE(lunout,*)'YAMADA4 1'
693 END IF !(prt_level>1) THEN
694
695
696 !======================================================
697 ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane
698 !======================================================
699 !
700 ! Reference: A New Second-Order Turbulence Closure Scheme for the Planetary Boundary Layer
701 ! Abdella K and McFarlane N
702 ! J. Atmos. Sci., 54, 1850-1867, 1997
703
704 ! Diagnostique pour stokage
705 !..........................
706
707 IF (1==0) THEN
708 rino = rif
709 smyam(1:ngrid, 1) = 0.
710 styam(1:ngrid, 1) = 0.
711 lyam(1:ngrid, 1) = 0.
712 knyam(1:ngrid, 1) = 0.
713 w2yam(1:ngrid, 1) = 0.
714 t2yam(1:ngrid, 1) = 0.
715
716 smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)
717 styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev)
718 lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev)
719 knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)
720
721
722 ! Calcul de w'2 et T'2
723 !.......................
724
725 w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &
726 lyam(1:ngrid, 2:klev)*5.17*kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ &
727 sqrt(q2(1:ngrid,2:klev))
728
729 t2yam(1:ngrid, 2:klev) = 9.1*kn(1:ngrid, 2:klev)* &
730 dtetadz(1:ngrid, 2:klev)**2/sqrt(q2(1:ngrid,2:klev))* &
731 lyam(1:ngrid, 2:klev)
732 END IF
733
734
735
736 !============================================================================
737 ! Mise a jour de la tke
738 !============================================================================
739
740
1/2
✓ Branch 0 taken 1920 times.
✗ Branch 1 not taken.
 1920 IF (new_yamada4) THEN
741
2/2
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
 78720 DO k=1,klev+1
742
2/2
✓ Branch 0 taken 31538080 times.
✓ Branch 1 taken 76800 times.
 31616800 tke(1:ngrid,k)=q2(1:ngrid,k)/ydeux
743 ENDDO
744 ELSE
745 ✗ DO k=1,klev+1
746 ✗ tke(1:ngrid,k)=q2(1:ngrid,k)
747 ENDDO
748 ENDIF
749
750
751 !============================================================================
752 ! Diagnostique de la dissipation
753 !============================================================================
754
755 ! Diagnostics
756
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 31538080 times.
✓ Branch 3 taken 76800 times.
 31616800 tke_dissip(1:ngrid,:,nsrf)=0.
757 ! DO k=2,klev
758 ! DO ig=1,ngrid
759 ! jg=ni(ig)
760 ! tke_dissip(jg,k,nsrf)=dissip(ig,k)
761 ! ENDDO
762 ! ENDDO
763
764 !=============================================================================
765
766 1920 RETURN
767
768
769 END SUBROUTINE yamada4
770
771 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
772
773 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
774 1920 SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
775
776 USE dimphy
777 IMPLICIT NONE
778
779 ! vdif_q2: subroutine qui calcule la diffusion de la TKE par la TKE
780 ! avec un schema implicite en temps avec
781 ! inversion d'un syst??me tridiagonal
782 !
783 ! Reference: Description of the interface with the surface and
784 ! the computation of the turbulet diffusion in LMDZ
785 ! Technical note on LMDZ
786 ! Dufresne, J-L, Ghattas, J. and Grandpeix, J-Y
787 !
788 !============================================================================
789 ! Declarations
790 !============================================================================
791
792 REAL plev(klon, klev+1)
793 REAL temp(klon, klev)
794 REAL timestep
795 REAL gravity, rconst
796 3840 REAL kstar(klon, klev+1), zz
797 REAL kmy(klon, klev+1)
798 REAL q2(klon, klev+1)
799 3840 REAL deltap(klon, klev+1)
800 1920 REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
801 INTEGER ngrid
802
803 INTEGER i, k
804
805
806 !=========================================================================
807 ! Calcul
808 !=========================================================================
809
810
2/2
✓ Branch 0 taken 1920 times.
✓ Branch 1 taken 74880 times.
 76800 DO k = 1, klev
811
2/2
✓ Branch 0 taken 30749628 times.
✓ Branch 1 taken 74880 times.
 30826428 DO i = 1, ngrid
812 30749628 zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
813 kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
814 30824508 (plev(i,k)-plev(i,k+1))*timestep
815 END DO
816 END DO
817
818
2/2
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
 74880 DO k = 2, klev
819
2/2
✓ Branch 0 taken 29961176 times.
✓ Branch 1 taken 72960 times.
 30036056 DO i = 1, ngrid
820 30034136 deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
821 END DO
822 END DO
823
2/2
✓ Branch 0 taken 1920 times.
✓ Branch 1 taken 788452 times.
 790372 DO i = 1, ngrid
824 788452 deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
825 788452 deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
826 788452 denom(i, klev+1) = deltap(i, klev+1) + kstar(i, klev)
827 788452 alpha(i, klev+1) = deltap(i, klev+1)*q2(i, klev+1)/denom(i, klev+1)
828 790372 beta(i, klev+1) = kstar(i, klev)/denom(i, klev+1)
829 END DO
830
831
2/2
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
 74880 DO k = klev, 2, -1
832
2/2
✓ Branch 0 taken 29961176 times.
✓ Branch 1 taken 72960 times.
 30036056 DO i = 1, ngrid
833 denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + &
834 29961176 kstar(i, k-1)
835 29961176 alpha(i, k) = (q2(i,k)*deltap(i,k)+kstar(i,k)*alpha(i,k+1))/denom(i, k)
836 30034136 beta(i, k) = kstar(i, k-1)/denom(i, k)
837 END DO
838 END DO
839
840 ! Si on recalcule q2(1)
841 !.......................
842 IF (1==0) THEN
843 DO i = 1, ngrid
844 denom(i, 1) = deltap(i, 1) + (1-beta(i,2))*kstar(i, 1)
845 q2(i, 1) = (q2(i,1)*deltap(i,1)+kstar(i,1)*alpha(i,2))/denom(i, 1)
846 END DO
847 END IF
848
849
850
2/2
✓ Branch 0 taken 74880 times.
✓ Branch 1 taken 1920 times.
 76800 DO k = 2, klev + 1
851
2/2
✓ Branch 0 taken 30749628 times.
✓ Branch 1 taken 74880 times.
 30826428 DO i = 1, ngrid
852 30824508 q2(i, k) = alpha(i, k) + beta(i, k)*q2(i, k-1)
853 END DO
854 END DO
855
856 1920 RETURN
857 END SUBROUTINE vdif_q2
858 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
859
860
861
862 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
863 ✗ SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
864
865 USE dimphy
866 IMPLICIT NONE
867
868 ! vdif_q2e: subroutine qui calcule la diffusion de TKE par la TKE
869 ! avec un schema explicite en temps
870
871
872 !====================================================
873 ! Declarations
874 !====================================================
875
876 REAL plev(klon, klev+1)
877 REAL temp(klon, klev)
878 REAL timestep
879 REAL gravity, rconst
880 ✗ REAL kstar(klon, klev+1), zz
881 REAL kmy(klon, klev+1)
882 REAL q2(klon, klev+1)
883 ✗ REAL deltap(klon, klev+1)
884 REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
885 INTEGER ngrid
886 INTEGER i, k
887
888
889 !==================================================
890 ! Calcul
891 !==================================================
892
893 ✗ DO k = 1, klev
894 ✗ DO i = 1, ngrid
895 ✗ zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
896 kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
897 ✗ (plev(i,k)-plev(i,k+1))*timestep
898 END DO
899 END DO
900
901 ✗ DO k = 2, klev
902 ✗ DO i = 1, ngrid
903 ✗ deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
904 END DO
905 END DO
906 ✗ DO i = 1, ngrid
907 ✗ deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
908 ✗ deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
909 END DO
910
911 ✗ DO k = klev, 2, -1
912 ✗ DO i = 1, ngrid
913 q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, &
914 ✗ k))-kstar(i,k-1)*(q2(i,k)-q2(i,k-1)))/deltap(i, k)
915 END DO
916 END DO
917
918 ✗ DO i = 1, ngrid
919 ✗ q2(i, 1) = q2(i, 1) + (kstar(i,1)*(q2(i,2)-q2(i,1)))/deltap(i, 1)
920 q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, &
921 ✗ klev)))/deltap(i, klev+1)
922 END DO
923
924 ✗ RETURN
925 END SUBROUTINE vdif_q2e
926
927 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
928
929
930 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
931
932 1920 SUBROUTINE mixinglength(ni, nsrf, ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, lmix)
933
934
935
936 USE dimphy
937 USE phys_state_var_mod, only: zstd, zsig, zmea
938 USE phys_local_var_mod, only: l_mixmin, l_mix
939
940 ! zstd: ecart type de la'altitud e sous-maille
941 ! zmea: altitude moyenne sous maille
942 ! zsig: pente moyenne de le maille
943
944 USE geometry_mod, only: cell_area
945 ! aire_cell: aire de la maille
946
947 IMPLICIT NONE
948 !*************************************************************************
949 ! Subrourine qui calcule la longueur de m??lange dans le sch??ma de turbulence
950 ! avec la formule de Blackadar
951 ! Calcul d'un minimum en fonction de l'orographie sous-maille:
952 ! L'id??e est la suivante: plus il y a de relief, plus il y a du m??lange
953 ! induit par les circulations meso et submeso ??chelles.
954 !
955 ! References: * The vertical distribution of wind and turbulent exchange in a neutral atmosphere
956 ! Blackadar A.K.
957 ! J. Geophys. Res., 64, No 8, 1962
958 !
959 ! * An evaluation of neutral and convective planetary boundary-layer parametrisations relative
960 ! to large eddy simulations
961 ! Ayotte K et al
962 ! Boundary Layer Meteorology, 79, 131-175, 1996
963 !
964 !
965 ! * Local Similarity in the Stable Boundary Layer and Mixing length Approaches: consistency of concepts
966 ! Van de Wiel B.J.H et al
967 ! Boundary-Lay Meteorol, 128, 103-166, 2008
968 !
969 !
970 ! Histoire:
971 !----------
972 ! * premi??re r??daction, Etienne et Frederic, 09/06/2016
973 !
974 ! ***********************************************************************
975
976 !==================================================================
977 ! Declarations
978 !==================================================================
979
980 ! Inputs
981 !-------
982 INTEGER ni(klon) ! indice sur la grille original (non restreinte)
983 INTEGER nsrf ! Type de surface
984 INTEGER ngrid ! Nombre de points concern??s sur l'horizontal
985 INTEGER iflag_pbl ! Choix du sch??ma de turbulence
986 REAL pbl_lmixmin_alpha ! on active ou non le calcul de la longueur de melange minimum
987 REAL lmixmin ! Minimum absolu de la longueur de m??lange
988 REAL zlay(klon, klev) ! altitude du centre de la couche
989 REAL zlev(klon, klev+1) ! atitude de l'interface inf??rieure de la couche
990 REAL u(klon, klev) ! vitesse du vent zonal
991 REAL v(klon, klev) ! vitesse du vent meridional
992 REAL q2(klon, klev+1) ! energie cin??tique turbulente
993 REAL n2(klon, klev+1) ! frequence de Brunt-Vaisala
994
995 !In/out
996 !-------
997
998 LOGICAL, SAVE :: firstcall = .TRUE.
999 !$OMP THREADPRIVATE(firstcall)
1000
1001 ! Outputs
1002 !---------
1003
1004 REAL lmix(klon, klev+1) ! Longueur de melange
1005
1006
1007 ! Local
1008 !-------
1009
1010 INTEGER ig,jg, k
1011 3840 REAL h_oro(klon)
1012 3840 REAL hlim(klon)
1013 REAL, SAVE :: kap=0.4,kapb=0.4
1014 !$OMP THREADPRIVATE(kap,kapb)
1015 REAL zq
1016 3840 REAL sq(klon), sqz(klon)
1017 REAL, ALLOCATABLE, SAVE :: l0(:)
1018 !$OMP THREADPRIVATE(l0)
1019 REAL fl, zzz, zl0, zq2, zn2
1020 REAL famorti, zzzz, zh_oro, zhlim
1021 3840 REAL l1(klon, klev+1), l2(klon,klev+1)
1022 3840 REAL winds(klon, klev)
1023 REAL xcell
1024 REAL zstdslope(klon)
1025 REAL lmax
1026 REAL l2strat, l2neutre, extent
1027 ✗ REAL l2limit(klon)
1028 !===============================================================
1029 ! Fonctions utiles
1030 !===============================================================
1031
1032 ! Calcul de l suivant la formule de Blackadar 1962 adapt??e par Ayotte 1996
1033 !..........................................................................
1034
1035 fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, &
1036 k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( &
1037 max(n2(ig,k),1.E-10))), 1.E-5)
1038
1039 ! Fonction d'amortissement de la turbulence au dessus de la montagne
1040 ! On augmente l'amortissement en diminuant la valeur de hlim (extent) dans le code
1041 !.....................................................................
1042
1043 famorti(zzzz, zh_oro, zhlim)=(-1.)*ATAN((zzzz-zh_oro)/(zhlim-zh_oro))*2./3.1416+1.
1044
1045
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 IF (ngrid==0) RETURN
1046
1047
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1919 times.
 1920 IF (firstcall) THEN
1048
3/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
 1 ALLOCATE (l0(klon))
1049 1 firstcall = .FALSE.
1050 END IF
1051
1052
1053 !=====================================================================
1054 ! CALCUL de la LONGUEUR de m??lange suivant BLACKADAR: l1
1055 !=====================================================================
1056
1057
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 31538080 times.
✓ Branch 3 taken 76800 times.
 31616800 l1(1:ngrid,:)=0.
1058
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
 1920 IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN
1059
1060
1061 ! Iterative computation of l0
1062 ! This version is kept for iflag_pbl only for convergence
1063 ! with NPv3.1 Cmip5 simulations
1064 !...................................................................
1065
1066 ✗ DO ig = 1, ngrid
1067 ✗ sq(ig) = 1.E-10
1068 ✗ sqz(ig) = 1.E-10
1069 END DO
1070 ✗ DO k = 2, klev - 1
1071 ✗ DO ig = 1, ngrid
1072 ✗ zq = sqrt(q2(ig,k))
1073 ✗ sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
1074 ✗ sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
1075 END DO
1076 END DO
1077 ✗ DO ig = 1, ngrid
1078 ✗ l0(ig) = 0.2*sqz(ig)/sq(ig)
1079 END DO
1080 ✗ DO k = 2, klev
1081 ✗ DO ig = 1, ngrid
1082 ✗ l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
1083 END DO
1084 END DO
1085
1086 ELSE
1087
1088
1089 ! In all other case, the assymptotic mixing length l0 is imposed (150m)
1090 !......................................................................
1091
1092
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 l0(1:ngrid) = 150.
1093
2/2
✓ Branch 0 taken 1920 times.
✓ Branch 1 taken 72960 times.
 74880 DO k = 2, klev
1094
2/2
✓ Branch 0 taken 29961176 times.
✓ Branch 1 taken 72960 times.
 30036056 DO ig = 1, ngrid
1095 30034136 l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
1096 END DO
1097 END DO
1098
1099 END IF
1100
1101 !=================================================================================
1102 ! CALCUL d'une longueur de melange en fonctions de la topographie sous maille: l2
1103 ! si plb_lmixmin_alpha=TRUE et si on se trouve sur de la terre ( pas actif sur les
1104 ! glacier, la glace de mer et les oc??ans)
1105 !=================================================================================
1106
1107
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 31538080 times.
✓ Branch 3 taken 76800 times.
 31616800 l2(1:ngrid,:)=0.0
1108
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 31538080 times.
✓ Branch 3 taken 76800 times.
 31616800 l_mixmin(1:ngrid,:,nsrf)=0.
1109
4/4
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
✓ Branch 2 taken 31538080 times.
✓ Branch 3 taken 76800 times.
 31616800 l_mix(1:ngrid,:,nsrf)=0.
1110
1111
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 1440 times.
 1920 IF (nsrf .EQ. 1) THEN
1112
1113 ! coefficients
1114 !--------------
1115
1116 extent=2. ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1.
1117 lmax=150. ! Longueur de m??lange max dans l'absolu
1118
1119 ! calculs
1120 !---------
1121
1122
2/2
✓ Branch 0 taken 247680 times.
✓ Branch 1 taken 480 times.
 248160 DO ig=1,ngrid
1123
1124 ! On calcule la hauteur du relief
1125 !.................................
1126 ! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille
1127 ! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon)
1128 ! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief
1129 ! (en gros, une maille de taille xcell avec une pente constante zstdslope)
1130 247680 jg=ni(ig)
1131 ! IF (zsig(jg) .EQ. 0.) THEN
1132 ! zstdslope(ig)=0.
1133 ! ELSE
1134 ! xcell=sqrt(cell_area(jg))
1135 ! zstdslope(ig)=max((xcell*zsig(jg)-zmea(jg))**3 /(3.*zsig(jg)),0.)
1136 ! zstdslope(ig)=sqrt(zstdslope(ig))
1137 ! END IF
1138
1139 ! h_oro(ig)=max(zstd(jg)-zstdslope(ig),0.) ! Hauteur du relief
1140 247680 h_oro(ig)=zstd(jg)
1141 248160 hlim(ig)=extent*h_oro(ig)
1142 ENDDO
1143
1144
2/2
✓ Branch 0 taken 247680 times.
✓ Branch 1 taken 480 times.
 248160 l2limit(1:ngrid)=0.
1145
1146
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 18240 times.
 18720 DO k=2,klev
1147
2/2
✓ Branch 0 taken 9411840 times.
✓ Branch 1 taken 18240 times.
 9430560 DO ig=1,ngrid
1148 9411840 winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2)
1149
2/2
✓ Branch 0 taken 453017 times.
✓ Branch 1 taken 8958823 times.
 9430080 IF (zlev(ig,k) .LE. h_oro(ig)) THEN ! sous l'orographie
1150 453017 l2strat= kapb*pbl_lmixmin_alpha*winds(ig,k)/sqrt(max(n2(ig,k),1.E-10)) ! si stratifi??, amplitude d'oscillation * kappab (voir Van de Wiel et al 2008)
1151 453017 l2neutre=kap*zlev(ig,k)*h_oro(ig)/(kap*zlev(ig,k)+h_oro(ig)) ! Dans le cas neutre, formule de blackadar. tend asymptotiquement vers h
1152 453017 l2neutre=MIN(l2neutre,lmax) ! On majore par lmax
1153 453017 l2limit(ig)=MIN(l2neutre,l2strat) ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e)
1154 453017 l2(ig,k)=l2limit(ig)
1155
1156
2/2
✓ Branch 0 taken 255452 times.
✓ Branch 1 taken 8703371 times.
 8958823 ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles
1157
1158 ! Au dessus des montagnes, on prend la l2limit au sommet des montagnes
1159 ! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z)
1160 ! et on multiplie l2limit par une fonction qui d??croit entre h et hlim
1161 255452 l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig))
1162 ELSE ! Au dessus de extent*h, on prend l2=l0
1163 8703371 l2(ig,k)=0.
1164 END IF
1165 ENDDO
1166 ENDDO
1167 ENDIF ! pbl_lmixmin_alpha
1168
1169 !==================================================================================
1170 ! On prend le max entre la longueur de melange de blackadar et celle calcul??e
1171 ! en fonction de la topographie
1172 !===================================================================================
1173
1174
1175
2/2
✓ Branch 0 taken 76800 times.
✓ Branch 1 taken 1920 times.
 78720 DO k=1,klev+1
1176
2/2
✓ Branch 0 taken 31538080 times.
✓ Branch 1 taken 76800 times.
 31616800 DO ig=1,ngrid
1177 31614880 lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin)
1178 ENDDO
1179 ENDDO
1180
1181 ! Diagnostics
1182
1183
2/2
✓ Branch 0 taken 72960 times.
✓ Branch 1 taken 1920 times.
 74880 DO k=2,klev
1184
2/2
✓ Branch 0 taken 29961176 times.
✓ Branch 1 taken 72960 times.
 30036056 DO ig=1,ngrid
1185 29961176 jg=ni(ig)
1186 29961176 l_mix(jg,k,nsrf)=lmix(ig,k)
1187 30034136 l_mixmin(jg,k,nsrf)=l2(ig,k)
1188 ENDDO
1189 ENDDO
1190
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
 790372 DO ig=1,ngrid
1191 788452 jg=ni(ig)
1192 790372 l_mix(jg,1,nsrf)=hlim(ig)
1193 ENDDO
1194
1195
1196
1197 1920 END SUBROUTINE mixinglength
1198