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 |