GCC Code Coverage Report


Directory: ./
File: phys/thermcell_old.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 1771 0.0%
Branches: 0 1518 0.0%

Line Branch Exec Source
1 ✗ SUBROUTINE thermcell_2002(ngrid, nlay, ptimestep, iflag_thermals, pplay, &
2 ✗ pplev, pphi, pu, pv, pt, po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0, &
3 ✗ fraca, wa_moy, r_aspect, l_mix, w2di, tho)
4
5 USE dimphy
6 USE write_field_phy
7 IMPLICIT NONE
8
9 ! =======================================================================
10
11 ! Calcul du transport verticale dans la couche limite en presence
12 ! de "thermiques" explicitement representes
13
14 ! R��criture � partir d'un listing papier � Habas, le 14/02/00
15
16 ! le thermique est suppos� homog�ne et dissip� par m�lange avec
17 ! son environnement. la longueur l_mix contr�le l'efficacit� du
18 ! m�lange
19
20 ! Le calcul du transport des diff�rentes esp�ces se fait en prenant
21 ! en compte:
22 ! 1. un flux de masse montant
23 ! 2. un flux de masse descendant
24 ! 3. un entrainement
25 ! 4. un detrainement
26
27 ! =======================================================================
28
29 ! -----------------------------------------------------------------------
30 ! declarations:
31 ! -------------
32
33 include "YOMCST.h"
34
35 ! arguments:
36 ! ----------
37
38 INTEGER ngrid, nlay, w2di, iflag_thermals
39 REAL tho
40 REAL ptimestep, l_mix, r_aspect
41 REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
42 REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
43 REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
44 REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
45 REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
46 REAL pphi(ngrid, nlay)
47 ✗ REAL fraca(ngrid, nlay+1), zw2(ngrid, nlay+1)
48
49 INTEGER, SAVE :: idetr = 3, lev_out = 1
50 !$OMP THREADPRIVATE(idetr,lev_out)
51
52 ! local:
53 ! ------
54
55 INTEGER, SAVE :: dvdq = 0, flagdq = 0, dqimpl = 1
56 LOGICAL, SAVE :: debut = .TRUE.
57 !$OMP THREADPRIVATE(dvdq,flagdq,debut,dqimpl)
58
59 ✗ INTEGER ig, k, l, lmax(klon, klev+1), lmaxa(klon), lmix(klon)
60 ✗ REAL zmax(klon), zw, zz, ztva(klon, klev), zzz
61
62 ✗ REAL zlev(klon, klev+1), zlay(klon, klev)
63 ✗ REAL zh(klon, klev), zdhadj(klon, klev)
64 ✗ REAL ztv(klon, klev)
65 ✗ REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
66 REAL wh(klon, klev+1)
67 REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
68 REAL zla(klon, klev+1)
69 REAL zwa(klon, klev+1)
70 REAL zld(klon, klev+1)
71 REAL zwd(klon, klev+1)
72 REAL zsortie(klon, klev)
73 ✗ REAL zva(klon, klev)
74 ✗ REAL zua(klon, klev)
75 ✗ REAL zoa(klon, klev)
76
77 ✗ REAL zha(klon, klev)
78 REAL wa_moy(klon, klev+1)
79 ✗ REAL fracc(klon, klev+1)
80 REAL zf, zf2
81 ✗ REAL thetath2(klon, klev), wth2(klon, klev)
82 ! common/comtherm/thetath2,wth2
83
84 REAL count_time
85
86 LOGICAL sorties
87 ✗ REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
88 ✗ REAL zpspsk(klon, klev)
89
90 ✗ REAL wmax(klon, klev), wmaxa(klon)
91
92 ✗ REAL wa(klon, klev, klev+1)
93 ✗ REAL wd(klon, klev+1)
94 REAL larg_part(klon, klev, klev+1)
95 ✗ REAL fracd(klon, klev+1)
96 ✗ REAL xxx(klon, klev+1)
97 ✗ REAL larg_cons(klon, klev+1)
98 ✗ REAL larg_detr(klon, klev+1)
99 ✗ REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
100 ✗ REAL pu_therm(klon, klev), pv_therm(klon, klev)
101 ✗ REAL fm(klon, klev+1), entr(klon, klev)
102 ✗ REAL fmc(klon, klev+1)
103
104 CHARACTER (LEN=2) :: str2
105 CHARACTER (LEN=10) :: str10
106
107 CHARACTER (LEN=20) :: modname = 'thermcell2002'
108 CHARACTER (LEN=80) :: abort_message
109
110 LOGICAL vtest(klon), down
111
112 EXTERNAL scopy
113
114 INTEGER ncorrec, ll
115 SAVE ncorrec
116 DATA ncorrec/0/
117 !$OMP THREADPRIVATE(ncorrec)
118
119
120 ! -----------------------------------------------------------------------
121 ! initialisation:
122 ! ---------------
123
124 sorties = .TRUE.
125 ✗ IF (ngrid/=klon) THEN
126 ✗ PRINT *
127 ✗ PRINT *, 'STOP dans convadj'
128 ✗ PRINT *, 'ngrid =', ngrid
129 ✗ PRINT *, 'klon =', klon
130 END IF
131
132 ! -----------------------------------------------------------------------
133 ! incrementation eventuelle de tendances precedentes:
134 ! ---------------------------------------------------
135
136 ! print*,'0 OK convect8'
137
138 ✗ DO l = 1, nlay
139 ✗ DO ig = 1, ngrid
140 ✗ zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa
141 ✗ zh(ig, l) = pt(ig, l)/zpspsk(ig, l)
142 ✗ zu(ig, l) = pu(ig, l)
143 ✗ zv(ig, l) = pv(ig, l)
144 ✗ zo(ig, l) = po(ig, l)
145 ✗ ztv(ig, l) = zh(ig, l)*(1.+0.61*zo(ig,l))
146 END DO
147 END DO
148
149 ! print*,'1 OK convect8'
150 ! --------------------
151
152
153 ! + + + + + + + + + + +
154
155
156 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
157 ! wh,wt,wo ...
158
159 ! + + + + + + + + + + + zh,zu,zv,zo,rho
160
161
162 ! -------------------- zlev(1)
163 ! \\\\\\\\\\\\\\\\\\\
164
165
166 ! -----------------------------------------------------------------------
167 ! Calcul des altitudes des couches
168 ! -----------------------------------------------------------------------
169
170 ✗ IF (debut) THEN
171 ✗ flagdq = (iflag_thermals-1000)/100
172 ✗ dvdq = (iflag_thermals-(1000+flagdq*100))/10
173 ✗ IF (flagdq==2) dqimpl = -1
174 ✗ IF (flagdq==3) dqimpl = 1
175 ✗ debut = .FALSE.
176 END IF
177 ✗ PRINT *, 'TH flag th ', iflag_thermals, flagdq, dvdq, dqimpl
178
179 ✗ DO l = 2, nlay
180 ✗ DO ig = 1, ngrid
181 ✗ zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
182 END DO
183 END DO
184 ✗ DO ig = 1, ngrid
185 ✗ zlev(ig, 1) = 0.
186 ✗ zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
187 END DO
188 ✗ DO l = 1, nlay
189 ✗ DO ig = 1, ngrid
190 ✗ zlay(ig, l) = pphi(ig, l)/rg
191 END DO
192 END DO
193
194 ! print*,'2 OK convect8'
195 ! -----------------------------------------------------------------------
196 ! Calcul des densites
197 ! -----------------------------------------------------------------------
198
199 ✗ DO l = 1, nlay
200 ✗ DO ig = 1, ngrid
201 ✗ rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*zh(ig,l))
202 END DO
203 END DO
204
205 ✗ DO l = 2, nlay
206 ✗ DO ig = 1, ngrid
207 ✗ rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
208 END DO
209 END DO
210
211 ✗ DO k = 1, nlay
212 ✗ DO l = 1, nlay + 1
213 ✗ DO ig = 1, ngrid
214 ✗ wa(ig, k, l) = 0.
215 END DO
216 END DO
217 END DO
218
219 ! print*,'3 OK convect8'
220 ! ------------------------------------------------------------------
221 ! Calcul de w2, quarre de w a partir de la cape
222 ! a partir de w2, on calcule wa, vitesse de l'ascendance
223
224 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
225 ! w2 est stoke dans wa
226
227 ! ATTENTION: dans convect8, on n'utilise le calcule des wa
228 ! independants par couches que pour calculer l'entrainement
229 ! a la base et la hauteur max de l'ascendance.
230
231 ! Indicages:
232 ! l'ascendance provenant du niveau k traverse l'interface l avec
233 ! une vitesse wa(k,l).
234
235 ! --------------------
236
237 ! + + + + + + + + + +
238
239 ! wa(k,l) ---- -------------------- l
240 ! / ! /||\ + + + + + + + + + +
241 ! ||
242 ! || --------------------
243 ! ||
244 ! || + + + + + + + + + +
245 ! ||
246 ! || --------------------
247 ! ||__
248 ! |___ + + + + + + + + + + k
249
250 ! --------------------
251
252
253
254 ! ------------------------------------------------------------------
255
256
257 ✗ DO k = 1, nlay - 1
258 ✗ DO ig = 1, ngrid
259 ✗ wa(ig, k, k) = 0.
260 wa(ig, k, k+1) = 2.*rg*(ztv(ig,k)-ztv(ig,k+1))/ztv(ig, k+1)* &
261 ✗ (zlev(ig,k+1)-zlev(ig,k))
262 END DO
263 ✗ DO l = k + 1, nlay - 1
264 ✗ DO ig = 1, ngrid
265 wa(ig, k, l+1) = wa(ig, k, l) + 2.*rg*(ztv(ig,k)-ztv(ig,l))/ztv(ig, l &
266 ✗ )*(zlev(ig,l+1)-zlev(ig,l))
267 END DO
268 END DO
269 ✗ DO ig = 1, ngrid
270 ✗ wa(ig, k, nlay+1) = 0.
271 END DO
272 END DO
273
274 ! print*,'4 OK convect8'
275 ! Calcul de la couche correspondant a la hauteur du thermique
276 ✗ DO k = 1, nlay - 1
277 ✗ DO ig = 1, ngrid
278 ✗ lmax(ig, k) = k
279 END DO
280 ✗ DO l = nlay, k + 1, -1
281 ✗ DO ig = 1, ngrid
282 ✗ IF (wa(ig,k,l)<=1.E-10) lmax(ig, k) = l - 1
283 END DO
284 END DO
285 END DO
286
287 ! print*,'5 OK convect8'
288 ! Calcule du w max du thermique
289 ✗ DO k = 1, nlay
290 ✗ DO ig = 1, ngrid
291 ✗ wmax(ig, k) = 0.
292 END DO
293 END DO
294
295 ✗ DO k = 1, nlay - 1
296 ✗ DO l = k, nlay
297 ✗ DO ig = 1, ngrid
298 ✗ IF (l<=lmax(ig,k)) THEN
299 ✗ wa(ig, k, l) = sqrt(wa(ig,k,l))
300 ✗ wmax(ig, k) = max(wmax(ig,k), wa(ig,k,l))
301 ELSE
302 ✗ wa(ig, k, l) = 0.
303 END IF
304 END DO
305 END DO
306 END DO
307
308 ✗ DO k = 1, nlay - 1
309 ✗ DO ig = 1, ngrid
310 ✗ pu_therm(ig, k) = sqrt(wmax(ig,k))
311 ✗ pv_therm(ig, k) = sqrt(wmax(ig,k))
312 END DO
313 END DO
314
315 ! print*,'6 OK convect8'
316 ! Longueur caracteristique correspondant a la hauteur des thermiques.
317 ✗ DO ig = 1, ngrid
318 ✗ zmax(ig) = 500.
319 END DO
320 ! print*,'LMAX LMAX LMAX '
321 ✗ DO k = 1, nlay - 1
322 ✗ DO ig = 1, ngrid
323 ✗ zmax(ig) = max(zmax(ig), zlev(ig,lmax(ig,k))-zlev(ig,k))
324 END DO
325 ! print*,k,lmax(1,k)
326 END DO
327 ! print*,'ZMAX ZMAX ZMAX ',zmax
328 ! call dump2d(iim,jjm-1,zmax(2:ngrid-1),'ZMAX ')
329
330 ! print*,'OKl336'
331 ! Calcul de l'entrainement.
332 ! Le rapport d'aspect relie la largeur de l'ascendance a l'epaisseur
333 ! de la couche d'alimentation en partant du principe que la vitesse
334 ! maximum dans l'ascendance est la vitesse d'entrainement horizontale.
335 ✗ DO k = 1, nlay
336 ✗ DO ig = 1, ngrid
337 zzz = rho(ig, k)*wmax(ig, k)*(zlev(ig,k+1)-zlev(ig,k))/ &
338 ✗ (zmax(ig)*r_aspect)
339 ✗ IF (w2di==2) THEN
340 ✗ entr(ig, k) = entr(ig, k) + ptimestep*(zzz-entr(ig,k))/tho
341 ELSE
342 ✗ entr(ig, k) = zzz
343 END IF
344 ✗ ztva(ig, k) = ztv(ig, k)
345 END DO
346 END DO
347
348
349 ! print*,'7 OK convect8'
350 ✗ DO k = 1, klev + 1
351 ✗ DO ig = 1, ngrid
352 ✗ zw2(ig, k) = 0.
353 ✗ fmc(ig, k) = 0.
354 ✗ larg_cons(ig, k) = 0.
355 ✗ larg_detr(ig, k) = 0.
356 ✗ wa_moy(ig, k) = 0.
357 END DO
358 END DO
359
360 ! print*,'8 OK convect8'
361 ✗ DO ig = 1, ngrid
362 ✗ lmaxa(ig) = 1
363 ✗ lmix(ig) = 1
364 ✗ wmaxa(ig) = 0.
365 END DO
366
367
368 ! print*,'OKl372'
369 ✗ DO l = 1, nlay - 2
370 ✗ DO ig = 1, ngrid
371 ! if (zw2(ig,l).lt.1.e-10.and.ztv(ig,l).gt.ztv(ig,l+1)) then
372 ! print*,'COUCOU ',l,zw2(ig,l),ztv(ig,l),ztv(ig,l+1)
373 ✗ IF (zw2(ig,l)<1.E-10 .AND. ztv(ig,l)>ztv(ig,l+1) .AND. &
374 entr(ig,l)>1.E-10) THEN
375 ! print*,'COUCOU cas 1'
376 ! Initialisation de l'ascendance
377 ! lmix(ig)=1
378 ✗ ztva(ig, l) = ztv(ig, l)
379 ✗ fmc(ig, l) = 0.
380 ✗ fmc(ig, l+1) = entr(ig, l)
381 ✗ zw2(ig, l) = 0.
382 ! if (.not.ztv(ig,l+1).gt.150.) then
383 ! print*,'ig,l+1,ztv(ig,l+1)'
384 ! print*, ig,l+1,ztv(ig,l+1)
385 ! endif
386 zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
387 ✗ (zlev(ig,l+1)-zlev(ig,l))
388 ✗ larg_detr(ig, l) = 0.
389 ✗ ELSE IF (zw2(ig,l)>=1.E-10 .AND. fmc(ig,l)+entr(ig,l)>1.E-10) THEN
390 ! Incrementation...
391 ✗ fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
392 ! if (.not.fmc(ig,l+1).gt.1.e-15) then
393 ! print*,'ig,l+1,fmc(ig,l+1)'
394 ! print*, ig,l+1,fmc(ig,l+1)
395 ! print*,'Fmc ',(fmc(ig,ll),ll=1,klev+1)
396 ! print*,'W2 ',(zw2(ig,ll),ll=1,klev+1)
397 ! print*,'Tv ',(ztv(ig,ll),ll=1,klev)
398 ! print*,'Entr ',(entr(ig,ll),ll=1,klev)
399 ! endif
400 ztva(ig, l) = (fmc(ig,l)*ztva(ig,l-1)+entr(ig,l)*ztv(ig,l))/ &
401 ✗ fmc(ig, l+1)
402 ! mise a jour de la vitesse ascendante (l'air entraine de la couche
403 ! consideree commence avec une vitesse nulle).
404 zw2(ig, l+1) = zw2(ig, l)*(fmc(ig,l)/fmc(ig,l+1))**2 + &
405 ✗ 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
406 END IF
407 ✗ IF (zw2(ig,l+1)<0.) THEN
408 ✗ zw2(ig, l+1) = 0.
409 ✗ lmaxa(ig) = l
410 ELSE
411 ✗ wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
412 END IF
413 ✗ IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
414 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
415 ✗ lmix(ig) = l + 1
416 ✗ wmaxa(ig) = wa_moy(ig, l+1)
417 END IF
418 ! print*,'COUCOU cas 2 LMIX=',lmix(ig),wa_moy(ig,l+1),wmaxa(ig)
419 END DO
420 END DO
421
422 ! print*,'9 OK convect8'
423 ! print*,'WA1 ',wa_moy
424
425 ! determination de l'indice du debut de la mixed layer ou w decroit
426
427 ! calcul de la largeur de chaque ascendance dans le cas conservatif.
428 ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
429 ! d'une couche est �gale � la hauteur de la couche alimentante.
430 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
431 ! de la vitesse d'entrainement horizontal dans la couche alimentante.
432
433 ! print*,'OKl439'
434 ✗ DO l = 2, nlay
435 ✗ DO ig = 1, ngrid
436 ✗ IF (l<=lmaxa(ig)) THEN
437 ✗ zw = max(wa_moy(ig,l), 1.E-10)
438 ✗ larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw)
439 END IF
440 END DO
441 END DO
442
443 ✗ DO l = 2, nlay
444 ✗ DO ig = 1, ngrid
445 ✗ IF (l<=lmaxa(ig)) THEN
446 ! if (idetr.eq.0) then
447 ! cette option est finalement en dur.
448 ✗ larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
449 ! else if (idetr.eq.1) then
450 ! larg_detr(ig,l)=larg_cons(ig,l)
451 ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
452 ! else if (idetr.eq.2) then
453 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
454 ! s *sqrt(wa_moy(ig,l))
455 ! else if (idetr.eq.4) then
456 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
457 ! s *wa_moy(ig,l)
458 ! endif
459 END IF
460 END DO
461 END DO
462
463 ! print*,'10 OK convect8'
464 ! print*,'WA2 ',wa_moy
465 ! calcul de la fraction de la maille concern�e par l'ascendance en tenant
466 ! compte de l'epluchage du thermique.
467
468 ✗ DO l = 2, nlay
469 ✗ DO ig = 1, ngrid
470 ✗ IF (larg_cons(ig,l)>1.) THEN
471 ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
472 ✗ fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
473 ✗ IF (l>lmix(ig)) THEN
474 ✗ xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
475 ✗ IF (idetr==0) THEN
476 ✗ fraca(ig, l) = fraca(ig, lmix(ig))
477 ✗ ELSE IF (idetr==1) THEN
478 ✗ fraca(ig, l) = fraca(ig, lmix(ig))*xxx(ig, l)
479 ✗ ELSE IF (idetr==2) THEN
480 ✗ fraca(ig, l) = fraca(ig, lmix(ig))*(1.-(1.-xxx(ig,l))**2)
481 ELSE
482 ✗ fraca(ig, l) = fraca(ig, lmix(ig))*xxx(ig, l)**2
483 END IF
484 END IF
485 ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
486 ✗ fraca(ig, l) = max(fraca(ig,l), 0.)
487 ✗ fraca(ig, l) = min(fraca(ig,l), 0.5)
488 ✗ fracd(ig, l) = 1. - fraca(ig, l)
489 ✗ fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
490 ELSE
491 ! wa_moy(ig,l)=0.
492 ✗ fraca(ig, l) = 0.
493 ✗ fracc(ig, l) = 0.
494 ✗ fracd(ig, l) = 1.
495 END IF
496 END DO
497 END DO
498
499 ! print*,'11 OK convect8'
500 ! print*,'Ea3 ',wa_moy
501 ! ------------------------------------------------------------------
502 ! Calcul de fracd, wd
503 ! somme wa - wd = 0
504 ! ------------------------------------------------------------------
505
506
507 ✗ DO ig = 1, ngrid
508 ✗ fm(ig, 1) = 0.
509 ✗ fm(ig, nlay+1) = 0.
510 END DO
511
512 ✗ DO l = 2, nlay
513 ✗ DO ig = 1, ngrid
514 ✗ fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
515 END DO
516 ✗ DO ig = 1, ngrid
517 ✗ IF (fracd(ig,l)<0.1) THEN
518 ✗ abort_message = 'fracd trop petit'
519 ✗ CALL abort_physic(modname, abort_message, 1)
520 ELSE
521 ! vitesse descendante "diagnostique"
522 ✗ wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
523 END IF
524 END DO
525 END DO
526
527 ✗ DO l = 1, nlay
528 ✗ DO ig = 1, ngrid
529 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
530 ✗ masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
531 END DO
532 END DO
533
534 ! print*,'12 OK convect8'
535 ! print*,'WA4 ',wa_moy
536 ! c------------------------------------------------------------------
537 ! calcul du transport vertical
538 ! ------------------------------------------------------------------
539
540 GO TO 4444
541 ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
542 DO l = 2, nlay - 1
543 DO ig = 1, ngrid
544 IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( &
545 ig,l+1)) THEN
546 ! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
547 ! s ,fm(ig,l+1)*ptimestep
548 ! s ,' M=',masse(ig,l),masse(ig,l+1)
549 END IF
550 END DO
551 END DO
552
553 DO l = 1, nlay
554 DO ig = 1, ngrid
555 IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
556 ! print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
557 ! s ,entr(ig,l)*ptimestep
558 ! s ,' M=',masse(ig,l)
559 END IF
560 END DO
561 END DO
562
563 DO l = 1, nlay
564 DO ig = 1, ngrid
565 IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
566 ! print*,'WARN!!! fm exagere ig=',ig,' l=',l
567 ! s ,' FM=',fm(ig,l)
568 END IF
569 IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
570 ! print*,'WARN!!! masse exagere ig=',ig,' l=',l
571 ! s ,' M=',masse(ig,l)
572 ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
573 ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
574 ! print*,'zlev(ig,l+1),zlev(ig,l)'
575 ! s ,zlev(ig,l+1),zlev(ig,l)
576 ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
577 ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
578 END IF
579 IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
580 ! print*,'WARN!!! entr exagere ig=',ig,' l=',l
581 ! s ,' E=',entr(ig,l)
582 END IF
583 END DO
584 END DO
585
586 4444 CONTINUE
587 ! print*,'OK 444 '
588
589 ✗ IF (w2di==1) THEN
590 ✗ fm0 = fm0 + ptimestep*(fm-fm0)/tho
591 ✗ entr0 = entr0 + ptimestep*(entr-entr0)/tho
592 ELSE
593 ✗ fm0 = fm
594 ✗ entr0 = entr
595 END IF
596
597 ✗ IF (flagdq==0) THEN
598 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, &
599 ✗ zha)
600 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, &
601 ✗ zoa)
602 ✗ PRINT *, 'THERMALS OPT 1'
603 ✗ ELSE IF (flagdq==1) THEN
604 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
605 ✗ zdhadj, zha)
606 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
607 ✗ pdoadj, zoa)
608 ✗ PRINT *, 'THERMALS OPT 2'
609 ELSE
610 CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zh, &
611 ✗ zdhadj, zha, lev_out)
612 CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zo, &
613 ✗ pdoadj, zoa, lev_out)
614 ✗ PRINT *, 'THERMALS OPT 3', dqimpl
615 END IF
616
617 ✗ PRINT *, 'TH VENT ', dvdq
618 ✗ IF (dvdq==0) THEN
619 ! print*,'TH VENT OK ',dvdq
620 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
621 ✗ zua)
622 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
623 ✗ zva)
624 ✗ ELSE IF (dvdq==1) THEN
625 CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
626 ✗ zu, zv, pduadj, pdvadj, zua, zva)
627 ✗ ELSE IF (dvdq==2) THEN
628 CALL thermcell_dv2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, &
629 ✗ zmax, zu, zv, pduadj, pdvadj, zua, zva, lev_out)
630 ✗ ELSE IF (dvdq==3) THEN
631 CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zu, &
632 ✗ pduadj, zua, lev_out)
633 CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zv, &
634 ✗ pdvadj, zva, lev_out)
635 END IF
636
637 ! CALL writefield_phy('duadj',pduadj,klev)
638
639 ✗ DO l = 1, nlay
640 ✗ DO ig = 1, ngrid
641 ✗ zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
642 ✗ zf2 = zf/(1.-zf)
643 ✗ thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l))**2
644 ✗ wth2(ig, l) = zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
645 END DO
646 END DO
647
648
649
650 ! print*,'13 OK convect8'
651 ! print*,'WA5 ',wa_moy
652 ✗ DO l = 1, nlay
653 ✗ DO ig = 1, ngrid
654 ✗ pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l)
655 END DO
656 END DO
657
658
659 ! do l=1,nlay
660 ! do ig=1,ngrid
661 ! if(abs(pdtadj(ig,l))*86400..gt.500.) then
662 ! print*,'WARN!!! ig=',ig,' l=',l
663 ! s ,' pdtadj=',pdtadj(ig,l)
664 ! endif
665 ! if(abs(pdoadj(ig,l))*86400..gt.1.) then
666 ! print*,'WARN!!! ig=',ig,' l=',l
667 ! s ,' pdoadj=',pdoadj(ig,l)
668 ! endif
669 ! enddo
670 ! enddo
671
672 ! print*,'14 OK convect8'
673 ! ------------------------------------------------------------------
674 ! Calculs pour les sorties
675 ! ------------------------------------------------------------------
676
677 IF (sorties) THEN
678 DO l = 1, nlay
679 DO ig = 1, ngrid
680 zla(ig, l) = (1.-fracd(ig,l))*zmax(ig)
681 zld(ig, l) = fracd(ig, l)*zmax(ig)
682 IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ &
683 (1.-fracd(ig,l))
684 END DO
685 END DO
686
687 DO l = 1, nlay
688 DO ig = 1, ngrid
689 detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1)
690 IF (detr(ig,l)<0.) THEN
691 entr(ig, l) = entr(ig, l) - detr(ig, l)
692 detr(ig, l) = 0.
693 ! print*,'WARNING !!! detrainement negatif ',ig,l
694 END IF
695 END DO
696 END DO
697 END IF
698
699 ! print*,'15 OK convect8'
700
701
702 ! if(wa_moy(1,4).gt.1.e-10) stop
703
704 ! print*,'19 OK convect8'
705 ✗ RETURN
706 END SUBROUTINE thermcell_2002
707
708 ✗ SUBROUTINE thermcell_cld(ngrid, nlay, ptimestep, pplay, pplev, pphi, zlev, &
709 ✗ debut, pu, pv, pt, po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0, zqla, &
710 ✗ lmax, zmax_sec, wmax_sec, zw_sec, lmix_sec, ratqscth, ratqsdiff & ! s
711 ! ,pu_therm,pv_therm
712 , r_aspect, l_mix, w2di, tho)
713
714 USE dimphy
715 IMPLICIT NONE
716
717 ! =======================================================================
718
719 ! Calcul du transport verticale dans la couche limite en presence
720 ! de "thermiques" explicitement representes
721
722 ! R��criture � partir d'un listing papier � Habas, le 14/02/00
723
724 ! le thermique est suppos� homog�ne et dissip� par m�lange avec
725 ! son environnement. la longueur l_mix contr�le l'efficacit� du
726 ! m�lange
727
728 ! Le calcul du transport des diff�rentes esp�ces se fait en prenant
729 ! en compte:
730 ! 1. un flux de masse montant
731 ! 2. un flux de masse descendant
732 ! 3. un entrainement
733 ! 4. un detrainement
734
735 ! =======================================================================
736
737 ! -----------------------------------------------------------------------
738 ! declarations:
739 ! -------------
740
741 include "YOMCST.h"
742 include "YOETHF.h"
743 include "FCTTRE.h"
744
745 ! arguments:
746 ! ----------
747
748 INTEGER ngrid, nlay, w2di
749 REAL tho
750 REAL ptimestep, l_mix, r_aspect
751 REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
752 REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
753 REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
754 REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
755 REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
756 REAL pphi(ngrid, nlay)
757
758 INTEGER idetr
759 SAVE idetr
760 DATA idetr/3/
761 !$OMP THREADPRIVATE(idetr)
762
763 ! local:
764 ! ------
765
766 ✗ INTEGER ig, k, l, lmaxa(klon), lmix(klon)
767 REAL zsortie1d(klon)
768 ! CR: on remplace lmax(klon,klev+1)
769 ✗ INTEGER lmax(klon), lmin(klon), lentr(klon)
770 ✗ REAL linter(klon)
771 ✗ REAL zmix(klon), fracazmix(klon)
772 REAL alpha
773 SAVE alpha
774 DATA alpha/1./
775 !$OMP THREADPRIVATE(alpha)
776
777 ! RC
778 ✗ REAL zmax(klon), zw, zz, zw2(klon, klev+1), ztva(klon, klev), zzz
779 REAL zmax_sec(klon)
780 ✗ REAL zmax_sec2(klon)
781 REAL zw_sec(klon, klev+1)
782 INTEGER lmix_sec(klon)
783 ✗ REAL w_est(klon, klev+1)
784 ! on garde le zmax du pas de temps precedent
785 ! real zmax0(klon)
786 ! save zmax0
787 ! real zmix0(klon)
788 ! save zmix0
789 REAL, SAVE, ALLOCATABLE :: zmax0(:), zmix0(:)
790 !$OMP THREADPRIVATE(zmax0, zmix0)
791
792 ✗ REAL zlev(klon, klev+1), zlay(klon, klev)
793 ✗ REAL deltaz(klon, klev)
794 ✗ REAL zh(klon, klev), zdhadj(klon, klev)
795 ✗ REAL zthl(klon, klev), zdthladj(klon, klev)
796 ✗ REAL ztv(klon, klev)
797 ✗ REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
798 ✗ REAL zl(klon, klev)
799 REAL wh(klon, klev+1)
800 REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
801 REAL zla(klon, klev+1)
802 REAL zwa(klon, klev+1)
803 REAL zld(klon, klev+1)
804 REAL zwd(klon, klev+1)
805 REAL zsortie(klon, klev)
806 ✗ REAL zva(klon, klev)
807 ✗ REAL zua(klon, klev)
808 ✗ REAL zoa(klon, klev)
809
810 ✗ REAL zta(klon, klev)
811 ✗ REAL zha(klon, klev)
812 ✗ REAL wa_moy(klon, klev+1)
813 ✗ REAL fraca(klon, klev+1)
814 ✗ REAL fracc(klon, klev+1)
815 REAL zf, zf2
816 REAL thetath2(klon, klev), wth2(klon, klev), wth3(klon, klev)
817 REAL q2(klon, klev)
818 ✗ REAL dtheta(klon, klev)
819 ! common/comtherm/thetath2,wth2
820
821 REAL ratqscth(klon, klev)
822 REAL sum
823 REAL sumdiff
824 REAL ratqsdiff(klon, klev)
825 REAL count_time
826 INTEGER ialt
827
828 LOGICAL sorties
829 ✗ REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
830 ✗ REAL zpspsk(klon, klev)
831
832 ! real wmax(klon,klev),wmaxa(klon)
833 ✗ REAL wmax(klon), wmaxa(klon)
834 REAL wmax_sec(klon)
835 ✗ REAL wmax_sec2(klon)
836 ✗ REAL wa(klon, klev, klev+1)
837 ✗ REAL wd(klon, klev+1)
838 REAL larg_part(klon, klev, klev+1)
839 ✗ REAL fracd(klon, klev+1)
840 ✗ REAL xxx(klon, klev+1)
841 ✗ REAL larg_cons(klon, klev+1)
842 ✗ REAL larg_detr(klon, klev+1)
843 ✗ REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
844 REAL massetot(klon, klev)
845 ✗ REAL detr0(klon, klev)
846 ✗ REAL alim0(klon, klev)
847 REAL pu_therm(klon, klev), pv_therm(klon, klev)
848 ✗ REAL fm(klon, klev+1), entr(klon, klev)
849 ✗ REAL fmc(klon, klev+1)
850
851 REAL zcor, zdelta, zcvm5, qlbef
852 ✗ REAL tbef(klon), qsatbef(klon)
853 REAL dqsat_dt, dt, num, denom
854 REAL reps, rlvcp, ddt0
855 ✗ REAL ztla(klon, klev), zqla(klon, klev), zqta(klon, klev)
856 ! CR niveau de condensation
857 REAL nivcon(klon)
858 REAL zcon(klon)
859 ✗ REAL zqsat(klon, klev)
860 ✗ REAL zqsatth(klon, klev)
861 PARAMETER (ddt0=.01)
862
863
864 ! CR:nouvelles variables
865 ✗ REAL f_star(klon, klev+1), entr_star(klon, klev)
866 ✗ REAL detr_star(klon, klev)
867 ✗ REAL alim_star_tot(klon), alim_star2(klon)
868 ✗ REAL entr_star_tot(klon)
869 ✗ REAL detr_star_tot(klon)
870 ✗ REAL alim_star(klon, klev)
871 ✗ REAL alim(klon, klev)
872 ✗ REAL nu(klon, klev)
873 ✗ REAL nu_e(klon, klev)
874 REAL nu_min
875 REAL nu_max
876 REAL nu_r
877 ✗ REAL f(klon)
878 ! real f(klon), f0(klon)
879 ! save f0
880 REAL, SAVE, ALLOCATABLE :: f0(:)
881 !$OMP THREADPRIVATE(f0)
882
883 REAL f_old
884 ✗ REAL zlevinter(klon)
885 LOGICAL, SAVE :: first = .TRUE.
886 !$OMP THREADPRIVATE(first)
887 ! data first /.false./
888 ! save first
889 LOGICAL nuage
890 ! save nuage
891 LOGICAL boucle
892 LOGICAL therm
893 LOGICAL debut
894 LOGICAL rale
895 ✗ INTEGER test(klon)
896 INTEGER signe_zw2
897 ! RC
898
899 CHARACTER *2 str2
900 CHARACTER *10 str10
901
902 CHARACTER (LEN=20) :: modname = 'thermcell_cld'
903 CHARACTER (LEN=80) :: abort_message
904
905 LOGICAL vtest(klon), down
906 ✗ LOGICAL zsat(klon)
907
908 EXTERNAL scopy
909
910 INTEGER ncorrec, ll
911 SAVE ncorrec
912 DATA ncorrec/0/
913 !$OMP THREADPRIVATE(ncorrec)
914
915
916
917 ! -----------------------------------------------------------------------
918 ! initialisation:
919 ! ---------------
920
921 ✗ IF (first) THEN
922 ✗ ALLOCATE (zmix0(klon))
923 ✗ ALLOCATE (zmax0(klon))
924 ✗ ALLOCATE (f0(klon))
925 ✗ first = .FALSE.
926 END IF
927
928 sorties = .FALSE.
929 ! print*,'NOUVEAU DETR PLUIE '
930 ✗ IF (ngrid/=klon) THEN
931 ✗ PRINT *
932 ✗ PRINT *, 'STOP dans convadj'
933 ✗ PRINT *, 'ngrid =', ngrid
934 ✗ PRINT *, 'klon =', klon
935 END IF
936
937 ! Initialisation
938 ✗ rlvcp = rlvtt/rcpd
939 reps = rd/rv
940 ! initialisations de zqsat
941 ✗ DO ll = 1, nlay
942 ✗ DO ig = 1, ngrid
943 ✗ zqsat(ig, ll) = 0.
944 ✗ zqsatth(ig, ll) = 0.
945 END DO
946 END DO
947
948 ! on met le first a true pour le premier passage de la journ�e
949 ✗ DO ig = 1, klon
950 ✗ test(ig) = 0
951 END DO
952 ✗ IF (debut) THEN
953 ✗ DO ig = 1, klon
954 ✗ test(ig) = 1
955 ✗ f0(ig) = 0.
956 ✗ zmax0(ig) = 0.
957 END DO
958 END IF
959 ✗ DO ig = 1, klon
960 ✗ IF ((.NOT. debut) .AND. (f0(ig)<1.E-10)) THEN
961 ✗ test(ig) = 1
962 END IF
963 END DO
964 ! do ig=1,klon
965 ! print*,'test(ig)',test(ig),zmax0(ig)
966 ! enddo
967 nuage = .FALSE.
968 ! -----------------------------------------------------------------------
969 ! AM Calcul de T,q,ql a partir de Tl et qT
970 ! ---------------------------------------------------
971
972 ! Pr Tprec=Tl calcul de qsat
973 ! Si qsat>qT T=Tl, q=qT
974 ! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
975 ! On cherche DDT < DDT0
976
977 ! defaut
978 ✗ DO ll = 1, nlay
979 ✗ DO ig = 1, ngrid
980 ✗ zo(ig, ll) = po(ig, ll)
981 ✗ zl(ig, ll) = 0.
982 ✗ zh(ig, ll) = pt(ig, ll)
983 END DO
984 END DO
985 ✗ DO ig = 1, ngrid
986 ✗ zsat(ig) = .FALSE.
987 END DO
988
989
990 ✗ DO ll = 1, nlay
991 ! les points insatures sont definitifs
992 ✗ DO ig = 1, ngrid
993 ✗ tbef(ig) = pt(ig, ll)
994 ✗ zdelta = max(0., sign(1.,rtt-tbef(ig)))
995 ✗ qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, ll)
996 ✗ qsatbef(ig) = min(0.5, qsatbef(ig))
997 ✗ zcor = 1./(1.-retv*qsatbef(ig))
998 ✗ qsatbef(ig) = qsatbef(ig)*zcor
999 ✗ zsat(ig) = (max(0.,po(ig,ll)-qsatbef(ig))>1.E-10)
1000 END DO
1001
1002 ✗ DO ig = 1, ngrid
1003 ✗ IF (zsat(ig) .AND. (1==1)) THEN
1004 ✗ qlbef = max(0., po(ig,ll)-qsatbef(ig))
1005 ! si sature: ql est surestime, d'ou la sous-relax
1006 ✗ dt = 0.5*rlvcp*qlbef
1007 ! write(18,*),'DT0=',DT
1008 ! on pourra enchainer 2 ou 3 calculs sans Do while
1009 ✗ DO WHILE (abs(dt)>ddt0)
1010 ! il faut verifier si c,a conserve quand on repasse en insature ...
1011 ✗ tbef(ig) = tbef(ig) + dt
1012 ✗ zdelta = max(0., sign(1.,rtt-tbef(ig)))
1013 ✗ qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, ll)
1014 ✗ qsatbef(ig) = min(0.5, qsatbef(ig))
1015 ✗ zcor = 1./(1.-retv*qsatbef(ig))
1016 ✗ qsatbef(ig) = qsatbef(ig)*zcor
1017 ! on veut le signe de qlbef
1018 ✗ qlbef = po(ig, ll) - qsatbef(ig)
1019 zdelta = max(0., sign(1.,rtt-tbef(ig)))
1020 ✗ zcvm5 = r5les*(1.-zdelta) + r5ies*zdelta
1021 ✗ zcor = 1./(1.-retv*qsatbef(ig))
1022 ✗ dqsat_dt = foede(tbef(ig), zdelta, zcvm5, qsatbef(ig), zcor)
1023 ✗ num = -tbef(ig) + pt(ig, ll) + rlvcp*qlbef
1024 ✗ denom = 1. + rlvcp*dqsat_dt
1025 ✗ IF (denom<1.E-10) THEN
1026 ✗ PRINT *, 'pb denom'
1027 END IF
1028 ✗ dt = num/denom
1029 END DO
1030 ! on ecrit de maniere conservative (sat ou non)
1031 ✗ zl(ig, ll) = max(0., qlbef)
1032 ! T = Tl +Lv/Cp ql
1033 ✗ zh(ig, ll) = pt(ig, ll) + rlvcp*zl(ig, ll)
1034 ✗ zo(ig, ll) = po(ig, ll) - zl(ig, ll)
1035 END IF
1036 ! on ecrit zqsat
1037 ✗ zqsat(ig, ll) = qsatbef(ig)
1038 END DO
1039 END DO
1040 ! AM fin
1041
1042 ! -----------------------------------------------------------------------
1043 ! incrementation eventuelle de tendances precedentes:
1044 ! ---------------------------------------------------
1045
1046 ! print*,'0 OK convect8'
1047
1048 ✗ DO l = 1, nlay
1049 ✗ DO ig = 1, ngrid
1050 ✗ zpspsk(ig, l) = (pplay(ig,l)/100000.)**rkappa
1051 ! zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
1052 ! zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
1053 ✗ zu(ig, l) = pu(ig, l)
1054 ✗ zv(ig, l) = pv(ig, l)
1055 ! zo(ig,l)=po(ig,l)
1056 ! ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
1057 ! AM attention zh est maintenant le profil de T et plus le profil de
1058 ! theta !
1059
1060 ! T-> Theta
1061 ✗ ztv(ig, l) = zh(ig, l)/zpspsk(ig, l)
1062 ! AM Theta_v
1063 ✗ ztv(ig, l) = ztv(ig, l)*(1.+retv*(zo(ig,l))-zl(ig,l))
1064 ! AM Thetal
1065 ✗ zthl(ig, l) = pt(ig, l)/zpspsk(ig, l)
1066
1067 END DO
1068 END DO
1069
1070 ! print*,'1 OK convect8'
1071 ! --------------------
1072
1073
1074 ! + + + + + + + + + + +
1075
1076
1077 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
1078 ! wh,wt,wo ...
1079
1080 ! + + + + + + + + + + + zh,zu,zv,zo,rho
1081
1082
1083 ! -------------------- zlev(1)
1084 ! \\\\\\\\\\\\\\\\\\\
1085
1086
1087 ! -----------------------------------------------------------------------
1088 ! Calcul des altitudes des couches
1089 ! -----------------------------------------------------------------------
1090
1091 ✗ DO l = 2, nlay
1092 ✗ DO ig = 1, ngrid
1093 ✗ zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
1094 END DO
1095 END DO
1096 ✗ DO ig = 1, ngrid
1097 ✗ zlev(ig, 1) = 0.
1098 ✗ zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
1099 END DO
1100 ✗ DO l = 1, nlay
1101 ✗ DO ig = 1, ngrid
1102 ✗ zlay(ig, l) = pphi(ig, l)/rg
1103 END DO
1104 END DO
1105 ! calcul de deltaz
1106 ✗ DO l = 1, nlay
1107 ✗ DO ig = 1, ngrid
1108 ✗ deltaz(ig, l) = zlev(ig, l+1) - zlev(ig, l)
1109 END DO
1110 END DO
1111
1112 ! print*,'2 OK convect8'
1113 ! -----------------------------------------------------------------------
1114 ! Calcul des densites
1115 ! -----------------------------------------------------------------------
1116
1117 ✗ DO l = 1, nlay
1118 ✗ DO ig = 1, ngrid
1119 ! rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
1120 ✗ rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*ztv(ig,l))
1121 END DO
1122 END DO
1123
1124 ✗ DO l = 2, nlay
1125 ✗ DO ig = 1, ngrid
1126 ✗ rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
1127 END DO
1128 END DO
1129
1130 ✗ DO k = 1, nlay
1131 ✗ DO l = 1, nlay + 1
1132 ✗ DO ig = 1, ngrid
1133 ✗ wa(ig, k, l) = 0.
1134 END DO
1135 END DO
1136 END DO
1137 ! Cr:ajout:calcul de la masse
1138 ✗ DO l = 1, nlay
1139 ✗ DO ig = 1, ngrid
1140 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1141 ✗ masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
1142 END DO
1143 END DO
1144 ! print*,'3 OK convect8'
1145 ! ------------------------------------------------------------------
1146 ! Calcul de w2, quarre de w a partir de la cape
1147 ! a partir de w2, on calcule wa, vitesse de l'ascendance
1148
1149 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
1150 ! w2 est stoke dans wa
1151
1152 ! ATTENTION: dans convect8, on n'utilise le calcule des wa
1153 ! independants par couches que pour calculer l'entrainement
1154 ! a la base et la hauteur max de l'ascendance.
1155
1156 ! Indicages:
1157 ! l'ascendance provenant du niveau k traverse l'interface l avec
1158 ! une vitesse wa(k,l).
1159
1160 ! --------------------
1161
1162 ! + + + + + + + + + +
1163
1164 ! wa(k,l) ---- -------------------- l
1165 ! / ! /||\ + + + + + + + + + +
1166 ! ||
1167 ! || --------------------
1168 ! ||
1169 ! || + + + + + + + + + +
1170 ! ||
1171 ! || --------------------
1172 ! ||__
1173 ! |___ + + + + + + + + + + k
1174
1175 ! --------------------
1176
1177
1178
1179 ! ------------------------------------------------------------------
1180
1181 ! CR: ponderation entrainement des couches instables
1182 ! def des alim_star tels que alim=f*alim_star
1183 ✗ DO l = 1, klev
1184 ✗ DO ig = 1, ngrid
1185 ✗ alim_star(ig, l) = 0.
1186 ✗ alim(ig, l) = 0.
1187 END DO
1188 END DO
1189 ! determination de la longueur de la couche d entrainement
1190 ✗ DO ig = 1, ngrid
1191 ✗ lentr(ig) = 1
1192 END DO
1193
1194 ! on ne considere que les premieres couches instables
1195 therm = .FALSE.
1196 ✗ DO k = nlay - 2, 1, -1
1197 ✗ DO ig = 1, ngrid
1198 ✗ IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN
1199 ✗ lentr(ig) = k + 1
1200 therm = .TRUE.
1201 END IF
1202 END DO
1203 END DO
1204
1205 ! determination du lmin: couche d ou provient le thermique
1206 ✗ DO ig = 1, ngrid
1207 ✗ lmin(ig) = 1
1208 END DO
1209 ✗ DO ig = 1, ngrid
1210 ✗ DO l = nlay, 2, -1
1211 ✗ IF (ztv(ig,l-1)>ztv(ig,l)) THEN
1212 ✗ lmin(ig) = l - 1
1213 END IF
1214 END DO
1215 END DO
1216
1217 ! definition de l'entrainement des couches
1218 ✗ DO l = 1, klev - 1
1219 ✗ DO ig = 1, ngrid
1220 ✗ IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<lentr(ig)) THEN
1221 ! def possibles pour alim_star: zdthetadz, dthetadz, zdtheta
1222 alim_star(ig, l) = max((ztv(ig,l)-ztv(ig,l+1)), 0.) & ! s
1223 ! *(zlev(ig,l+1)-zlev(ig,l))
1224 ✗ *sqrt(zlev(ig,l+1))
1225 ! alim_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
1226 ! s /zlev(ig,lentr(ig)+2)))**(3./2.)
1227 END IF
1228 END DO
1229 END DO
1230
1231 ! pas de thermique si couche 1 stable
1232 ✗ DO ig = 1, ngrid
1233 ! if (lmin(ig).gt.1) then
1234 ! CRnouveau test
1235 ✗ IF (alim_star(ig,1)<1.E-10) THEN
1236 ✗ DO l = 1, klev
1237 ✗ alim_star(ig, l) = 0.
1238 END DO
1239 END IF
1240 END DO
1241 ! calcul de l entrainement total
1242 ✗ DO ig = 1, ngrid
1243 ✗ alim_star_tot(ig) = 0.
1244 ✗ entr_star_tot(ig) = 0.
1245 ✗ detr_star_tot(ig) = 0.
1246 END DO
1247 ✗ DO ig = 1, ngrid
1248 ✗ DO k = 1, klev
1249 ✗ alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, k)
1250 END DO
1251 END DO
1252
1253 ! Calcul entrainement normalise
1254 ✗ DO ig = 1, ngrid
1255 ✗ IF (alim_star_tot(ig)>1.E-10) THEN
1256 ! do l=1,lentr(ig)
1257 ✗ DO l = 1, klev
1258 ! def possibles pour entr_star: zdthetadz, dthetadz, zdtheta
1259 ✗ alim_star(ig, l) = alim_star(ig, l)/alim_star_tot(ig)
1260 END DO
1261 END IF
1262 END DO
1263
1264 ! print*,'fin calcul alim_star'
1265
1266 ! AM:initialisations
1267 ✗ DO k = 1, nlay
1268 ✗ DO ig = 1, ngrid
1269 ✗ ztva(ig, k) = ztv(ig, k)
1270 ✗ ztla(ig, k) = zthl(ig, k)
1271 ✗ zqla(ig, k) = 0.
1272 ✗ zqta(ig, k) = po(ig, k)
1273 ✗ zsat(ig) = .FALSE.
1274 END DO
1275 END DO
1276 ✗ DO k = 1, klev
1277 ✗ DO ig = 1, ngrid
1278 ✗ detr_star(ig, k) = 0.
1279 ✗ entr_star(ig, k) = 0.
1280 ✗ detr(ig, k) = 0.
1281 ✗ entr(ig, k) = 0.
1282 END DO
1283 END DO
1284 ! print*,'7 OK convect8'
1285 ✗ DO k = 1, klev + 1
1286 ✗ DO ig = 1, ngrid
1287 ✗ zw2(ig, k) = 0.
1288 ✗ fmc(ig, k) = 0.
1289 ! CR
1290 ✗ f_star(ig, k) = 0.
1291 ! RC
1292 ✗ larg_cons(ig, k) = 0.
1293 ✗ larg_detr(ig, k) = 0.
1294 ✗ wa_moy(ig, k) = 0.
1295 END DO
1296 END DO
1297
1298 ! n print*,'8 OK convect8'
1299 ✗ DO ig = 1, ngrid
1300 ✗ linter(ig) = 1.
1301 ✗ lmaxa(ig) = 1
1302 ✗ lmix(ig) = 1
1303 ✗ wmaxa(ig) = 0.
1304 END DO
1305
1306 ✗ nu_min = l_mix
1307 ✗ nu_max = 1000.
1308 ! do ig=1,ngrid
1309 ! nu_max=wmax_sec(ig)
1310 ! enddo
1311 ✗ DO ig = 1, ngrid
1312 ✗ DO k = 1, klev
1313 ✗ nu(ig, k) = 0.
1314 ✗ nu_e(ig, k) = 0.
1315 END DO
1316 END DO
1317 ! Calcul de l'exc�s de temp�rature du � la diffusion turbulente
1318 ✗ DO ig = 1, ngrid
1319 ✗ DO l = 1, klev
1320 ✗ dtheta(ig, l) = 0.
1321 END DO
1322 END DO
1323 ✗ DO ig = 1, ngrid
1324 ✗ DO l = 1, lentr(ig) - 1
1325 dtheta(ig, l) = sqrt(10.*0.4*zlev(ig,l+1)**2*1.*((ztv(ig,l+1)- &
1326 ✗ ztv(ig,l))/(zlev(ig,l+1)-zlev(ig,l)))**2)
1327 END DO
1328 END DO
1329 ! do l=1,nlay-2
1330 ✗ DO l = 1, klev - 1
1331 ✗ DO ig = 1, ngrid
1332 ✗ IF (ztv(ig,l)>ztv(ig,l+1) .AND. alim_star(ig,l)>1.E-10 .AND. &
1333 zw2(ig,l)<1E-10) THEN
1334 ! AM
1335 ! test:on rajoute un exc�s de T dans couche alim
1336 ! ztla(ig,l)=zthl(ig,l)+dtheta(ig,l)
1337 ✗ ztla(ig, l) = zthl(ig, l)
1338 ! test: on rajoute un exc�s de q dans la couche alim
1339 ! zqta(ig,l)=po(ig,l)+0.001
1340 ✗ zqta(ig, l) = po(ig, l)
1341 ✗ zqla(ig, l) = zl(ig, l)
1342 ! AM
1343 ✗ f_star(ig, l+1) = alim_star(ig, l)
1344 ! test:calcul de dteta
1345 zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
1346 ✗ (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
1347 ✗ w_est(ig, l+1) = zw2(ig, l+1)
1348 ✗ larg_detr(ig, l) = 0.
1349 ! print*,'coucou boucle 1'
1350 ✗ ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+alim_star(ig, &
1351 l))>1.E-10) THEN
1352 ! print*,'coucou boucle 2'
1353 ! estimation du detrainement a partir de la geometrie du pas
1354 ! precedent
1355 ✗ IF ((test(ig)==1) .OR. ((.NOT. debut) .AND. (f0(ig)<1.E-10))) THEN
1356 ✗ detr_star(ig, l) = 0.
1357 ✗ entr_star(ig, l) = 0.
1358 ! print*,'coucou test(ig)',test(ig),f0(ig),zmax0(ig)
1359 ELSE
1360 ! print*,'coucou debut detr'
1361 ! tests sur la definition du detr
1362 ✗ IF (zqla(ig,l-1)>1.E-10) THEN
1363 nuage = .TRUE.
1364 END IF
1365
1366 w_est(ig, l+1) = zw2(ig, l)*((f_star(ig,l))**2)/(f_star(ig,l)+ &
1367 alim_star(ig,l))**2 + 2.*rg*(ztva(ig,l-1)-ztv(ig,l))/ztv(ig, l)*( &
1368 ✗ zlev(ig,l+1)-zlev(ig,l))
1369 ✗ IF (w_est(ig,l+1)<0.) THEN
1370 ✗ w_est(ig, l+1) = zw2(ig, l)
1371 END IF
1372 ✗ IF (l>2) THEN
1373 IF ((w_est(ig,l+1)>w_est(ig,l)) .AND. (zlev(ig, &
1374 ✗ l+1)<zmax_sec(ig)) .AND. (zqla(ig,l-1)<1.E-10)) THEN
1375 detr_star(ig, l) = max(0., (rhobarz(ig, &
1376 l+1)*sqrt(w_est(ig,l+1))*sqrt(nu(ig,l)* &
1377 zlev(ig,l+1))-rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(nu(ig,l)* &
1378 ✗ zlev(ig,l)))/(r_aspect*zmax_sec(ig)))
1379 ✗ ELSE IF ((zlev(ig,l+1)<zmax_sec(ig)) .AND. (zqla(ig, &
1380 l-1)<1.E-10)) THEN
1381 detr_star(ig, l) = -f0(ig)*f_star(ig, lmix(ig))/(rhobarz(ig, &
1382 lmix(ig))*wmaxa(ig))*(rhobarz(ig,l+1)*sqrt(w_est(ig, &
1383 l+1))*((zmax_sec(ig)-zlev(ig,l+1))/((zmax_sec(ig)-zlev(ig, &
1384 lmix(ig)))))**2.-rhobarz(ig,l)*sqrt(w_est(ig, &
1385 l))*((zmax_sec(ig)-zlev(ig,l))/((zmax_sec(ig)-zlev(ig,lmix(ig &
1386 ✗ )))))**2.)
1387 ELSE
1388 detr_star(ig, l) = 0.002*f0(ig)*f_star(ig, l)* &
1389 ✗ (zlev(ig,l+1)-zlev(ig,l))
1390
1391 END IF
1392 ELSE
1393 ✗ detr_star(ig, l) = 0.
1394 END IF
1395
1396 ✗ detr_star(ig, l) = detr_star(ig, l)/f0(ig)
1397 ✗ IF (nuage) THEN
1398 ✗ entr_star(ig, l) = 0.4*detr_star(ig, l)
1399 ELSE
1400 ✗ entr_star(ig, l) = 0.4*detr_star(ig, l)
1401 END IF
1402
1403 ✗ IF ((detr_star(ig,l))>f_star(ig,l)) THEN
1404 ✗ detr_star(ig, l) = f_star(ig, l)
1405 ! entr_star(ig,l)=0.
1406 END IF
1407
1408 ✗ IF ((l<lentr(ig))) THEN
1409 ✗ entr_star(ig, l) = 0.
1410 ! detr_star(ig,l)=0.
1411 END IF
1412
1413 ! print*,'ok detr_star'
1414 END IF
1415 ! prise en compte du detrainement dans le calcul du flux
1416 f_star(ig, l+1) = f_star(ig, l) + alim_star(ig, l) + &
1417 ✗ entr_star(ig, l) - detr_star(ig, l)
1418 ! test
1419 ! if (f_star(ig,l+1).lt.0.) then
1420 ! f_star(ig,l+1)=0.
1421 ! entr_star(ig,l)=0.
1422 ! detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)
1423 ! endif
1424 ! test sur le signe de f_star
1425 ✗ IF (f_star(ig,l+1)>1.E-10) THEN
1426 ! then
1427 ! test
1428 ! if (((f_star(ig,l+1)+detr_star(ig,l)).gt.1.e-10)) then
1429 ! AM on melange Tl et qt du thermique
1430 ! on rajoute un exc�s de T dans la couche alim
1431 ! if (l.lt.lentr(ig)) then
1432 ! ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+
1433 ! s
1434 ! (alim_star(ig,l)+entr_star(ig,l))*(zthl(ig,l)+dtheta(ig,l)))
1435 ! s /(f_star(ig,l+1)+detr_star(ig,l))
1436 ! else
1437 ztla(ig, l) = (f_star(ig,l)*ztla(ig,l-1)+(alim_star(ig, &
1438 ✗ l)+entr_star(ig,l))*zthl(ig,l))/(f_star(ig,l+1)+detr_star(ig,l))
1439 ! s /(f_star(ig,l+1))
1440 ! endif
1441 ! on rajoute un exc�s de q dans la couche alim
1442 ! if (l.lt.lentr(ig)) then
1443 ! zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+
1444 ! s (alim_star(ig,l)+entr_star(ig,l))*(po(ig,l)+0.001))
1445 ! s /(f_star(ig,l+1)+detr_star(ig,l))
1446 ! else
1447 zqta(ig, l) = (f_star(ig,l)*zqta(ig,l-1)+(alim_star(ig, &
1448 ✗ l)+entr_star(ig,l))*po(ig,l))/(f_star(ig,l+1)+detr_star(ig,l))
1449 ! s /(f_star(ig,l+1))
1450 ! endif
1451 ! AM on en deduit thetav et ql du thermique
1452 ! CR test
1453 ! Tbef(ig)=ztla(ig,l)*zpspsk(ig,l)
1454 ✗ tbef(ig) = ztla(ig, l)*zpspsk(ig, l)
1455 ✗ zdelta = max(0., sign(1.,rtt-tbef(ig)))
1456 ✗ qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, l)
1457 ✗ qsatbef(ig) = min(0.5, qsatbef(ig))
1458 ✗ zcor = 1./(1.-retv*qsatbef(ig))
1459 ✗ qsatbef(ig) = qsatbef(ig)*zcor
1460 ✗ zsat(ig) = (max(0.,zqta(ig,l)-qsatbef(ig))>1.E-10)
1461
1462 ✗ IF (zsat(ig) .AND. (1==1)) THEN
1463 qlbef = max(0., zqta(ig,l)-qsatbef(ig))
1464 ✗ dt = 0.5*rlvcp*qlbef
1465 ! write(17,*)'DT0=',DT
1466 ✗ DO WHILE (abs(dt)>ddt0)
1467 ! print*,'aie'
1468 ✗ tbef(ig) = tbef(ig) + dt
1469 ✗ zdelta = max(0., sign(1.,rtt-tbef(ig)))
1470 ✗ qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, l)
1471 ✗ qsatbef(ig) = min(0.5, qsatbef(ig))
1472 ✗ zcor = 1./(1.-retv*qsatbef(ig))
1473 ✗ qsatbef(ig) = qsatbef(ig)*zcor
1474 ✗ qlbef = zqta(ig, l) - qsatbef(ig)
1475
1476 zdelta = max(0., sign(1.,rtt-tbef(ig)))
1477 ✗ zcvm5 = r5les*(1.-zdelta) + r5ies*zdelta
1478 ✗ zcor = 1./(1.-retv*qsatbef(ig))
1479 ✗ dqsat_dt = foede(tbef(ig), zdelta, zcvm5, qsatbef(ig), zcor)
1480 ✗ num = -tbef(ig) + ztla(ig, l)*zpspsk(ig, l) + rlvcp*qlbef
1481 ✗ denom = 1. + rlvcp*dqsat_dt
1482 ✗ IF (denom<1.E-10) THEN
1483 ✗ PRINT *, 'pb denom'
1484 END IF
1485 ✗ dt = num/denom
1486 ! write(17,*)'DT=',DT
1487 END DO
1488 zqla(ig, l) = max(0., zqta(ig,l)-qsatbef(ig))
1489 ✗ zqla(ig, l) = max(0., qlbef)
1490 ! zqla(ig,l)=0.
1491 END IF
1492 ! zqla(ig,l) = max(0.,zqta(ig,l)-qsatbef(ig))
1493
1494 ! on ecrit de maniere conservative (sat ou non)
1495 ! T = Tl +Lv/Cp ql
1496 ! CR rq utilisation de humidite specifique ou rapport de melange?
1497 ✗ ztva(ig, l) = ztla(ig, l)*zpspsk(ig, l) + rlvcp*zqla(ig, l)
1498 ✗ ztva(ig, l) = ztva(ig, l)/zpspsk(ig, l)
1499 ! on rajoute le calcul de zha pour diagnostiques (temp potentielle)
1500 ✗ zha(ig, l) = ztva(ig, l)
1501 ! if (l.lt.lentr(ig)) then
1502 ! ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)
1503 ! s -zqla(ig,l))-zqla(ig,l)) + 0.1
1504 ! else
1505 ztva(ig, l) = ztva(ig, l)*(1.+retv*(zqta(ig,l)-zqla(ig, &
1506 ✗ l))-zqla(ig,l))
1507 ! endif
1508 ! ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
1509 ! s /(1.-retv*zqla(ig,l))
1510 ! ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1511 ! ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)
1512 ! s /(1.-retv*zqta(ig,l))
1513 ! s -zqla(ig,l)/(1.-retv*zqla(ig,l)))
1514 ! s -zqla(ig,l)/(1.-retv*zqla(ig,l)))
1515 ! write(13,*)zqla(ig,l),zqla(ig,l)/(1.-retv*zqla(ig,l))
1516 ! on ecrit zqsat
1517 ✗ zqsatth(ig, l) = qsatbef(ig)
1518 ! enddo
1519 ! DO ig=1,ngrid
1520 ! if (zw2(ig,l).ge.1.e-10.and.
1521 ! s f_star(ig,l)+entr_star(ig,l).gt.1.e-10) then
1522 ! mise a jour de la vitesse ascendante (l'air entraine de la couche
1523 ! consideree commence avec une vitesse nulle).
1524
1525 ! if (f_star(ig,l+1).gt.1.e-10) then
1526 zw2(ig, l+1) = zw2(ig, l)* & ! s
1527 ! ((f_star(ig,l)-detr_star(ig,l))**2)
1528 ! s /f_star(ig,l+1)**2+
1529 ((f_star(ig,l))**2)/(f_star(ig,l+1)+detr_star(ig,l))**2 + & ! s
1530 ! /(f_star(ig,l+1))**2+
1531 ✗ 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
1532 ! s *(f_star(ig,l)/f_star(ig,l+1))**2
1533
1534 END IF
1535 END IF
1536
1537 ✗ IF (zw2(ig,l+1)<0.) THEN
1538 linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
1539 ✗ ig,l))
1540 ✗ zw2(ig, l+1) = 0.
1541 ! print*,'linter=',linter(ig)
1542 ! else if ((zw2(ig,l+1).lt.1.e-10).and.(zw2(ig,l+1).ge.0.)) then
1543 ! linter(ig)=l+1
1544 ! print*,'linter=l',zw2(ig,l),zw2(ig,l+1)
1545 ELSE
1546 ✗ wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
1547 ! wa_moy(ig,l+1)=zw2(ig,l+1)
1548 END IF
1549 ✗ IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
1550 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
1551 ✗ lmix(ig) = l + 1
1552 ✗ wmaxa(ig) = wa_moy(ig, l+1)
1553 END IF
1554 END DO
1555 END DO
1556 ✗ PRINT *, 'fin calcul zw2'
1557
1558 ! Calcul de la couche correspondant a la hauteur du thermique
1559 ✗ DO ig = 1, ngrid
1560 ✗ lmax(ig) = lentr(ig)
1561 END DO
1562 ✗ DO ig = 1, ngrid
1563 ✗ DO l = nlay, lentr(ig) + 1, -1
1564 ✗ IF (zw2(ig,l)<=1.E-10) THEN
1565 ✗ lmax(ig) = l - 1
1566 END IF
1567 END DO
1568 END DO
1569 ! pas de thermique si couche 1 stable
1570 ✗ DO ig = 1, ngrid
1571 ✗ IF (lmin(ig)>1) THEN
1572 ✗ lmax(ig) = 1
1573 ✗ lmin(ig) = 1
1574 ✗ lentr(ig) = 1
1575 END IF
1576 END DO
1577
1578 ! Determination de zw2 max
1579 ✗ DO ig = 1, ngrid
1580 ✗ wmax(ig) = 0.
1581 END DO
1582
1583 ✗ DO l = 1, nlay
1584 ✗ DO ig = 1, ngrid
1585 ✗ IF (l<=lmax(ig)) THEN
1586 ✗ IF (zw2(ig,l)<0.) THEN
1587 ✗ PRINT *, 'pb2 zw2<0'
1588 END IF
1589 ✗ zw2(ig, l) = sqrt(zw2(ig,l))
1590 ✗ wmax(ig) = max(wmax(ig), zw2(ig,l))
1591 ELSE
1592 ✗ zw2(ig, l) = 0.
1593 END IF
1594 END DO
1595 END DO
1596
1597 ! Longueur caracteristique correspondant a la hauteur des thermiques.
1598 ✗ DO ig = 1, ngrid
1599 ✗ zmax(ig) = 0.
1600 ✗ zlevinter(ig) = zlev(ig, 1)
1601 END DO
1602 ✗ DO ig = 1, ngrid
1603 ! calcul de zlevinter
1604 zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
1605 ✗ zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
1606 ! pour le cas ou on prend tjs lmin=1
1607 ! zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
1608 ✗ zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,1))
1609 ✗ zmax0(ig) = zmax(ig)
1610 ✗ WRITE (11, *) 'ig,lmax,linter', ig, lmax(ig), linter(ig)
1611 ✗ WRITE (12, *) 'ig,zlevinter,zmax', ig, zmax(ig), zlevinter(ig)
1612 END DO
1613
1614 ! Calcul de zmax_sec et wmax_sec
1615 CALL fermeture_seche(ngrid, nlay, pplay, pplev, pphi, zlev, rhobarz, f0, &
1616 zpspsk, alim, zh, zo, lentr, lmin, nu_min, nu_max, r_aspect, zmax_sec2, &
1617 ✗ wmax_sec2)
1618
1619 ✗ PRINT *, 'avant fermeture'
1620 ! Fermeture,determination de f
1621 ! en lmax f=d-e
1622 ✗ DO ig = 1, ngrid
1623 ! entr_star(ig,lmax(ig))=0.
1624 ! f_star(ig,lmax(ig)+1)=0.
1625 ! detr_star(ig,lmax(ig))=f_star(ig,lmax(ig))+entr_star(ig,lmax(ig))
1626 ! s +alim_star(ig,lmax(ig))
1627 END DO
1628
1629 ✗ DO ig = 1, ngrid
1630 ✗ alim_star2(ig) = 0.
1631 END DO
1632 ! calcul de entr_star_tot
1633 ✗ DO ig = 1, ngrid
1634 ✗ DO k = 1, lmix(ig)
1635 entr_star_tot(ig) = entr_star_tot(ig) & ! s
1636 ! +entr_star(ig,k)
1637 ✗ +alim_star(ig, k)
1638 ! s -detr_star(ig,k)
1639 detr_star_tot(ig) = detr_star_tot(ig) & ! s
1640 ! +alim_star(ig,k)
1641 ✗ -detr_star(ig, k) + entr_star(ig, k)
1642 END DO
1643 END DO
1644
1645 ✗ DO ig = 1, ngrid
1646 ✗ IF (alim_star_tot(ig)<1.E-10) THEN
1647 ✗ f(ig) = 0.
1648 ELSE
1649 ! do k=lmin(ig),lentr(ig)
1650 ✗ DO k = 1, lentr(ig)
1651 alim_star2(ig) = alim_star2(ig) + alim_star(ig, k)**2/(rho(ig,k)*( &
1652 ✗ zlev(ig,k+1)-zlev(ig,k)))
1653 END DO
1654 ✗ IF ((zmax_sec(ig)>1.E-10) .AND. (1==1)) THEN
1655 ✗ f(ig) = wmax_sec(ig)/(max(500.,zmax_sec(ig))*r_aspect*alim_star2(ig))
1656 f(ig) = f(ig) + (f0(ig)-f(ig))*exp((-ptimestep/zmax_sec(ig))*wmax_sec &
1657 ✗ (ig))
1658 ELSE
1659 ✗ f(ig) = wmax(ig)/(max(500.,zmax(ig))*r_aspect*alim_star2(ig))
1660 ✗ f(ig) = f(ig) + (f0(ig)-f(ig))*exp((-ptimestep/zmax(ig))*wmax(ig))
1661 END IF
1662 END IF
1663 ✗ f0(ig) = f(ig)
1664 END DO
1665 ✗ PRINT *, 'apres fermeture'
1666 ! Calcul de l'entrainement
1667 ✗ DO ig = 1, ngrid
1668 ✗ DO k = 1, klev
1669 ✗ alim(ig, k) = f(ig)*alim_star(ig, k)
1670 END DO
1671 END DO
1672 ! CR:test pour entrainer moins que la masse
1673 ! do ig=1,ngrid
1674 ! do l=1,lentr(ig)
1675 ! if ((alim(ig,l)*ptimestep).gt.(0.9*masse(ig,l))) then
1676 ! alim(ig,l+1)=alim(ig,l+1)+alim(ig,l)
1677 ! s -0.9*masse(ig,l)/ptimestep
1678 ! alim(ig,l)=0.9*masse(ig,l)/ptimestep
1679 ! endif
1680 ! enddo
1681 ! enddo
1682 ! calcul du d�trainement
1683 ✗ DO ig = 1, klon
1684 ✗ DO k = 1, klev
1685 ✗ detr(ig, k) = f(ig)*detr_star(ig, k)
1686 ✗ IF (detr(ig,k)<0.) THEN
1687 ! print*,'detr1<0!!!'
1688 END IF
1689 END DO
1690 ✗ DO k = 1, klev
1691 ✗ entr(ig, k) = f(ig)*entr_star(ig, k)
1692 ✗ IF (entr(ig,k)<0.) THEN
1693 ! print*,'entr1<0!!!'
1694 END IF
1695 END DO
1696 END DO
1697
1698 ! do ig=1,ngrid
1699 ! do l=1,klev
1700 ! if (((detr(ig,l)+entr(ig,l)+alim(ig,l))*ptimestep).gt.
1701 ! s (masse(ig,l))) then
1702 ! print*,'d2+e2+a2>m2','ig=',ig,'l=',l,'lmax(ig)=',lmax(ig),'d+e+a='
1703 ! s,(detr(ig,l)+entr(ig,l)+alim(ig,l))*ptimestep,'m=',masse(ig,l)
1704 ! endif
1705 ! enddo
1706 ! enddo
1707 ! Calcul des flux
1708
1709 ✗ DO ig = 1, ngrid
1710 ✗ DO l = 1, lmax(ig)
1711 ! do l=1,klev
1712 ! fmc(ig,l+1)=f(ig)*f_star(ig,l+1)
1713 ✗ fmc(ig, l+1) = fmc(ig, l) + alim(ig, l) + entr(ig, l) - detr(ig, l)
1714 ! print*,'??!!','ig=',ig,'l=',l,'lmax=',lmax(ig),'lmix=',lmix(ig),
1715 ! s 'e=',entr(ig,l),'d=',detr(ig,l),'a=',alim(ig,l),'f=',fmc(ig,l),
1716 ! s 'f+1=',fmc(ig,l+1)
1717 ✗ IF (fmc(ig,l+1)<0.) THEN
1718 ✗ PRINT *, 'fmc1<0', l + 1, lmax(ig), fmc(ig, l+1)
1719 ✗ fmc(ig, l+1) = fmc(ig, l)
1720 ✗ detr(ig, l) = alim(ig, l) + entr(ig, l)
1721 ! fmc(ig,l+1)=0.
1722 ! print*,'fmc1<0',l+1,lmax(ig),fmc(ig,l+1)
1723 END IF
1724 ! if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1725 ! f_old=fmc(ig,l+1)
1726 ! fmc(ig,l+1)=fmc(ig,l)
1727 ! detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l+1)
1728 ! endif
1729
1730 ! if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1731 ! f_old=fmc(ig,l+1)
1732 ! fmc(ig,l+1)=fmc(ig,l)
1733 ! detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l)
1734 ! endif
1735 ! rajout du test sur alpha croissant
1736 ! if test
1737 ! if (1.eq.0) then
1738
1739 ✗ IF (l==klev) THEN
1740 ✗ PRINT *, 'THERMCELL PB ig=', ig, ' l=', l
1741 ✗ abort_message = 'THERMCELL PB'
1742 ✗ CALL abort_physic(modname, abort_message, 1)
1743 END IF
1744 ! if ((zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10).and.
1745 ! s (l.ge.lentr(ig)).and.
1746 ✗ IF ((zw2(ig,l+1)>1.E-10) .AND. (zw2(ig,l)>1.E-10) .AND. (l>=lentr(ig))) &
1747 THEN
1748 ✗ IF (((fmc(ig,l+1)/(rhobarz(ig,l+1)*zw2(ig,l+1)))>(fmc(ig,l)/ &
1749 (rhobarz(ig,l)*zw2(ig,l))))) THEN
1750 f_old = fmc(ig, l+1)
1751 fmc(ig, l+1) = fmc(ig, l)*rhobarz(ig, l+1)*zw2(ig, l+1)/ &
1752 ✗ (rhobarz(ig,l)*zw2(ig,l))
1753 ✗ detr(ig, l) = detr(ig, l) + f_old - fmc(ig, l+1)
1754 ! detr(ig,l)=(fmc(ig,l+1)-fmc(ig,l))/(0.4-1.)
1755 ! entr(ig,l)=0.4*detr(ig,l)
1756 ! entr(ig,l)=fmc(ig,l+1)-fmc(ig,l)+detr(ig,l)
1757 END IF
1758 END IF
1759 ✗ IF ((fmc(ig,l+1)>fmc(ig,l)) .AND. (l>lentr(ig))) THEN
1760 f_old = fmc(ig, l+1)
1761 ✗ fmc(ig, l+1) = fmc(ig, l)
1762 ✗ detr(ig, l) = detr(ig, l) + f_old - fmc(ig, l+1)
1763 END IF
1764 ✗ IF (detr(ig,l)>fmc(ig,l)) THEN
1765 ✗ detr(ig, l) = fmc(ig, l)
1766 ✗ entr(ig, l) = fmc(ig, l+1) - alim(ig, l)
1767 END IF
1768 ✗ IF (fmc(ig,l+1)<0.) THEN
1769 ✗ detr(ig, l) = detr(ig, l) + fmc(ig, l+1)
1770 ✗ fmc(ig, l+1) = 0.
1771 ✗ PRINT *, 'fmc2<0', l + 1, lmax(ig)
1772 END IF
1773
1774 ! test pour ne pas avoir f=0 et d=e/=0
1775 ! if (fmc(ig,l+1).lt.1.e-10) then
1776 ! detr(ig,l+1)=0.
1777 ! entr(ig,l+1)=0.
1778 ! zqla(ig,l+1)=0.
1779 ! zw2(ig,l+1)=0.
1780 ! lmax(ig)=l+1
1781 ! zmax(ig)=zlev(ig,lmax(ig))
1782 ! endif
1783 ✗ IF (zw2(ig,l+1)>1.E-10) THEN
1784 ✗ IF ((((fmc(ig,l+1))/(rhobarz(ig,l+1)*zw2(ig,l+1)))>1.)) THEN
1785 f_old = fmc(ig, l+1)
1786 ✗ fmc(ig, l+1) = rhobarz(ig, l+1)*zw2(ig, l+1)
1787 ✗ zw2(ig, l+1) = 0.
1788 ✗ zqla(ig, l+1) = 0.
1789 ✗ detr(ig, l) = detr(ig, l) + f_old - fmc(ig, l+1)
1790 ✗ lmax(ig) = l + 1
1791 ✗ zmax(ig) = zlev(ig, lmax(ig))
1792 ✗ PRINT *, 'alpha>1', l + 1, lmax(ig)
1793 END IF
1794 END IF
1795 ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
1796 ! endif test
1797 ! endif
1798 END DO
1799 END DO
1800 ✗ DO ig = 1, ngrid
1801 ! if (fmc(ig,lmax(ig)+1).ne.0.) then
1802 ✗ fmc(ig, lmax(ig)+1) = 0.
1803 ✗ entr(ig, lmax(ig)) = 0.
1804 detr(ig, lmax(ig)) = fmc(ig, lmax(ig)) + entr(ig, lmax(ig)) + &
1805 ✗ alim(ig, lmax(ig))
1806 ! endif
1807 END DO
1808 ! test sur le signe de fmc
1809 ✗ DO ig = 1, ngrid
1810 ✗ DO l = 1, klev + 1
1811 ✗ IF (fmc(ig,l)<0.) THEN
1812 ✗ PRINT *, 'fm1<0!!!', 'ig=', ig, 'l=', l, 'a=', alim(ig, l-1), 'e=', &
1813 ✗ entr(ig, l-1), 'f=', fmc(ig, l-1), 'd=', detr(ig, l-1), 'f+1=', &
1814 ✗ fmc(ig, l)
1815 END IF
1816 END DO
1817 END DO
1818 ! test de verification
1819 ✗ DO ig = 1, ngrid
1820 ✗ DO l = 1, lmax(ig)
1821 IF ((abs(fmc(ig,l+1)-fmc(ig,l)-alim(ig,l)-entr(ig,l)+ &
1822 detr(ig,l)))>1.E-4) THEN
1823 ! print*,'pbcm!!','ig=',ig,'l=',l,'lmax=',lmax(ig),'lmix=',lmix(ig),
1824 ! s 'e=',entr(ig,l),'d=',detr(ig,l),'a=',alim(ig,l),'f=',fmc(ig,l),
1825 ! s 'f+1=',fmc(ig,l+1)
1826 END IF
1827 ✗ IF (detr(ig,l)<0.) THEN
1828 ✗ PRINT *, 'detrdemi<0!!!'
1829 END IF
1830 END DO
1831 END DO
1832
1833 ! RC
1834 ! CR def de zmix continu (profil parabolique des vitesses)
1835 ✗ DO ig = 1, ngrid
1836 ✗ IF (lmix(ig)>1.) THEN
1837 ! test
1838 ✗ IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
1839 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
1840 zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- &
1841 (zlev(ig,lmix(ig)))))>1E-10) THEN
1842
1843 zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) &
1844 )**2-(zlev(ig,lmix(ig)+1))**2)-(zw2(ig,lmix(ig))-zw2(ig, &
1845 lmix(ig)+1))*((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))/ &
1846 (2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
1847 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
1848 ✗ zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
1849 ELSE
1850 ✗ zmix(ig) = zlev(ig, lmix(ig))
1851 ✗ PRINT *, 'pb zmix'
1852 END IF
1853 ELSE
1854 ✗ zmix(ig) = 0.
1855 END IF
1856 ! test
1857 ✗ IF ((zmax(ig)-zmix(ig))<=0.) THEN
1858 ✗ zmix(ig) = 0.9*zmax(ig)
1859 ! print*,'pb zmix>zmax'
1860 END IF
1861 END DO
1862 ✗ DO ig = 1, klon
1863 ✗ zmix0(ig) = zmix(ig)
1864 END DO
1865
1866 ! calcul du nouveau lmix correspondant
1867 ✗ DO ig = 1, ngrid
1868 ✗ DO l = 1, klev
1869 ✗ IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
1870 ✗ lmix(ig) = l
1871 END IF
1872 END DO
1873 END DO
1874
1875 ! ne devrait pas arriver!!!!!
1876 ✗ DO ig = 1, ngrid
1877 ✗ DO l = 1, klev
1878 ✗ IF (detr(ig,l)>(fmc(ig,l)+alim(ig,l))+entr(ig,l)) THEN
1879 ✗ PRINT *, 'detr2>fmc2!!!', 'ig=', ig, 'l=', l, 'd=', detr(ig, l), &
1880 ✗ 'f=', fmc(ig, l), 'lmax=', lmax(ig)
1881 ! detr(ig,l)=fmc(ig,l)+alim(ig,l)+entr(ig,l)
1882 ! entr(ig,l)=0.
1883 ! fmc(ig,l+1)=0.
1884 ! zw2(ig,l+1)=0.
1885 ! zqla(ig,l+1)=0.
1886 ✗ PRINT *, 'pb!fm=0 et f_star>0', l, lmax(ig)
1887 ! lmax(ig)=l
1888 END IF
1889 END DO
1890 END DO
1891 ✗ DO ig = 1, ngrid
1892 ✗ DO l = lmax(ig) + 1, klev + 1
1893 ! fmc(ig,l)=0.
1894 ! detr(ig,l)=0.
1895 ! entr(ig,l)=0.
1896 ! zw2(ig,l)=0.
1897 ! zqla(ig,l)=0.
1898 END DO
1899 END DO
1900
1901 ! Calcul du detrainement lors du premier passage
1902 ! print*,'9 OK convect8'
1903 ! print*,'WA1 ',wa_moy
1904
1905 ! determination de l'indice du debut de la mixed layer ou w decroit
1906
1907 ! calcul de la largeur de chaque ascendance dans le cas conservatif.
1908 ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
1909 ! d'une couche est �gale � la hauteur de la couche alimentante.
1910 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
1911 ! de la vitesse d'entrainement horizontal dans la couche alimentante.
1912
1913 ✗ DO l = 2, nlay
1914 ✗ DO ig = 1, ngrid
1915 ✗ IF (l<=lmax(ig) .AND. (test(ig)==1)) THEN
1916 ✗ zw = max(wa_moy(ig,l), 1.E-10)
1917 ✗ larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw)
1918 END IF
1919 END DO
1920 END DO
1921
1922 ✗ DO l = 2, nlay
1923 ✗ DO ig = 1, ngrid
1924 ✗ IF (l<=lmax(ig) .AND. (test(ig)==1)) THEN
1925 ! if (idetr.eq.0) then
1926 ! cette option est finalement en dur.
1927 ✗ IF ((l_mix*zlev(ig,l))<0.) THEN
1928 ✗ PRINT *, 'pb l_mix*zlev<0'
1929 END IF
1930 ! CR: test: nouvelle def de lambda
1931 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
1932 ✗ IF (zw2(ig,l)>1.E-10) THEN
1933 ✗ larg_detr(ig, l) = sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
1934 ELSE
1935 ✗ larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
1936 END IF
1937 ! else if (idetr.eq.1) then
1938 ! larg_detr(ig,l)=larg_cons(ig,l)
1939 ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
1940 ! else if (idetr.eq.2) then
1941 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
1942 ! s *sqrt(wa_moy(ig,l))
1943 ! else if (idetr.eq.4) then
1944 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
1945 ! s *wa_moy(ig,l)
1946 ! endif
1947 END IF
1948 END DO
1949 END DO
1950
1951 ! print*,'10 OK convect8'
1952 ! print*,'WA2 ',wa_moy
1953 ! cal1cul de la fraction de la maille concern�e par l'ascendance en tenant
1954 ! compte de l'epluchage du thermique.
1955
1956
1957 ✗ DO l = 2, nlay
1958 ✗ DO ig = 1, ngrid
1959 ✗ IF (larg_cons(ig,l)>1. .AND. (test(ig)==1)) THEN
1960 ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
1961 ✗ fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
1962 ! test
1963 ✗ fraca(ig, l) = max(fraca(ig,l), 0.)
1964 ✗ fraca(ig, l) = min(fraca(ig,l), 0.5)
1965 ✗ fracd(ig, l) = 1. - fraca(ig, l)
1966 ✗ fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
1967 ELSE
1968 ! wa_moy(ig,l)=0.
1969 ✗ fraca(ig, l) = 0.
1970 ✗ fracc(ig, l) = 0.
1971 ✗ fracd(ig, l) = 1.
1972 END IF
1973 END DO
1974 END DO
1975 ! CR: calcul de fracazmix
1976 ✗ DO ig = 1, ngrid
1977 ✗ IF (test(ig)==1) THEN
1978 fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
1979 (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
1980 fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca( &
1981 ✗ ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
1982 END IF
1983 END DO
1984
1985 ✗ DO l = 2, nlay
1986 ✗ DO ig = 1, ngrid
1987 ✗ IF (larg_cons(ig,l)>1. .AND. (test(ig)==1)) THEN
1988 ✗ IF (l>lmix(ig)) THEN
1989 ! test
1990 ✗ IF (zmax(ig)-zmix(ig)<1.E-10) THEN
1991 ! print*,'pb xxx'
1992 ✗ xxx(ig, l) = (lmax(ig)+1.-l)/(lmax(ig)+1.-lmix(ig))
1993 ELSE
1994 ✗ xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
1995 END IF
1996 ✗ IF (idetr==0) THEN
1997 ✗ fraca(ig, l) = fracazmix(ig)
1998 ✗ ELSE IF (idetr==1) THEN
1999 ✗ fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
2000 ✗ ELSE IF (idetr==2) THEN
2001 ✗ fraca(ig, l) = fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
2002 ELSE
2003 ✗ fraca(ig, l) = fracazmix(ig)*xxx(ig, l)**2
2004 END IF
2005 ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
2006 ✗ fraca(ig, l) = max(fraca(ig,l), 0.)
2007 ✗ fraca(ig, l) = min(fraca(ig,l), 0.5)
2008 ✗ fracd(ig, l) = 1. - fraca(ig, l)
2009 ✗ fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
2010 END IF
2011 END IF
2012 END DO
2013 END DO
2014
2015 ✗ PRINT *, 'fin calcul fraca'
2016 ! print*,'11 OK convect8'
2017 ! print*,'Ea3 ',wa_moy
2018 ! ------------------------------------------------------------------
2019 ! Calcul de fracd, wd
2020 ! somme wa - wd = 0
2021 ! ------------------------------------------------------------------
2022
2023
2024 ✗ DO ig = 1, ngrid
2025 ✗ fm(ig, 1) = 0.
2026 ✗ fm(ig, nlay+1) = 0.
2027 END DO
2028
2029 ✗ DO l = 2, nlay
2030 ✗ DO ig = 1, ngrid
2031 ✗ IF (test(ig)==1) THEN
2032 ✗ fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
2033 ! CR:test
2034 ✗ IF (alim(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) &
2035 THEN
2036 ✗ fm(ig, l) = fm(ig, l-1)
2037 ! write(1,*)'ajustement fm, l',l
2038 END IF
2039 ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
2040 ! RC
2041 END IF
2042 END DO
2043 ✗ DO ig = 1, ngrid
2044 ✗ IF (fracd(ig,l)<0.1 .AND. (test(ig)==1)) THEN
2045 ✗ abort_message = 'fracd trop petit'
2046 ✗ CALL abort_physic(modname, abort_message, 1)
2047 ELSE
2048 ! vitesse descendante "diagnostique"
2049 ✗ wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
2050 END IF
2051 END DO
2052 END DO
2053
2054 ✗ DO l = 1, nlay + 1
2055 ✗ DO ig = 1, ngrid
2056 ✗ IF (test(ig)==0) THEN
2057 ✗ fm(ig, l) = fmc(ig, l)
2058 END IF
2059 END DO
2060 END DO
2061
2062 ! fin du first
2063 ✗ DO l = 1, nlay
2064 ✗ DO ig = 1, ngrid
2065 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
2066 ✗ masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
2067 END DO
2068 END DO
2069
2070 ! print*,'12 OK convect8'
2071 ! print*,'WA4 ',wa_moy
2072 ! c------------------------------------------------------------------
2073 ! calcul du transport vertical
2074 ! ------------------------------------------------------------------
2075
2076 GO TO 4444
2077 ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
2078 DO l = 2, nlay - 1
2079 DO ig = 1, ngrid
2080 IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( &
2081 ig,l+1)) THEN
2082 PRINT *, 'WARN!!! FM>M ig=', ig, ' l=', l, ' FM=', &
2083 fm(ig, l+1)*ptimestep, ' M=', masse(ig, l), masse(ig, l+1)
2084 END IF
2085 END DO
2086 END DO
2087
2088 DO l = 1, nlay
2089 DO ig = 1, ngrid
2090 IF ((alim(ig,l)+entr(ig,l))*ptimestep>masse(ig,l)) THEN
2091 PRINT *, 'WARN!!! E>M ig=', ig, ' l=', l, ' E==', &
2092 (entr(ig,l)+alim(ig,l))*ptimestep, ' M=', masse(ig, l)
2093 END IF
2094 END DO
2095 END DO
2096
2097 DO l = 1, nlay
2098 DO ig = 1, ngrid
2099 IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
2100 ! print*,'WARN!!! fm exagere ig=',ig,' l=',l
2101 ! s ,' FM=',fm(ig,l)
2102 END IF
2103 IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
2104 ! print*,'WARN!!! masse exagere ig=',ig,' l=',l
2105 ! s ,' M=',masse(ig,l)
2106 ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
2107 ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
2108 ! print*,'zlev(ig,l+1),zlev(ig,l)'
2109 ! s ,zlev(ig,l+1),zlev(ig,l)
2110 ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
2111 ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
2112 END IF
2113 IF (.NOT. alim(ig,l)>=0. .OR. .NOT. alim(ig,l)<=10.) THEN
2114 ! print*,'WARN!!! entr exagere ig=',ig,' l=',l
2115 ! s ,' E=',entr(ig,l)
2116 END IF
2117 END DO
2118 END DO
2119
2120 4444 CONTINUE
2121
2122 ! CR:redefinition du entr
2123 ! CR:test:on ne change pas la def du entr mais la def du fm
2124 ✗ DO l = 1, nlay
2125 ✗ DO ig = 1, ngrid
2126 ✗ IF (test(ig)==1) THEN
2127 ✗ detr(ig, l) = fm(ig, l) + alim(ig, l) - fm(ig, l+1)
2128 ✗ IF (detr(ig,l)<0.) THEN
2129 ! entr(ig,l)=entr(ig,l)-detr(ig,l)
2130 ✗ fm(ig, l+1) = fm(ig, l) + alim(ig, l)
2131 ✗ detr(ig, l) = 0.
2132 ! write(11,*)'l,ig,entr',l,ig,entr(ig,l)
2133 ! print*,'WARNING !!! detrainement negatif ',ig,l
2134 END IF
2135 END IF
2136 END DO
2137 END DO
2138 ! RC
2139
2140 ✗ IF (w2di==1) THEN
2141 ✗ fm0 = fm0 + ptimestep*(fm-fm0)/tho
2142 ✗ entr0 = entr0 + ptimestep*(alim+entr-entr0)/tho
2143 ELSE
2144 ✗ fm0 = fm
2145 ✗ entr0 = alim + entr
2146 ✗ detr0 = detr
2147 ✗ alim0 = alim
2148 ! zoa=zqta
2149 ! entr0=alim
2150 END IF
2151
2152 IF (1==1) THEN
2153 ! call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2154 ! . ,zh,zdhadj,zha)
2155 ! call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2156 ! . ,zo,pdoadj,zoa)
2157 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zthl, &
2158 ✗ zdthladj, zta)
2159 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, po, pdoadj, &
2160 ✗ zoa)
2161 ELSE
2162 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
2163 zdhadj, zha)
2164 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
2165 pdoadj, zoa)
2166 END IF
2167
2168 IF (1==0) THEN
2169 CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
2170 zu, zv, pduadj, pdvadj, zua, zva)
2171 ELSE
2172 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
2173 ✗ zua)
2174 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
2175 ✗ zva)
2176 END IF
2177
2178 ! Calcul des moments
2179 ! do l=1,nlay
2180 ! do ig=1,ngrid
2181 ! zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
2182 ! zf2=zf/(1.-zf)
2183 ! thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
2184 ! wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
2185 ! enddo
2186 ! enddo
2187
2188
2189
2190
2191
2192
2193 ! print*,'13 OK convect8'
2194 ! print*,'WA5 ',wa_moy
2195 ✗ DO l = 1, nlay
2196 ✗ DO ig = 1, ngrid
2197 ! pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
2198 ✗ pdtadj(ig, l) = zdthladj(ig, l)*zpspsk(ig, l)
2199 END DO
2200 END DO
2201
2202
2203 ! do l=1,nlay
2204 ! do ig=1,ngrid
2205 ! if(abs(pdtadj(ig,l))*86400..gt.500.) then
2206 ! print*,'WARN!!! ig=',ig,' l=',l
2207 ! s ,' pdtadj=',pdtadj(ig,l)
2208 ! endif
2209 ! if(abs(pdoadj(ig,l))*86400..gt.1.) then
2210 ! print*,'WARN!!! ig=',ig,' l=',l
2211 ! s ,' pdoadj=',pdoadj(ig,l)
2212 ! endif
2213 ! enddo
2214 ! enddo
2215
2216 ! print*,'14 OK convect8'
2217 ! ------------------------------------------------------------------
2218 ! Calculs pour les sorties
2219 ! ------------------------------------------------------------------
2220 ! calcul de fraca pour les sorties
2221 ✗ DO l = 2, klev
2222 ✗ DO ig = 1, klon
2223 ✗ IF (zw2(ig,l)>1.E-10) THEN
2224 ✗ fraca(ig, l) = fm(ig, l)/(rhobarz(ig,l)*zw2(ig,l))
2225 ELSE
2226 ✗ fraca(ig, l) = 0.
2227 END IF
2228 END DO
2229 END DO
2230 IF (sorties) THEN
2231 DO l = 1, nlay
2232 DO ig = 1, ngrid
2233 zla(ig, l) = (1.-fracd(ig,l))*zmax(ig)
2234 zld(ig, l) = fracd(ig, l)*zmax(ig)
2235 IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ &
2236 (1.-fracd(ig,l))
2237 END DO
2238 END DO
2239 ! CR calcul du niveau de condensation
2240 ! initialisation
2241 DO ig = 1, ngrid
2242 nivcon(ig) = 0.
2243 zcon(ig) = 0.
2244 END DO
2245 DO k = nlay, 1, -1
2246 DO ig = 1, ngrid
2247 IF (zqla(ig,k)>1E-10) THEN
2248 nivcon(ig) = k
2249 zcon(ig) = zlev(ig, k)
2250 END IF
2251 ! if (zcon(ig).gt.1.e-10) then
2252 ! nuage=.true.
2253 ! else
2254 ! nuage=.false.
2255 ! endif
2256 END DO
2257 END DO
2258
2259 DO l = 1, nlay
2260 DO ig = 1, ngrid
2261 zf = fraca(ig, l)
2262 zf2 = zf/(1.-zf)
2263 thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
2264 wth2(ig, l) = zf2*(zw2(ig,l))**2
2265 ! print*,'wth2=',wth2(ig,l)
2266 wth3(ig, l) = zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))*zw2(ig, l)* &
2267 zw2(ig, l)*zw2(ig, l)
2268 q2(ig, l) = zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
2269 ! test: on calcul q2/po=ratqsc
2270 ! if (nuage) then
2271 ratqscth(ig, l) = sqrt(q2(ig,l))/(po(ig,l)*1000.)
2272 ! else
2273 ! ratqscth(ig,l)=0.
2274 ! endif
2275 END DO
2276 END DO
2277 ! calcul du ratqscdiff
2278 sum = 0.
2279 sumdiff = 0.
2280 ratqsdiff(:, :) = 0.
2281 DO ig = 1, ngrid
2282 DO l = 1, lentr(ig)
2283 sum = sum + alim_star(ig, l)*zqta(ig, l)*1000.
2284 END DO
2285 END DO
2286 DO ig = 1, ngrid
2287 DO l = 1, lentr(ig)
2288 zf = fraca(ig, l)
2289 zf2 = zf/(1.-zf)
2290 sumdiff = sumdiff + alim_star(ig, l)*(zqta(ig,l)*1000.-sum)**2
2291 ! ratqsdiff=ratqsdiff+alim_star(ig,l)*
2292 ! s (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
2293 END DO
2294 END DO
2295 DO l = 1, klev
2296 DO ig = 1, ngrid
2297 ratqsdiff(ig, l) = sqrt(sumdiff)/(po(ig,l)*1000.)
2298 ! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
2299 END DO
2300 END DO
2301
2302 END IF
2303
2304 ! print*,'19 OK convect8'
2305 ✗ RETURN
2306 END SUBROUTINE thermcell_cld
2307
2308 ✗ SUBROUTINE thermcell_eau(ngrid, nlay, ptimestep, pplay, pplev, pphi, pu, pv, &
2309 ✗ pt, po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0 & ! s
2310 ! ,pu_therm,pv_therm
2311 , r_aspect, l_mix, w2di, tho)
2312
2313 USE dimphy
2314 IMPLICIT NONE
2315
2316 ! =======================================================================
2317
2318 ! Calcul du transport verticale dans la couche limite en presence
2319 ! de "thermiques" explicitement representes
2320
2321 ! R��criture � partir d'un listing papier � Habas, le 14/02/00
2322
2323 ! le thermique est suppos� homog�ne et dissip� par m�lange avec
2324 ! son environnement. la longueur l_mix contr�le l'efficacit� du
2325 ! m�lange
2326
2327 ! Le calcul du transport des diff�rentes esp�ces se fait en prenant
2328 ! en compte:
2329 ! 1. un flux de masse montant
2330 ! 2. un flux de masse descendant
2331 ! 3. un entrainement
2332 ! 4. un detrainement
2333
2334 ! =======================================================================
2335
2336 ! -----------------------------------------------------------------------
2337 ! declarations:
2338 ! -------------
2339
2340 include "YOMCST.h"
2341 include "YOETHF.h"
2342 include "FCTTRE.h"
2343
2344 ! arguments:
2345 ! ----------
2346
2347 INTEGER ngrid, nlay, w2di
2348 REAL tho
2349 REAL ptimestep, l_mix, r_aspect
2350 REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
2351 REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
2352 REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
2353 REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
2354 REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
2355 REAL pphi(ngrid, nlay)
2356
2357 INTEGER idetr
2358 SAVE idetr
2359 DATA idetr/3/
2360 !$OMP THREADPRIVATE(idetr)
2361
2362 ! local:
2363 ! ------
2364
2365 ✗ INTEGER ig, k, l, lmaxa(klon), lmix(klon)
2366 REAL zsortie1d(klon)
2367 ! CR: on remplace lmax(klon,klev+1)
2368 ✗ INTEGER lmax(klon), lmin(klon), lentr(klon)
2369 ✗ REAL linter(klon)
2370 ✗ REAL zmix(klon), fracazmix(klon)
2371 ! RC
2372 ✗ REAL zmax(klon), zw, zz, zw2(klon, klev+1), ztva(klon, klev), zzz
2373
2374 ✗ REAL zlev(klon, klev+1), zlay(klon, klev)
2375 ✗ REAL zh(klon, klev), zdhadj(klon, klev)
2376 ✗ REAL zthl(klon, klev), zdthladj(klon, klev)
2377 ✗ REAL ztv(klon, klev)
2378 ✗ REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
2379 ✗ REAL zl(klon, klev)
2380 REAL wh(klon, klev+1)
2381 REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
2382 REAL zla(klon, klev+1)
2383 REAL zwa(klon, klev+1)
2384 REAL zld(klon, klev+1)
2385 REAL zwd(klon, klev+1)
2386 REAL zsortie(klon, klev)
2387 ✗ REAL zva(klon, klev)
2388 ✗ REAL zua(klon, klev)
2389 ✗ REAL zoa(klon, klev)
2390
2391 ✗ REAL zta(klon, klev)
2392 ✗ REAL zha(klon, klev)
2393 ✗ REAL wa_moy(klon, klev+1)
2394 ✗ REAL fraca(klon, klev+1)
2395 ✗ REAL fracc(klon, klev+1)
2396 REAL zf, zf2
2397 ✗ REAL thetath2(klon, klev), wth2(klon, klev)
2398 ! common/comtherm/thetath2,wth2
2399
2400 REAL count_time
2401 INTEGER ialt
2402
2403 LOGICAL sorties
2404 ✗ REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
2405 ✗ REAL zpspsk(klon, klev)
2406
2407 ! real wmax(klon,klev),wmaxa(klon)
2408 ✗ REAL wmax(klon), wmaxa(klon)
2409 ✗ REAL wa(klon, klev, klev+1)
2410 ✗ REAL wd(klon, klev+1)
2411 REAL larg_part(klon, klev, klev+1)
2412 ✗ REAL fracd(klon, klev+1)
2413 ✗ REAL xxx(klon, klev+1)
2414 ✗ REAL larg_cons(klon, klev+1)
2415 ✗ REAL larg_detr(klon, klev+1)
2416 REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
2417 REAL pu_therm(klon, klev), pv_therm(klon, klev)
2418 ✗ REAL fm(klon, klev+1), entr(klon, klev)
2419 ✗ REAL fmc(klon, klev+1)
2420
2421 REAL zcor, zdelta, zcvm5, qlbef
2422 ✗ REAL tbef(klon), qsatbef(klon)
2423 REAL dqsat_dt, dt, num, denom
2424 REAL reps, rlvcp, ddt0
2425 ✗ REAL ztla(klon, klev), zqla(klon, klev), zqta(klon, klev)
2426
2427 PARAMETER (ddt0=.01)
2428
2429 ! CR:nouvelles variables
2430 ✗ REAL f_star(klon, klev+1), entr_star(klon, klev)
2431 ✗ REAL entr_star_tot(klon), entr_star2(klon)
2432 ✗ REAL f(klon), f0(klon)
2433 ✗ REAL zlevinter(klon)
2434 LOGICAL first
2435 DATA first/.FALSE./
2436 SAVE first
2437 !$OMP THREADPRIVATE(first)
2438
2439 ! RC
2440
2441 CHARACTER *2 str2
2442 CHARACTER *10 str10
2443
2444 CHARACTER (LEN=20) :: modname = 'thermcell_eau'
2445 CHARACTER (LEN=80) :: abort_message
2446
2447 LOGICAL vtest(klon), down
2448 ✗ LOGICAL zsat(klon)
2449
2450 EXTERNAL scopy
2451
2452 INTEGER ncorrec, ll
2453 SAVE ncorrec
2454 DATA ncorrec/0/
2455 !$OMP THREADPRIVATE(ncorrec)
2456
2457
2458
2459 ! -----------------------------------------------------------------------
2460 ! initialisation:
2461 ! ---------------
2462
2463 sorties = .TRUE.
2464 ✗ IF (ngrid/=klon) THEN
2465 ✗ PRINT *
2466 ✗ PRINT *, 'STOP dans convadj'
2467 ✗ PRINT *, 'ngrid =', ngrid
2468 ✗ PRINT *, 'klon =', klon
2469 END IF
2470
2471 ! Initialisation
2472 ✗ rlvcp = rlvtt/rcpd
2473 reps = rd/rv
2474
2475 ! -----------------------------------------------------------------------
2476 ! AM Calcul de T,q,ql a partir de Tl et qT
2477 ! ---------------------------------------------------
2478
2479 ! Pr Tprec=Tl calcul de qsat
2480 ! Si qsat>qT T=Tl, q=qT
2481 ! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
2482 ! On cherche DDT < DDT0
2483
2484 ! defaut
2485 ✗ DO ll = 1, nlay
2486 ✗ DO ig = 1, ngrid
2487 ✗ zo(ig, ll) = po(ig, ll)
2488 ✗ zl(ig, ll) = 0.
2489 ✗ zh(ig, ll) = pt(ig, ll)
2490 END DO
2491 END DO
2492 ✗ DO ig = 1, ngrid
2493 ✗ zsat(ig) = .FALSE.
2494 END DO
2495
2496
2497 ✗ DO ll = 1, nlay
2498 ! les points insatures sont definitifs
2499 ✗ DO ig = 1, ngrid
2500 ✗ tbef(ig) = pt(ig, ll)
2501 ✗ zdelta = max(0., sign(1.,rtt-tbef(ig)))
2502 ✗ qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, ll)
2503 ✗ qsatbef(ig) = min(0.5, qsatbef(ig))
2504 ✗ zcor = 1./(1.-retv*qsatbef(ig))
2505 ✗ qsatbef(ig) = qsatbef(ig)*zcor
2506 ✗ zsat(ig) = (max(0.,po(ig,ll)-qsatbef(ig))>0.00001)
2507 END DO
2508
2509 ✗ DO ig = 1, ngrid
2510 ✗ IF (zsat(ig)) THEN
2511 ✗ qlbef = max(0., po(ig,ll)-qsatbef(ig))
2512 ! si sature: ql est surestime, d'ou la sous-relax
2513 ✗ dt = 0.5*rlvcp*qlbef
2514 ! on pourra enchainer 2 ou 3 calculs sans Do while
2515 ✗ DO WHILE (dt>ddt0)
2516 ! il faut verifier si c,a conserve quand on repasse en insature ...
2517 ✗ tbef(ig) = tbef(ig) + dt
2518 ✗ zdelta = max(0., sign(1.,rtt-tbef(ig)))
2519 ✗ qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, ll)
2520 ✗ qsatbef(ig) = min(0.5, qsatbef(ig))
2521 ✗ zcor = 1./(1.-retv*qsatbef(ig))
2522 ✗ qsatbef(ig) = qsatbef(ig)*zcor
2523 ! on veut le signe de qlbef
2524 ✗ qlbef = po(ig, ll) - qsatbef(ig)
2525 ! dqsat_dT
2526 zdelta = max(0., sign(1.,rtt-tbef(ig)))
2527 ✗ zcvm5 = r5les*(1.-zdelta) + r5ies*zdelta
2528 ✗ zcor = 1./(1.-retv*qsatbef(ig))
2529 ✗ dqsat_dt = foede(tbef(ig), zdelta, zcvm5, qsatbef(ig), zcor)
2530 ✗ num = -tbef(ig) + pt(ig, ll) + rlvcp*qlbef
2531 ✗ denom = 1. + rlvcp*dqsat_dt
2532 ✗ dt = num/denom
2533 END DO
2534 ! on ecrit de maniere conservative (sat ou non)
2535 ✗ zl(ig, ll) = max(0., qlbef)
2536 ! T = Tl +Lv/Cp ql
2537 ✗ zh(ig, ll) = pt(ig, ll) + rlvcp*zl(ig, ll)
2538 ✗ zo(ig, ll) = po(ig, ll) - zl(ig, ll)
2539 END IF
2540 END DO
2541 END DO
2542 ! AM fin
2543
2544 ! -----------------------------------------------------------------------
2545 ! incrementation eventuelle de tendances precedentes:
2546 ! ---------------------------------------------------
2547
2548 ! print*,'0 OK convect8'
2549
2550 ✗ DO l = 1, nlay
2551 ✗ DO ig = 1, ngrid
2552 ✗ zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa
2553 ! zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
2554 ✗ zu(ig, l) = pu(ig, l)
2555 ✗ zv(ig, l) = pv(ig, l)
2556 ! zo(ig,l)=po(ig,l)
2557 ! ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
2558 ! AM attention zh est maintenant le profil de T et plus le profil de
2559 ! theta !
2560
2561 ! T-> Theta
2562 ✗ ztv(ig, l) = zh(ig, l)/zpspsk(ig, l)
2563 ! AM Theta_v
2564 ✗ ztv(ig, l) = ztv(ig, l)*(1.+retv*(zo(ig,l))-zl(ig,l))
2565 ! AM Thetal
2566 ✗ zthl(ig, l) = pt(ig, l)/zpspsk(ig, l)
2567
2568 END DO
2569 END DO
2570
2571 ! print*,'1 OK convect8'
2572 ! --------------------
2573
2574
2575 ! + + + + + + + + + + +
2576
2577
2578 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
2579 ! wh,wt,wo ...
2580
2581 ! + + + + + + + + + + + zh,zu,zv,zo,rho
2582
2583
2584 ! -------------------- zlev(1)
2585 ! \\\\\\\\\\\\\\\\\\\
2586
2587
2588 ! -----------------------------------------------------------------------
2589 ! Calcul des altitudes des couches
2590 ! -----------------------------------------------------------------------
2591
2592 ✗ DO l = 2, nlay
2593 ✗ DO ig = 1, ngrid
2594 ✗ zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
2595 END DO
2596 END DO
2597 ✗ DO ig = 1, ngrid
2598 ✗ zlev(ig, 1) = 0.
2599 ✗ zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
2600 END DO
2601 ✗ DO l = 1, nlay
2602 ✗ DO ig = 1, ngrid
2603 ✗ zlay(ig, l) = pphi(ig, l)/rg
2604 END DO
2605 END DO
2606
2607 ! print*,'2 OK convect8'
2608 ! -----------------------------------------------------------------------
2609 ! Calcul des densites
2610 ! -----------------------------------------------------------------------
2611
2612 ✗ DO l = 1, nlay
2613 ✗ DO ig = 1, ngrid
2614 ! rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
2615 ✗ rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*ztv(ig,l))
2616 END DO
2617 END DO
2618
2619 ✗ DO l = 2, nlay
2620 ✗ DO ig = 1, ngrid
2621 ✗ rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
2622 END DO
2623 END DO
2624
2625 ✗ DO k = 1, nlay
2626 ✗ DO l = 1, nlay + 1
2627 ✗ DO ig = 1, ngrid
2628 ✗ wa(ig, k, l) = 0.
2629 END DO
2630 END DO
2631 END DO
2632
2633 ! print*,'3 OK convect8'
2634 ! ------------------------------------------------------------------
2635 ! Calcul de w2, quarre de w a partir de la cape
2636 ! a partir de w2, on calcule wa, vitesse de l'ascendance
2637
2638 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
2639 ! w2 est stoke dans wa
2640
2641 ! ATTENTION: dans convect8, on n'utilise le calcule des wa
2642 ! independants par couches que pour calculer l'entrainement
2643 ! a la base et la hauteur max de l'ascendance.
2644
2645 ! Indicages:
2646 ! l'ascendance provenant du niveau k traverse l'interface l avec
2647 ! une vitesse wa(k,l).
2648
2649 ! --------------------
2650
2651 ! + + + + + + + + + +
2652
2653 ! wa(k,l) ---- -------------------- l
2654 ! / ! /||\ + + + + + + + + + +
2655 ! ||
2656 ! || --------------------
2657 ! ||
2658 ! || + + + + + + + + + +
2659 ! ||
2660 ! || --------------------
2661 ! ||__
2662 ! |___ + + + + + + + + + + k
2663
2664 ! --------------------
2665
2666
2667
2668 ! ------------------------------------------------------------------
2669
2670 ! CR: ponderation entrainement des couches instables
2671 ! def des entr_star tels que entr=f*entr_star
2672 ✗ DO l = 1, klev
2673 ✗ DO ig = 1, ngrid
2674 ✗ entr_star(ig, l) = 0.
2675 END DO
2676 END DO
2677 ! determination de la longueur de la couche d entrainement
2678 ✗ DO ig = 1, ngrid
2679 ✗ lentr(ig) = 1
2680 END DO
2681
2682 ! on ne considere que les premieres couches instables
2683 ✗ DO k = nlay - 1, 1, -1
2684 ✗ DO ig = 1, ngrid
2685 ✗ IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<ztv(ig,k+2)) THEN
2686 ✗ lentr(ig) = k
2687 END IF
2688 END DO
2689 END DO
2690
2691 ! determination du lmin: couche d ou provient le thermique
2692 ✗ DO ig = 1, ngrid
2693 ✗ lmin(ig) = 1
2694 END DO
2695 ✗ DO ig = 1, ngrid
2696 ✗ DO l = nlay, 2, -1
2697 ✗ IF (ztv(ig,l-1)>ztv(ig,l)) THEN
2698 ✗ lmin(ig) = l - 1
2699 END IF
2700 END DO
2701 END DO
2702
2703 ! definition de l'entrainement des couches
2704 ✗ DO l = 1, klev - 1
2705 ✗ DO ig = 1, ngrid
2706 ✗ IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<=lentr(ig)) THEN
2707 ✗ entr_star(ig, l) = (ztv(ig,l)-ztv(ig,l+1))*(zlev(ig,l+1)-zlev(ig,l))
2708 END IF
2709 END DO
2710 END DO
2711 ! pas de thermique si couche 1 stable
2712 ✗ DO ig = 1, ngrid
2713 ✗ IF (lmin(ig)>1) THEN
2714 ✗ DO l = 1, klev
2715 ✗ entr_star(ig, l) = 0.
2716 END DO
2717 END IF
2718 END DO
2719 ! calcul de l entrainement total
2720 ✗ DO ig = 1, ngrid
2721 ✗ entr_star_tot(ig) = 0.
2722 END DO
2723 ✗ DO ig = 1, ngrid
2724 ✗ DO k = 1, klev
2725 ✗ entr_star_tot(ig) = entr_star_tot(ig) + entr_star(ig, k)
2726 END DO
2727 END DO
2728
2729 ✗ DO k = 1, klev
2730 ✗ DO ig = 1, ngrid
2731 ✗ ztva(ig, k) = ztv(ig, k)
2732 END DO
2733 END DO
2734 ! RC
2735 ! AM:initialisations
2736 ✗ DO k = 1, nlay
2737 ✗ DO ig = 1, ngrid
2738 ✗ ztva(ig, k) = ztv(ig, k)
2739 ✗ ztla(ig, k) = zthl(ig, k)
2740 ✗ zqla(ig, k) = 0.
2741 ✗ zqta(ig, k) = po(ig, k)
2742 ✗ zsat(ig) = .FALSE.
2743 END DO
2744 END DO
2745
2746 ! print*,'7 OK convect8'
2747 ✗ DO k = 1, klev + 1
2748 ✗ DO ig = 1, ngrid
2749 ✗ zw2(ig, k) = 0.
2750 ✗ fmc(ig, k) = 0.
2751 ! CR
2752 ✗ f_star(ig, k) = 0.
2753 ! RC
2754 ✗ larg_cons(ig, k) = 0.
2755 ✗ larg_detr(ig, k) = 0.
2756 ✗ wa_moy(ig, k) = 0.
2757 END DO
2758 END DO
2759
2760 ! print*,'8 OK convect8'
2761 ✗ DO ig = 1, ngrid
2762 ✗ linter(ig) = 1.
2763 ✗ lmaxa(ig) = 1
2764 ✗ lmix(ig) = 1
2765 ✗ wmaxa(ig) = 0.
2766 END DO
2767
2768 ! CR:
2769 ✗ DO l = 1, nlay - 2
2770 ✗ DO ig = 1, ngrid
2771 ✗ IF (ztv(ig,l)>ztv(ig,l+1) .AND. entr_star(ig,l)>1.E-10 .AND. &
2772 ✗ zw2(ig,l)<1E-10) THEN
2773 ! AM
2774 ✗ ztla(ig, l) = zthl(ig, l)
2775 ✗ zqta(ig, l) = po(ig, l)
2776 ✗ zqla(ig, l) = zl(ig, l)
2777 ! AM
2778 ✗ f_star(ig, l+1) = entr_star(ig, l)
2779 ! test:calcul de dteta
2780 zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
2781 ✗ (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
2782 ✗ larg_detr(ig, l) = 0.
2783 ✗ ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+entr_star(ig, &
2784 l)>1.E-10)) THEN
2785 ✗ f_star(ig, l+1) = f_star(ig, l) + entr_star(ig, l)
2786
2787 ! AM on melange Tl et qt du thermique
2788 ztla(ig, l) = (f_star(ig,l)*ztla(ig,l-1)+entr_star(ig,l)*zthl(ig,l))/ &
2789 ✗ f_star(ig, l+1)
2790 zqta(ig, l) = (f_star(ig,l)*zqta(ig,l-1)+entr_star(ig,l)*po(ig,l))/ &
2791 ✗ f_star(ig, l+1)
2792
2793 ! ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
2794 ! s *ztv(ig,l))/f_star(ig,l+1)
2795
2796 ! AM on en deduit thetav et ql du thermique
2797 ✗ tbef(ig) = ztla(ig, l)*zpspsk(ig, l)
2798 ✗ zdelta = max(0., sign(1.,rtt-tbef(ig)))
2799 ✗ qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, l)
2800 ✗ qsatbef(ig) = min(0.5, qsatbef(ig))
2801 ✗ zcor = 1./(1.-retv*qsatbef(ig))
2802 ✗ qsatbef(ig) = qsatbef(ig)*zcor
2803 ✗ zsat(ig) = (max(0.,zqta(ig,l)-qsatbef(ig))>0.00001)
2804 END IF
2805 END DO
2806 ✗ DO ig = 1, ngrid
2807 ✗ IF (zsat(ig)) THEN
2808 ✗ qlbef = max(0., zqta(ig,l)-qsatbef(ig))
2809 ✗ dt = 0.5*rlvcp*qlbef
2810 ✗ DO WHILE (dt>ddt0)
2811 ✗ tbef(ig) = tbef(ig) + dt
2812 ✗ zdelta = max(0., sign(1.,rtt-tbef(ig)))
2813 ✗ qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, l)
2814 ✗ qsatbef(ig) = min(0.5, qsatbef(ig))
2815 ✗ zcor = 1./(1.-retv*qsatbef(ig))
2816 ✗ qsatbef(ig) = qsatbef(ig)*zcor
2817 ✗ qlbef = zqta(ig, l) - qsatbef(ig)
2818
2819 zdelta = max(0., sign(1.,rtt-tbef(ig)))
2820 ✗ zcvm5 = r5les*(1.-zdelta) + r5ies*zdelta
2821 ✗ zcor = 1./(1.-retv*qsatbef(ig))
2822 ✗ dqsat_dt = foede(tbef(ig), zdelta, zcvm5, qsatbef(ig), zcor)
2823 ✗ num = -tbef(ig) + ztla(ig, l)*zpspsk(ig, l) + rlvcp*qlbef
2824 ✗ denom = 1. + rlvcp*dqsat_dt
2825 ✗ dt = num/denom
2826 END DO
2827 ✗ zqla(ig, l) = max(0., zqta(ig,l)-qsatbef(ig))
2828 END IF
2829 ! on ecrit de maniere conservative (sat ou non)
2830 ! T = Tl +Lv/Cp ql
2831 ✗ ztva(ig, l) = ztla(ig, l)*zpspsk(ig, l) + rlvcp*zqla(ig, l)
2832 ✗ ztva(ig, l) = ztva(ig, l)/zpspsk(ig, l)
2833 ✗ ztva(ig, l) = ztva(ig, l)*(1.+retv*(zqta(ig,l)-zqla(ig,l))-zqla(ig,l))
2834
2835 END DO
2836 ✗ DO ig = 1, ngrid
2837 ✗ IF (zw2(ig,l)>=1.E-10 .AND. f_star(ig,l)+entr_star(ig,l)>1.E-10) THEN
2838 ! mise a jour de la vitesse ascendante (l'air entraine de la couche
2839 ! consideree commence avec une vitesse nulle).
2840
2841 zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/f_star(ig,l+1))**2 + &
2842 ✗ 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
2843 END IF
2844 ! determination de zmax continu par interpolation lineaire
2845 ✗ IF (zw2(ig,l+1)<0.) THEN
2846 linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
2847 ✗ ig,l))
2848 ✗ zw2(ig, l+1) = 0.
2849 ✗ lmaxa(ig) = l
2850 ELSE
2851 ✗ wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
2852 END IF
2853 ✗ IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
2854 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
2855 ✗ lmix(ig) = l + 1
2856 ✗ wmaxa(ig) = wa_moy(ig, l+1)
2857 END IF
2858 END DO
2859 END DO
2860
2861 ! Calcul de la couche correspondant a la hauteur du thermique
2862 ✗ DO ig = 1, ngrid
2863 ✗ lmax(ig) = lentr(ig)
2864 END DO
2865 ✗ DO ig = 1, ngrid
2866 ✗ DO l = nlay, lentr(ig) + 1, -1
2867 ✗ IF (zw2(ig,l)<=1.E-10) THEN
2868 ✗ lmax(ig) = l - 1
2869 END IF
2870 END DO
2871 END DO
2872 ! pas de thermique si couche 1 stable
2873 ✗ DO ig = 1, ngrid
2874 ✗ IF (lmin(ig)>1) THEN
2875 ✗ lmax(ig) = 1
2876 ✗ lmin(ig) = 1
2877 END IF
2878 END DO
2879
2880 ! Determination de zw2 max
2881 ✗ DO ig = 1, ngrid
2882 ✗ wmax(ig) = 0.
2883 END DO
2884
2885 ✗ DO l = 1, nlay
2886 ✗ DO ig = 1, ngrid
2887 ✗ IF (l<=lmax(ig)) THEN
2888 ✗ zw2(ig, l) = sqrt(zw2(ig,l))
2889 ✗ wmax(ig) = max(wmax(ig), zw2(ig,l))
2890 ELSE
2891 ✗ zw2(ig, l) = 0.
2892 END IF
2893 END DO
2894 END DO
2895
2896 ! Longueur caracteristique correspondant a la hauteur des thermiques.
2897 ✗ DO ig = 1, ngrid
2898 ✗ zmax(ig) = 500.
2899 ✗ zlevinter(ig) = zlev(ig, 1)
2900 END DO
2901 ✗ DO ig = 1, ngrid
2902 ! calcul de zlevinter
2903 zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
2904 ✗ zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
2905 ✗ zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,lmin(ig)))
2906 END DO
2907
2908 ! Fermeture,determination de f
2909 ✗ DO ig = 1, ngrid
2910 ✗ entr_star2(ig) = 0.
2911 END DO
2912 ✗ DO ig = 1, ngrid
2913 ✗ IF (entr_star_tot(ig)<1.E-10) THEN
2914 ✗ f(ig) = 0.
2915 ELSE
2916 ✗ DO k = lmin(ig), lentr(ig)
2917 entr_star2(ig) = entr_star2(ig) + entr_star(ig, k)**2/(rho(ig,k)*( &
2918 ✗ zlev(ig,k+1)-zlev(ig,k)))
2919 END DO
2920 ! Nouvelle fermeture
2921 ✗ f(ig) = wmax(ig)/(zmax(ig)*r_aspect*entr_star2(ig))*entr_star_tot(ig)
2922 ! test
2923 ✗ IF (first) THEN
2924 ✗ f(ig) = f(ig) + (f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)*wmax(ig))
2925 END IF
2926 END IF
2927 ✗ f0(ig) = f(ig)
2928 ✗ first = .TRUE.
2929 END DO
2930
2931 ! Calcul de l'entrainement
2932 ✗ DO k = 1, klev
2933 ✗ DO ig = 1, ngrid
2934 ✗ entr(ig, k) = f(ig)*entr_star(ig, k)
2935 END DO
2936 END DO
2937 ! Calcul des flux
2938 ✗ DO ig = 1, ngrid
2939 ✗ DO l = 1, lmax(ig) - 1
2940 ✗ fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
2941 END DO
2942 END DO
2943
2944 ! RC
2945
2946
2947 ! print*,'9 OK convect8'
2948 ! print*,'WA1 ',wa_moy
2949
2950 ! determination de l'indice du debut de la mixed layer ou w decroit
2951
2952 ! calcul de la largeur de chaque ascendance dans le cas conservatif.
2953 ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
2954 ! d'une couche est �gale � la hauteur de la couche alimentante.
2955 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
2956 ! de la vitesse d'entrainement horizontal dans la couche alimentante.
2957
2958 ✗ DO l = 2, nlay
2959 ✗ DO ig = 1, ngrid
2960 ✗ IF (l<=lmaxa(ig)) THEN
2961 ✗ zw = max(wa_moy(ig,l), 1.E-10)
2962 ✗ larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw)
2963 END IF
2964 END DO
2965 END DO
2966
2967 ✗ DO l = 2, nlay
2968 ✗ DO ig = 1, ngrid
2969 ✗ IF (l<=lmaxa(ig)) THEN
2970 ! if (idetr.eq.0) then
2971 ! cette option est finalement en dur.
2972 ✗ larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
2973 ! else if (idetr.eq.1) then
2974 ! larg_detr(ig,l)=larg_cons(ig,l)
2975 ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
2976 ! else if (idetr.eq.2) then
2977 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2978 ! s *sqrt(wa_moy(ig,l))
2979 ! else if (idetr.eq.4) then
2980 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2981 ! s *wa_moy(ig,l)
2982 ! endif
2983 END IF
2984 END DO
2985 END DO
2986
2987 ! print*,'10 OK convect8'
2988 ! print*,'WA2 ',wa_moy
2989 ! calcul de la fraction de la maille concern�e par l'ascendance en tenant
2990 ! compte de l'epluchage du thermique.
2991
2992 ! CR def de zmix continu (profil parabolique des vitesses)
2993 ✗ DO ig = 1, ngrid
2994 ✗ IF (lmix(ig)>1.) THEN
2995 zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig))) &
2996 **2-(zlev(ig,lmix(ig)+1))**2)-(zw2(ig,lmix(ig))-zw2(ig, &
2997 lmix(ig)+1))*((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))/ &
2998 (2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
2999 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))*((zlev( &
3000 ✗ ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
3001 ELSE
3002 ✗ zmix(ig) = 0.
3003 END IF
3004 END DO
3005
3006 ! calcul du nouveau lmix correspondant
3007 ✗ DO ig = 1, ngrid
3008 ✗ DO l = 1, klev
3009 ✗ IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
3010 ✗ lmix(ig) = l
3011 END IF
3012 END DO
3013 END DO
3014
3015 ✗ DO l = 2, nlay
3016 ✗ DO ig = 1, ngrid
3017 ✗ IF (larg_cons(ig,l)>1.) THEN
3018 ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
3019 ✗ fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
3020 ! test
3021 ✗ fraca(ig, l) = max(fraca(ig,l), 0.)
3022 ✗ fraca(ig, l) = min(fraca(ig,l), 0.5)
3023 ✗ fracd(ig, l) = 1. - fraca(ig, l)
3024 ✗ fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
3025 ELSE
3026 ! wa_moy(ig,l)=0.
3027 ✗ fraca(ig, l) = 0.
3028 ✗ fracc(ig, l) = 0.
3029 ✗ fracd(ig, l) = 1.
3030 END IF
3031 END DO
3032 END DO
3033 ! CR: calcul de fracazmix
3034 ✗ DO ig = 1, ngrid
3035 fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
3036 (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
3037 fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca(ig &
3038 ✗ ,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
3039 END DO
3040
3041 ✗ DO l = 2, nlay
3042 ✗ DO ig = 1, ngrid
3043 ✗ IF (larg_cons(ig,l)>1.) THEN
3044 ✗ IF (l>lmix(ig)) THEN
3045 ✗ xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
3046 ✗ IF (idetr==0) THEN
3047 ✗ fraca(ig, l) = fracazmix(ig)
3048 ✗ ELSE IF (idetr==1) THEN
3049 ✗ fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
3050 ✗ ELSE IF (idetr==2) THEN
3051 ✗ fraca(ig, l) = fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
3052 ELSE
3053 ✗ fraca(ig, l) = fracazmix(ig)*xxx(ig, l)**2
3054 END IF
3055 ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
3056 ✗ fraca(ig, l) = max(fraca(ig,l), 0.)
3057 ✗ fraca(ig, l) = min(fraca(ig,l), 0.5)
3058 ✗ fracd(ig, l) = 1. - fraca(ig, l)
3059 ✗ fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
3060 END IF
3061 END IF
3062 END DO
3063 END DO
3064
3065 ! print*,'11 OK convect8'
3066 ! print*,'Ea3 ',wa_moy
3067 ! ------------------------------------------------------------------
3068 ! Calcul de fracd, wd
3069 ! somme wa - wd = 0
3070 ! ------------------------------------------------------------------
3071
3072
3073 ✗ DO ig = 1, ngrid
3074 ✗ fm(ig, 1) = 0.
3075 ✗ fm(ig, nlay+1) = 0.
3076 END DO
3077
3078 ✗ DO l = 2, nlay
3079 ✗ DO ig = 1, ngrid
3080 ✗ fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
3081 ! CR:test
3082 ✗ IF (entr(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) THEN
3083 ✗ fm(ig, l) = fm(ig, l-1)
3084 ! write(1,*)'ajustement fm, l',l
3085 END IF
3086 ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
3087 ! RC
3088 END DO
3089 ✗ DO ig = 1, ngrid
3090 ✗ IF (fracd(ig,l)<0.1) THEN
3091 ✗ abort_message = 'fracd trop petit'
3092 ✗ CALL abort_physic(modname, abort_message, 1)
3093 ELSE
3094 ! vitesse descendante "diagnostique"
3095 ✗ wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
3096 END IF
3097 END DO
3098 END DO
3099
3100 ✗ DO l = 1, nlay
3101 ✗ DO ig = 1, ngrid
3102 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
3103 ✗ masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
3104 END DO
3105 END DO
3106
3107 ! print*,'12 OK convect8'
3108 ! print*,'WA4 ',wa_moy
3109 ! c------------------------------------------------------------------
3110 ! calcul du transport vertical
3111 ! ------------------------------------------------------------------
3112
3113 GO TO 4444
3114 ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
3115 DO l = 2, nlay - 1
3116 DO ig = 1, ngrid
3117 IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( &
3118 ig,l+1)) THEN
3119 ! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
3120 ! s ,fm(ig,l+1)*ptimestep
3121 ! s ,' M=',masse(ig,l),masse(ig,l+1)
3122 END IF
3123 END DO
3124 END DO
3125
3126 DO l = 1, nlay
3127 DO ig = 1, ngrid
3128 IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
3129 ! print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
3130 ! s ,entr(ig,l)*ptimestep
3131 ! s ,' M=',masse(ig,l)
3132 END IF
3133 END DO
3134 END DO
3135
3136 DO l = 1, nlay
3137 DO ig = 1, ngrid
3138 IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
3139 ! print*,'WARN!!! fm exagere ig=',ig,' l=',l
3140 ! s ,' FM=',fm(ig,l)
3141 END IF
3142 IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
3143 ! print*,'WARN!!! masse exagere ig=',ig,' l=',l
3144 ! s ,' M=',masse(ig,l)
3145 ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
3146 ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
3147 ! print*,'zlev(ig,l+1),zlev(ig,l)'
3148 ! s ,zlev(ig,l+1),zlev(ig,l)
3149 ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
3150 ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
3151 END IF
3152 IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
3153 ! print*,'WARN!!! entr exagere ig=',ig,' l=',l
3154 ! s ,' E=',entr(ig,l)
3155 END IF
3156 END DO
3157 END DO
3158
3159 4444 CONTINUE
3160
3161 ✗ IF (w2di==1) THEN
3162 ✗ fm0 = fm0 + ptimestep*(fm-fm0)/tho
3163 ✗ entr0 = entr0 + ptimestep*(entr-entr0)/tho
3164 ELSE
3165 ✗ fm0 = fm
3166 ✗ entr0 = entr
3167 END IF
3168
3169 IF (1==1) THEN
3170 ! call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3171 ! . ,zh,zdhadj,zha)
3172 ! call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3173 ! . ,zo,pdoadj,zoa)
3174 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zthl, &
3175 ✗ zdthladj, zta)
3176 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, po, pdoadj, &
3177 ✗ zoa)
3178 ELSE
3179 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
3180 zdhadj, zha)
3181 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
3182 pdoadj, zoa)
3183 END IF
3184
3185 IF (1==0) THEN
3186 CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
3187 zu, zv, pduadj, pdvadj, zua, zva)
3188 ELSE
3189 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
3190 ✗ zua)
3191 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
3192 ✗ zva)
3193 END IF
3194
3195 ✗ DO l = 1, nlay
3196 ✗ DO ig = 1, ngrid
3197 ✗ zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
3198 ✗ zf2 = zf/(1.-zf)
3199 ✗ thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l))**2
3200 ✗ wth2(ig, l) = zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
3201 END DO
3202 END DO
3203
3204
3205
3206 ! print*,'13 OK convect8'
3207 ! print*,'WA5 ',wa_moy
3208 ✗ DO l = 1, nlay
3209 ✗ DO ig = 1, ngrid
3210 ! pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
3211 ✗ pdtadj(ig, l) = zdthladj(ig, l)*zpspsk(ig, l)
3212 END DO
3213 END DO
3214
3215
3216 ! do l=1,nlay
3217 ! do ig=1,ngrid
3218 ! if(abs(pdtadj(ig,l))*86400..gt.500.) then
3219 ! print*,'WARN!!! ig=',ig,' l=',l
3220 ! s ,' pdtadj=',pdtadj(ig,l)
3221 ! endif
3222 ! if(abs(pdoadj(ig,l))*86400..gt.1.) then
3223 ! print*,'WARN!!! ig=',ig,' l=',l
3224 ! s ,' pdoadj=',pdoadj(ig,l)
3225 ! endif
3226 ! enddo
3227 ! enddo
3228
3229 ! print*,'14 OK convect8'
3230 ! ------------------------------------------------------------------
3231 ! Calculs pour les sorties
3232 ! ------------------------------------------------------------------
3233
3234 ✗ RETURN
3235 END SUBROUTINE thermcell_eau
3236
3237 ✗ SUBROUTINE thermcell(ngrid, nlay, ptimestep, pplay, pplev, pphi, pu, pv, pt, &
3238 ✗ po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0 & ! s
3239 ! ,pu_therm,pv_therm
3240 , r_aspect, l_mix, w2di, tho)
3241
3242 USE dimphy
3243 IMPLICIT NONE
3244
3245 ! =======================================================================
3246
3247 ! Calcul du transport verticale dans la couche limite en presence
3248 ! de "thermiques" explicitement representes
3249
3250 ! R��criture � partir d'un listing papier � Habas, le 14/02/00
3251
3252 ! le thermique est suppos� homog�ne et dissip� par m�lange avec
3253 ! son environnement. la longueur l_mix contr�le l'efficacit� du
3254 ! m�lange
3255
3256 ! Le calcul du transport des diff�rentes esp�ces se fait en prenant
3257 ! en compte:
3258 ! 1. un flux de masse montant
3259 ! 2. un flux de masse descendant
3260 ! 3. un entrainement
3261 ! 4. un detrainement
3262
3263 ! =======================================================================
3264
3265 ! -----------------------------------------------------------------------
3266 ! declarations:
3267 ! -------------
3268
3269 include "YOMCST.h"
3270
3271 ! arguments:
3272 ! ----------
3273
3274 INTEGER ngrid, nlay, w2di
3275 REAL tho
3276 REAL ptimestep, l_mix, r_aspect
3277 REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
3278 REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
3279 REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
3280 REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
3281 REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
3282 REAL pphi(ngrid, nlay)
3283
3284 INTEGER idetr
3285 SAVE idetr
3286 DATA idetr/3/
3287 !$OMP THREADPRIVATE(idetr)
3288
3289 ! local:
3290 ! ------
3291
3292 ✗ INTEGER ig, k, l, lmaxa(klon), lmix(klon)
3293 REAL zsortie1d(klon)
3294 ! CR: on remplace lmax(klon,klev+1)
3295 ✗ INTEGER lmax(klon), lmin(klon), lentr(klon)
3296 ✗ REAL linter(klon)
3297 ✗ REAL zmix(klon), fracazmix(klon)
3298 ! RC
3299 ✗ REAL zmax(klon), zw, zz, zw2(klon, klev+1), ztva(klon, klev), zzz
3300
3301 ✗ REAL zlev(klon, klev+1), zlay(klon, klev)
3302 ✗ REAL zh(klon, klev), zdhadj(klon, klev)
3303 ✗ REAL ztv(klon, klev)
3304 ✗ REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
3305 REAL wh(klon, klev+1)
3306 REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
3307 REAL zla(klon, klev+1)
3308 REAL zwa(klon, klev+1)
3309 REAL zld(klon, klev+1)
3310 REAL zwd(klon, klev+1)
3311 REAL zsortie(klon, klev)
3312 ✗ REAL zva(klon, klev)
3313 ✗ REAL zua(klon, klev)
3314 ✗ REAL zoa(klon, klev)
3315
3316 ✗ REAL zha(klon, klev)
3317 ✗ REAL wa_moy(klon, klev+1)
3318 ✗ REAL fraca(klon, klev+1)
3319 ✗ REAL fracc(klon, klev+1)
3320 REAL zf, zf2
3321 ✗ REAL thetath2(klon, klev), wth2(klon, klev)
3322 ! common/comtherm/thetath2,wth2
3323
3324 REAL count_time
3325 INTEGER ialt
3326
3327 LOGICAL sorties
3328 ✗ REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
3329 ✗ REAL zpspsk(klon, klev)
3330
3331 ! real wmax(klon,klev),wmaxa(klon)
3332 ✗ REAL wmax(klon), wmaxa(klon)
3333 ✗ REAL wa(klon, klev, klev+1)
3334 ✗ REAL wd(klon, klev+1)
3335 REAL larg_part(klon, klev, klev+1)
3336 ✗ REAL fracd(klon, klev+1)
3337 ✗ REAL xxx(klon, klev+1)
3338 ✗ REAL larg_cons(klon, klev+1)
3339 ✗ REAL larg_detr(klon, klev+1)
3340 ✗ REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
3341 REAL pu_therm(klon, klev), pv_therm(klon, klev)
3342 ✗ REAL fm(klon, klev+1), entr(klon, klev)
3343 ✗ REAL fmc(klon, klev+1)
3344
3345 ! CR:nouvelles variables
3346 ✗ REAL f_star(klon, klev+1), entr_star(klon, klev)
3347 ✗ REAL entr_star_tot(klon), entr_star2(klon)
3348 ✗ REAL f(klon), f0(klon)
3349 ✗ REAL zlevinter(klon)
3350 LOGICAL first
3351 DATA first/.FALSE./
3352 SAVE first
3353 !$OMP THREADPRIVATE(first)
3354 ! RC
3355
3356 CHARACTER *2 str2
3357 CHARACTER *10 str10
3358
3359 CHARACTER (LEN=20) :: modname = 'thermcell'
3360 CHARACTER (LEN=80) :: abort_message
3361
3362 LOGICAL vtest(klon), down
3363
3364 EXTERNAL scopy
3365
3366 INTEGER ncorrec, ll
3367 SAVE ncorrec
3368 DATA ncorrec/0/
3369 !$OMP THREADPRIVATE(ncorrec)
3370
3371
3372 ! -----------------------------------------------------------------------
3373 ! initialisation:
3374 ! ---------------
3375
3376 sorties = .TRUE.
3377 ✗ IF (ngrid/=klon) THEN
3378 ✗ PRINT *
3379 ✗ PRINT *, 'STOP dans convadj'
3380 ✗ PRINT *, 'ngrid =', ngrid
3381 ✗ PRINT *, 'klon =', klon
3382 END IF
3383
3384 ! -----------------------------------------------------------------------
3385 ! incrementation eventuelle de tendances precedentes:
3386 ! ---------------------------------------------------
3387
3388 ! print*,'0 OK convect8'
3389
3390 ✗ DO l = 1, nlay
3391 ✗ DO ig = 1, ngrid
3392 ✗ zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa
3393 ✗ zh(ig, l) = pt(ig, l)/zpspsk(ig, l)
3394 ✗ zu(ig, l) = pu(ig, l)
3395 ✗ zv(ig, l) = pv(ig, l)
3396 ✗ zo(ig, l) = po(ig, l)
3397 ✗ ztv(ig, l) = zh(ig, l)*(1.+0.61*zo(ig,l))
3398 END DO
3399 END DO
3400
3401 ! print*,'1 OK convect8'
3402 ! --------------------
3403
3404
3405 ! + + + + + + + + + + +
3406
3407
3408 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
3409 ! wh,wt,wo ...
3410
3411 ! + + + + + + + + + + + zh,zu,zv,zo,rho
3412
3413
3414 ! -------------------- zlev(1)
3415 ! \\\\\\\\\\\\\\\\\\\
3416
3417
3418 ! -----------------------------------------------------------------------
3419 ! Calcul des altitudes des couches
3420 ! -----------------------------------------------------------------------
3421
3422 ✗ DO l = 2, nlay
3423 ✗ DO ig = 1, ngrid
3424 ✗ zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
3425 END DO
3426 END DO
3427 ✗ DO ig = 1, ngrid
3428 ✗ zlev(ig, 1) = 0.
3429 ✗ zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
3430 END DO
3431 ✗ DO l = 1, nlay
3432 ✗ DO ig = 1, ngrid
3433 ✗ zlay(ig, l) = pphi(ig, l)/rg
3434 END DO
3435 END DO
3436
3437 ! print*,'2 OK convect8'
3438 ! -----------------------------------------------------------------------
3439 ! Calcul des densites
3440 ! -----------------------------------------------------------------------
3441
3442 ✗ DO l = 1, nlay
3443 ✗ DO ig = 1, ngrid
3444 ✗ rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*zh(ig,l))
3445 END DO
3446 END DO
3447
3448 ✗ DO l = 2, nlay
3449 ✗ DO ig = 1, ngrid
3450 ✗ rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
3451 END DO
3452 END DO
3453
3454 ✗ DO k = 1, nlay
3455 ✗ DO l = 1, nlay + 1
3456 ✗ DO ig = 1, ngrid
3457 ✗ wa(ig, k, l) = 0.
3458 END DO
3459 END DO
3460 END DO
3461
3462 ! print*,'3 OK convect8'
3463 ! ------------------------------------------------------------------
3464 ! Calcul de w2, quarre de w a partir de la cape
3465 ! a partir de w2, on calcule wa, vitesse de l'ascendance
3466
3467 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
3468 ! w2 est stoke dans wa
3469
3470 ! ATTENTION: dans convect8, on n'utilise le calcule des wa
3471 ! independants par couches que pour calculer l'entrainement
3472 ! a la base et la hauteur max de l'ascendance.
3473
3474 ! Indicages:
3475 ! l'ascendance provenant du niveau k traverse l'interface l avec
3476 ! une vitesse wa(k,l).
3477
3478 ! --------------------
3479
3480 ! + + + + + + + + + +
3481
3482 ! wa(k,l) ---- -------------------- l
3483 ! / ! /||\ + + + + + + + + + +
3484 ! ||
3485 ! || --------------------
3486 ! ||
3487 ! || + + + + + + + + + +
3488 ! ||
3489 ! || --------------------
3490 ! ||__
3491 ! |___ + + + + + + + + + + k
3492
3493 ! --------------------
3494
3495
3496
3497 ! ------------------------------------------------------------------
3498
3499 ! CR: ponderation entrainement des couches instables
3500 ! def des entr_star tels que entr=f*entr_star
3501 ✗ DO l = 1, klev
3502 ✗ DO ig = 1, ngrid
3503 ✗ entr_star(ig, l) = 0.
3504 END DO
3505 END DO
3506 ! determination de la longueur de la couche d entrainement
3507 ✗ DO ig = 1, ngrid
3508 ✗ lentr(ig) = 1
3509 END DO
3510
3511 ! on ne considere que les premieres couches instables
3512 ✗ DO k = nlay - 2, 1, -1
3513 ✗ DO ig = 1, ngrid
3514 ✗ IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN
3515 ✗ lentr(ig) = k
3516 END IF
3517 END DO
3518 END DO
3519
3520 ! determination du lmin: couche d ou provient le thermique
3521 ✗ DO ig = 1, ngrid
3522 ✗ lmin(ig) = 1
3523 END DO
3524 ✗ DO ig = 1, ngrid
3525 ✗ DO l = nlay, 2, -1
3526 ✗ IF (ztv(ig,l-1)>ztv(ig,l)) THEN
3527 ✗ lmin(ig) = l - 1
3528 END IF
3529 END DO
3530 END DO
3531
3532 ! definition de l'entrainement des couches
3533 ✗ DO l = 1, klev - 1
3534 ✗ DO ig = 1, ngrid
3535 ✗ IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<=lentr(ig)) THEN
3536 ✗ entr_star(ig, l) = (ztv(ig,l)-ztv(ig,l+1))*(zlev(ig,l+1)-zlev(ig,l))
3537 END IF
3538 END DO
3539 END DO
3540 ! pas de thermique si couches 1->5 stables
3541 ✗ DO ig = 1, ngrid
3542 ✗ IF (lmin(ig)>5) THEN
3543 ✗ DO l = 1, klev
3544 ✗ entr_star(ig, l) = 0.
3545 END DO
3546 END IF
3547 END DO
3548 ! calcul de l entrainement total
3549 ✗ DO ig = 1, ngrid
3550 ✗ entr_star_tot(ig) = 0.
3551 END DO
3552 ✗ DO ig = 1, ngrid
3553 ✗ DO k = 1, klev
3554 ✗ entr_star_tot(ig) = entr_star_tot(ig) + entr_star(ig, k)
3555 END DO
3556 END DO
3557
3558 ✗ PRINT *, 'fin calcul entr_star'
3559 ✗ DO k = 1, klev
3560 ✗ DO ig = 1, ngrid
3561 ✗ ztva(ig, k) = ztv(ig, k)
3562 END DO
3563 END DO
3564 ! RC
3565 ! print*,'7 OK convect8'
3566 ✗ DO k = 1, klev + 1
3567 ✗ DO ig = 1, ngrid
3568 ✗ zw2(ig, k) = 0.
3569 ✗ fmc(ig, k) = 0.
3570 ! CR
3571 ✗ f_star(ig, k) = 0.
3572 ! RC
3573 ✗ larg_cons(ig, k) = 0.
3574 ✗ larg_detr(ig, k) = 0.
3575 ✗ wa_moy(ig, k) = 0.
3576 END DO
3577 END DO
3578
3579 ! print*,'8 OK convect8'
3580 ✗ DO ig = 1, ngrid
3581 ✗ linter(ig) = 1.
3582 ✗ lmaxa(ig) = 1
3583 ✗ lmix(ig) = 1
3584 ✗ wmaxa(ig) = 0.
3585 END DO
3586
3587 ! CR:
3588 ✗ DO l = 1, nlay - 2
3589 ✗ DO ig = 1, ngrid
3590 ✗ IF (ztv(ig,l)>ztv(ig,l+1) .AND. entr_star(ig,l)>1.E-10 .AND. &
3591 zw2(ig,l)<1E-10) THEN
3592 ✗ f_star(ig, l+1) = entr_star(ig, l)
3593 ! test:calcul de dteta
3594 zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
3595 ✗ (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
3596 ✗ larg_detr(ig, l) = 0.
3597 ✗ ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+entr_star(ig, &
3598 l)>1.E-10)) THEN
3599 ✗ f_star(ig, l+1) = f_star(ig, l) + entr_star(ig, l)
3600 ztva(ig, l) = (f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)*ztv(ig,l))/ &
3601 ✗ f_star(ig, l+1)
3602 zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/f_star(ig,l+1))**2 + &
3603 ✗ 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
3604 END IF
3605 ! determination de zmax continu par interpolation lineaire
3606 ✗ IF (zw2(ig,l+1)<0.) THEN
3607 ! test
3608 ✗ IF (abs(zw2(ig,l+1)-zw2(ig,l))<1E-10) THEN
3609 ✗ PRINT *, 'pb linter'
3610 END IF
3611 linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
3612 ✗ ig,l))
3613 ✗ zw2(ig, l+1) = 0.
3614 ✗ lmaxa(ig) = l
3615 ELSE
3616 IF (zw2(ig,l+1)<0.) THEN
3617 PRINT *, 'pb1 zw2<0'
3618 END IF
3619 ✗ wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
3620 END IF
3621 ✗ IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
3622 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
3623 ✗ lmix(ig) = l + 1
3624 ✗ wmaxa(ig) = wa_moy(ig, l+1)
3625 END IF
3626 END DO
3627 END DO
3628 ✗ PRINT *, 'fin calcul zw2'
3629
3630 ! Calcul de la couche correspondant a la hauteur du thermique
3631 ✗ DO ig = 1, ngrid
3632 ✗ lmax(ig) = lentr(ig)
3633 END DO
3634 ✗ DO ig = 1, ngrid
3635 ✗ DO l = nlay, lentr(ig) + 1, -1
3636 ✗ IF (zw2(ig,l)<=1.E-10) THEN
3637 ✗ lmax(ig) = l - 1
3638 END IF
3639 END DO
3640 END DO
3641 ! pas de thermique si couches 1->5 stables
3642 ✗ DO ig = 1, ngrid
3643 ✗ IF (lmin(ig)>5) THEN
3644 ✗ lmax(ig) = 1
3645 ✗ lmin(ig) = 1
3646 END IF
3647 END DO
3648
3649 ! Determination de zw2 max
3650 ✗ DO ig = 1, ngrid
3651 ✗ wmax(ig) = 0.
3652 END DO
3653
3654 ✗ DO l = 1, nlay
3655 ✗ DO ig = 1, ngrid
3656 ✗ IF (l<=lmax(ig)) THEN
3657 ✗ IF (zw2(ig,l)<0.) THEN
3658 ✗ PRINT *, 'pb2 zw2<0'
3659 END IF
3660 ✗ zw2(ig, l) = sqrt(zw2(ig,l))
3661 ✗ wmax(ig) = max(wmax(ig), zw2(ig,l))
3662 ELSE
3663 ✗ zw2(ig, l) = 0.
3664 END IF
3665 END DO
3666 END DO
3667
3668 ! Longueur caracteristique correspondant a la hauteur des thermiques.
3669 ✗ DO ig = 1, ngrid
3670 ✗ zmax(ig) = 0.
3671 ✗ zlevinter(ig) = zlev(ig, 1)
3672 END DO
3673 ✗ DO ig = 1, ngrid
3674 ! calcul de zlevinter
3675 zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
3676 ✗ zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
3677 ✗ zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,lmin(ig)))
3678 END DO
3679
3680 ✗ PRINT *, 'avant fermeture'
3681 ! Fermeture,determination de f
3682 ✗ DO ig = 1, ngrid
3683 ✗ entr_star2(ig) = 0.
3684 END DO
3685 ✗ DO ig = 1, ngrid
3686 ✗ IF (entr_star_tot(ig)<1.E-10) THEN
3687 ✗ f(ig) = 0.
3688 ELSE
3689 ✗ DO k = lmin(ig), lentr(ig)
3690 entr_star2(ig) = entr_star2(ig) + entr_star(ig, k)**2/(rho(ig,k)*( &
3691 ✗ zlev(ig,k+1)-zlev(ig,k)))
3692 END DO
3693 ! Nouvelle fermeture
3694 f(ig) = wmax(ig)/(max(500.,zmax(ig))*r_aspect*entr_star2(ig))* &
3695 ✗ entr_star_tot(ig)
3696 ! test
3697 ! if (first) then
3698 ! f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
3699 ! s *wmax(ig))
3700 ! endif
3701 END IF
3702 ! f0(ig)=f(ig)
3703 ! first=.true.
3704 END DO
3705 ✗ PRINT *, 'apres fermeture'
3706
3707 ! Calcul de l'entrainement
3708 ✗ DO k = 1, klev
3709 ✗ DO ig = 1, ngrid
3710 ✗ entr(ig, k) = f(ig)*entr_star(ig, k)
3711 END DO
3712 END DO
3713 ! Calcul des flux
3714 ✗ DO ig = 1, ngrid
3715 ✗ DO l = 1, lmax(ig) - 1
3716 ✗ fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
3717 END DO
3718 END DO
3719
3720 ! RC
3721
3722
3723 ! print*,'9 OK convect8'
3724 ! print*,'WA1 ',wa_moy
3725
3726 ! determination de l'indice du debut de la mixed layer ou w decroit
3727
3728 ! calcul de la largeur de chaque ascendance dans le cas conservatif.
3729 ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
3730 ! d'une couche est �gale � la hauteur de la couche alimentante.
3731 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
3732 ! de la vitesse d'entrainement horizontal dans la couche alimentante.
3733
3734 ✗ DO l = 2, nlay
3735 ✗ DO ig = 1, ngrid
3736 ✗ IF (l<=lmaxa(ig)) THEN
3737 ✗ zw = max(wa_moy(ig,l), 1.E-10)
3738 ✗ larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw)
3739 END IF
3740 END DO
3741 END DO
3742
3743 ✗ DO l = 2, nlay
3744 ✗ DO ig = 1, ngrid
3745 ✗ IF (l<=lmaxa(ig)) THEN
3746 ! if (idetr.eq.0) then
3747 ! cette option est finalement en dur.
3748 ✗ IF ((l_mix*zlev(ig,l))<0.) THEN
3749 ✗ PRINT *, 'pb l_mix*zlev<0'
3750 END IF
3751 ✗ larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
3752 ! else if (idetr.eq.1) then
3753 ! larg_detr(ig,l)=larg_cons(ig,l)
3754 ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
3755 ! else if (idetr.eq.2) then
3756 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3757 ! s *sqrt(wa_moy(ig,l))
3758 ! else if (idetr.eq.4) then
3759 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3760 ! s *wa_moy(ig,l)
3761 ! endif
3762 END IF
3763 END DO
3764 END DO
3765
3766 ! print*,'10 OK convect8'
3767 ! print*,'WA2 ',wa_moy
3768 ! calcul de la fraction de la maille concern�e par l'ascendance en tenant
3769 ! compte de l'epluchage du thermique.
3770
3771 ! CR def de zmix continu (profil parabolique des vitesses)
3772 ✗ DO ig = 1, ngrid
3773 ✗ IF (lmix(ig)>1.) THEN
3774 ! test
3775 ✗ IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
3776 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
3777 zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- &
3778 (zlev(ig,lmix(ig)))))>1E-10) THEN
3779
3780 zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) &
3781 )**2-(zlev(ig,lmix(ig)+1))**2)-(zw2(ig,lmix(ig))-zw2(ig, &
3782 lmix(ig)+1))*((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))/ &
3783 (2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
3784 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
3785 ✗ zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
3786 ELSE
3787 ✗ zmix(ig) = zlev(ig, lmix(ig))
3788 ✗ PRINT *, 'pb zmix'
3789 END IF
3790 ELSE
3791 ✗ zmix(ig) = 0.
3792 END IF
3793 ! test
3794 ✗ IF ((zmax(ig)-zmix(ig))<0.) THEN
3795 ✗ zmix(ig) = 0.99*zmax(ig)
3796 ! print*,'pb zmix>zmax'
3797 END IF
3798 END DO
3799
3800 ! calcul du nouveau lmix correspondant
3801 ✗ DO ig = 1, ngrid
3802 ✗ DO l = 1, klev
3803 ✗ IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
3804 ✗ lmix(ig) = l
3805 END IF
3806 END DO
3807 END DO
3808
3809 ✗ DO l = 2, nlay
3810 ✗ DO ig = 1, ngrid
3811 ✗ IF (larg_cons(ig,l)>1.) THEN
3812 ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
3813 ✗ fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
3814 ! test
3815 ✗ fraca(ig, l) = max(fraca(ig,l), 0.)
3816 ✗ fraca(ig, l) = min(fraca(ig,l), 0.5)
3817 ✗ fracd(ig, l) = 1. - fraca(ig, l)
3818 ✗ fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
3819 ELSE
3820 ! wa_moy(ig,l)=0.
3821 ✗ fraca(ig, l) = 0.
3822 ✗ fracc(ig, l) = 0.
3823 ✗ fracd(ig, l) = 1.
3824 END IF
3825 END DO
3826 END DO
3827 ! CR: calcul de fracazmix
3828 ✗ DO ig = 1, ngrid
3829 fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
3830 (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
3831 fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca(ig &
3832 ✗ ,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
3833 END DO
3834
3835 ✗ DO l = 2, nlay
3836 ✗ DO ig = 1, ngrid
3837 ✗ IF (larg_cons(ig,l)>1.) THEN
3838 ✗ IF (l>lmix(ig)) THEN
3839 ! test
3840 ✗ IF (zmax(ig)-zmix(ig)<1.E-10) THEN
3841 ! print*,'pb xxx'
3842 ✗ xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
3843 ELSE
3844 ✗ xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
3845 END IF
3846 ✗ IF (idetr==0) THEN
3847 ✗ fraca(ig, l) = fracazmix(ig)
3848 ✗ ELSE IF (idetr==1) THEN
3849 ✗ fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
3850 ✗ ELSE IF (idetr==2) THEN
3851 ✗ fraca(ig, l) = fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
3852 ELSE
3853 ✗ fraca(ig, l) = fracazmix(ig)*xxx(ig, l)**2
3854 END IF
3855 ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
3856 ✗ fraca(ig, l) = max(fraca(ig,l), 0.)
3857 ✗ fraca(ig, l) = min(fraca(ig,l), 0.5)
3858 ✗ fracd(ig, l) = 1. - fraca(ig, l)
3859 ✗ fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
3860 END IF
3861 END IF
3862 END DO
3863 END DO
3864
3865 ✗ PRINT *, 'fin calcul fraca'
3866 ! print*,'11 OK convect8'
3867 ! print*,'Ea3 ',wa_moy
3868 ! ------------------------------------------------------------------
3869 ! Calcul de fracd, wd
3870 ! somme wa - wd = 0
3871 ! ------------------------------------------------------------------
3872
3873
3874 ✗ DO ig = 1, ngrid
3875 ✗ fm(ig, 1) = 0.
3876 ✗ fm(ig, nlay+1) = 0.
3877 END DO
3878
3879 ✗ DO l = 2, nlay
3880 ✗ DO ig = 1, ngrid
3881 ✗ fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
3882 ! CR:test
3883 ✗ IF (entr(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) THEN
3884 ✗ fm(ig, l) = fm(ig, l-1)
3885 ! write(1,*)'ajustement fm, l',l
3886 END IF
3887 ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
3888 ! RC
3889 END DO
3890 ✗ DO ig = 1, ngrid
3891 ✗ IF (fracd(ig,l)<0.1) THEN
3892 ✗ abort_message = 'fracd trop petit'
3893 ✗ CALL abort_physic(modname, abort_message, 1)
3894 ELSE
3895 ! vitesse descendante "diagnostique"
3896 ✗ wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
3897 END IF
3898 END DO
3899 END DO
3900
3901 ✗ DO l = 1, nlay
3902 ✗ DO ig = 1, ngrid
3903 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
3904 ✗ masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
3905 END DO
3906 END DO
3907
3908 ! print*,'12 OK convect8'
3909 ! print*,'WA4 ',wa_moy
3910 ! c------------------------------------------------------------------
3911 ! calcul du transport vertical
3912 ! ------------------------------------------------------------------
3913
3914 GO TO 4444
3915 ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
3916 DO l = 2, nlay - 1
3917 DO ig = 1, ngrid
3918 IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( &
3919 ig,l+1)) THEN
3920 ! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
3921 ! s ,fm(ig,l+1)*ptimestep
3922 ! s ,' M=',masse(ig,l),masse(ig,l+1)
3923 END IF
3924 END DO
3925 END DO
3926
3927 DO l = 1, nlay
3928 DO ig = 1, ngrid
3929 IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
3930 ! print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
3931 ! s ,entr(ig,l)*ptimestep
3932 ! s ,' M=',masse(ig,l)
3933 END IF
3934 END DO
3935 END DO
3936
3937 DO l = 1, nlay
3938 DO ig = 1, ngrid
3939 IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
3940 ! print*,'WARN!!! fm exagere ig=',ig,' l=',l
3941 ! s ,' FM=',fm(ig,l)
3942 END IF
3943 IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
3944 ! print*,'WARN!!! masse exagere ig=',ig,' l=',l
3945 ! s ,' M=',masse(ig,l)
3946 ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
3947 ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
3948 ! print*,'zlev(ig,l+1),zlev(ig,l)'
3949 ! s ,zlev(ig,l+1),zlev(ig,l)
3950 ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
3951 ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
3952 END IF
3953 IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
3954 ! print*,'WARN!!! entr exagere ig=',ig,' l=',l
3955 ! s ,' E=',entr(ig,l)
3956 END IF
3957 END DO
3958 END DO
3959
3960 4444 CONTINUE
3961
3962 ! CR:redefinition du entr
3963 ✗ DO l = 1, nlay
3964 ✗ DO ig = 1, ngrid
3965 ✗ detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1)
3966 ✗ IF (detr(ig,l)<0.) THEN
3967 ✗ entr(ig, l) = entr(ig, l) - detr(ig, l)
3968 ✗ detr(ig, l) = 0.
3969 ! print*,'WARNING !!! detrainement negatif ',ig,l
3970 END IF
3971 END DO
3972 END DO
3973 ! RC
3974 ✗ IF (w2di==1) THEN
3975 ✗ fm0 = fm0 + ptimestep*(fm-fm0)/tho
3976 ✗ entr0 = entr0 + ptimestep*(entr-entr0)/tho
3977 ELSE
3978 ✗ fm0 = fm
3979 ✗ entr0 = entr
3980 END IF
3981
3982 IF (1==1) THEN
3983 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, &
3984 ✗ zha)
3985 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, &
3986 ✗ zoa)
3987 ELSE
3988 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
3989 zdhadj, zha)
3990 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
3991 pdoadj, zoa)
3992 END IF
3993
3994 IF (1==0) THEN
3995 CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
3996 zu, zv, pduadj, pdvadj, zua, zva)
3997 ELSE
3998 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
3999 ✗ zua)
4000 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
4001 ✗ zva)
4002 END IF
4003
4004 ✗ DO l = 1, nlay
4005 ✗ DO ig = 1, ngrid
4006 ✗ zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
4007 ✗ zf2 = zf/(1.-zf)
4008 ✗ thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l))**2
4009 ✗ wth2(ig, l) = zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
4010 END DO
4011 END DO
4012
4013
4014
4015 ! print*,'13 OK convect8'
4016 ! print*,'WA5 ',wa_moy
4017 ✗ DO l = 1, nlay
4018 ✗ DO ig = 1, ngrid
4019 ✗ pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l)
4020 END DO
4021 END DO
4022
4023
4024 ! do l=1,nlay
4025 ! do ig=1,ngrid
4026 ! if(abs(pdtadj(ig,l))*86400..gt.500.) then
4027 ! print*,'WARN!!! ig=',ig,' l=',l
4028 ! s ,' pdtadj=',pdtadj(ig,l)
4029 ! endif
4030 ! if(abs(pdoadj(ig,l))*86400..gt.1.) then
4031 ! print*,'WARN!!! ig=',ig,' l=',l
4032 ! s ,' pdoadj=',pdoadj(ig,l)
4033 ! endif
4034 ! enddo
4035 ! enddo
4036
4037 ! print*,'14 OK convect8'
4038 ! ------------------------------------------------------------------
4039 ! Calculs pour les sorties
4040 ! ------------------------------------------------------------------
4041
4042 IF (sorties) THEN
4043 DO l = 1, nlay
4044 DO ig = 1, ngrid
4045 zla(ig, l) = (1.-fracd(ig,l))*zmax(ig)
4046 zld(ig, l) = fracd(ig, l)*zmax(ig)
4047 IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ &
4048 (1.-fracd(ig,l))
4049 END DO
4050 END DO
4051
4052 ! deja fait
4053 ! do l=1,nlay
4054 ! do ig=1,ngrid
4055 ! detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
4056 ! if (detr(ig,l).lt.0.) then
4057 ! entr(ig,l)=entr(ig,l)-detr(ig,l)
4058 ! detr(ig,l)=0.
4059 ! print*,'WARNING !!! detrainement negatif ',ig,l
4060 ! endif
4061 ! enddo
4062 ! enddo
4063
4064 ! print*,'15 OK convect8'
4065
4066
4067 ! #define und
4068 GO TO 123
4069 123 CONTINUE
4070
4071 END IF
4072
4073 ! if(wa_moy(1,4).gt.1.e-10) stop
4074
4075 ! print*,'19 OK convect8'
4076 ✗ RETURN
4077 END SUBROUTINE thermcell
4078
4079 ✗ SUBROUTINE dqthermcell(ngrid, nlay, ptimestep, fm, entr, masse, q, dq, qa)
4080 USE dimphy
4081 IMPLICIT NONE
4082
4083 ! =======================================================================
4084
4085 ! Calcul du transport verticale dans la couche limite en presence
4086 ! de "thermiques" explicitement representes
4087 ! calcul du dq/dt une fois qu'on connait les ascendances
4088
4089 ! =======================================================================
4090
4091 INTEGER ngrid, nlay
4092
4093 REAL ptimestep
4094 REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4095 REAL entr(ngrid, nlay)
4096 REAL q(ngrid, nlay)
4097 REAL dq(ngrid, nlay)
4098
4099 ✗ REAL qa(klon, klev), detr(klon, klev), wqd(klon, klev+1)
4100
4101 INTEGER ig, k
4102
4103 ! calcul du detrainement
4104
4105 ✗ DO k = 1, nlay
4106 ✗ DO ig = 1, ngrid
4107 ✗ detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4108 ! test
4109 ✗ IF (detr(ig,k)<0.) THEN
4110 ✗ entr(ig, k) = entr(ig, k) - detr(ig, k)
4111 ✗ detr(ig, k) = 0.
4112 ! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
4113 ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
4114 END IF
4115 IF (fm(ig,k+1)<0.) THEN
4116 ! print*,'fm2<0!!!'
4117 END IF
4118 ✗ IF (entr(ig,k)<0.) THEN
4119 ! print*,'entr2<0!!!'
4120 END IF
4121 END DO
4122 END DO
4123
4124 ! calcul de la valeur dans les ascendances
4125 ✗ DO ig = 1, ngrid
4126 ✗ qa(ig, 1) = q(ig, 1)
4127 END DO
4128
4129 ✗ DO k = 2, nlay
4130 ✗ DO ig = 1, ngrid
4131 ✗ IF ((fm(ig,k+1)+detr(ig,k))*ptimestep>1.E-5*masse(ig,k)) THEN
4132 qa(ig, k) = (fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))/ &
4133 ✗ (fm(ig,k+1)+detr(ig,k))
4134 ELSE
4135 ✗ qa(ig, k) = q(ig, k)
4136 END IF
4137 IF (qa(ig,k)<0.) THEN
4138 ! print*,'qa<0!!!'
4139 END IF
4140 ✗ IF (q(ig,k)<0.) THEN
4141 ! print*,'q<0!!!'
4142 END IF
4143 END DO
4144 END DO
4145
4146 ✗ DO k = 2, nlay
4147 ✗ DO ig = 1, ngrid
4148 ! wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
4149 ✗ wqd(ig, k) = fm(ig, k)*q(ig, k)
4150 ✗ IF (wqd(ig,k)<0.) THEN
4151 ! print*,'wqd<0!!!'
4152 END IF
4153 END DO
4154 END DO
4155 ✗ DO ig = 1, ngrid
4156 ✗ wqd(ig, 1) = 0.
4157 ✗ wqd(ig, nlay+1) = 0.
4158 END DO
4159
4160 ✗ DO k = 1, nlay
4161 ✗ DO ig = 1, ngrid
4162 dq(ig, k) = (detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)-wqd(ig,k)+wqd(ig,k+ &
4163 ✗ 1))/masse(ig, k)
4164 ! if (dq(ig,k).lt.0.) then
4165 ! print*,'dq<0!!!'
4166 ! endif
4167 END DO
4168 END DO
4169
4170 ✗ RETURN
4171 END SUBROUTINE dqthermcell
4172 ✗ SUBROUTINE dvthermcell(ngrid, nlay, ptimestep, fm, entr, masse, fraca, larga, &
4173 u, v, du, dv, ua, va)
4174 USE dimphy
4175 IMPLICIT NONE
4176
4177 ! =======================================================================
4178
4179 ! Calcul du transport verticale dans la couche limite en presence
4180 ! de "thermiques" explicitement representes
4181 ! calcul du dq/dt une fois qu'on connait les ascendances
4182
4183 ! =======================================================================
4184
4185 INTEGER ngrid, nlay
4186
4187 REAL ptimestep
4188 REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4189 REAL fraca(ngrid, nlay+1)
4190 REAL larga(ngrid)
4191 REAL entr(ngrid, nlay)
4192 REAL u(ngrid, nlay)
4193 REAL ua(ngrid, nlay)
4194 REAL du(ngrid, nlay)
4195 REAL v(ngrid, nlay)
4196 REAL va(ngrid, nlay)
4197 REAL dv(ngrid, nlay)
4198
4199 ✗ REAL qa(klon, klev), detr(klon, klev)
4200 ✗ REAL wvd(klon, klev+1), wud(klon, klev+1)
4201 ✗ REAL gamma0, gamma(klon, klev+1)
4202 REAL dua, dva
4203 INTEGER iter
4204
4205 INTEGER ig, k
4206
4207 ! calcul du detrainement
4208
4209 ✗ DO k = 1, nlay
4210 ✗ DO ig = 1, ngrid
4211 ✗ detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4212 END DO
4213 END DO
4214
4215 ! calcul de la valeur dans les ascendances
4216 ✗ DO ig = 1, ngrid
4217 ✗ ua(ig, 1) = u(ig, 1)
4218 ✗ va(ig, 1) = v(ig, 1)
4219 END DO
4220
4221 ✗ DO k = 2, nlay
4222 ✗ DO ig = 1, ngrid
4223 ✗ IF ((fm(ig,k+1)+detr(ig,k))*ptimestep>1.E-5*masse(ig,k)) THEN
4224 ! On it�re sur la valeur du coeff de freinage.
4225 ! gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
4226 gamma0 = masse(ig, k)*sqrt(0.5*(fraca(ig,k+1)+fraca(ig, &
4227 ✗ k)))*0.5/larga(ig)
4228 ! gamma0=0.
4229 ! la premi�re fois on multiplie le coefficient de freinage
4230 ! par le module du vent dans la couche en dessous.
4231 ✗ dua = ua(ig, k-1) - u(ig, k-1)
4232 ✗ dva = va(ig, k-1) - v(ig, k-1)
4233 ✗ DO iter = 1, 5
4234 ✗ gamma(ig, k) = gamma0*sqrt(dua**2+dva**2)
4235 ua(ig, k) = (fm(ig,k)*ua(ig,k-1)+(entr(ig,k)+gamma(ig, &
4236 ✗ k))*u(ig,k))/(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
4237 va(ig, k) = (fm(ig,k)*va(ig,k-1)+(entr(ig,k)+gamma(ig, &
4238 ✗ k))*v(ig,k))/(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
4239 ! print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
4240 ✗ dua = ua(ig, k) - u(ig, k)
4241 ✗ dva = va(ig, k) - v(ig, k)
4242 END DO
4243 ELSE
4244 ✗ ua(ig, k) = u(ig, k)
4245 ✗ va(ig, k) = v(ig, k)
4246 ✗ gamma(ig, k) = 0.
4247 END IF
4248 END DO
4249 END DO
4250
4251 ✗ DO k = 2, nlay
4252 ✗ DO ig = 1, ngrid
4253 ✗ wud(ig, k) = fm(ig, k)*u(ig, k)
4254 ✗ wvd(ig, k) = fm(ig, k)*v(ig, k)
4255 END DO
4256 END DO
4257 ✗ DO ig = 1, ngrid
4258 ✗ wud(ig, 1) = 0.
4259 ✗ wud(ig, nlay+1) = 0.
4260 ✗ wvd(ig, 1) = 0.
4261 ✗ wvd(ig, nlay+1) = 0.
4262 END DO
4263
4264 ✗ DO k = 1, nlay
4265 ✗ DO ig = 1, ngrid
4266 du(ig, k) = ((detr(ig,k)+gamma(ig,k))*ua(ig,k)-(entr(ig,k)+gamma(ig, &
4267 ✗ k))*u(ig,k)-wud(ig,k)+wud(ig,k+1))/masse(ig, k)
4268 dv(ig, k) = ((detr(ig,k)+gamma(ig,k))*va(ig,k)-(entr(ig,k)+gamma(ig, &
4269 ✗ k))*v(ig,k)-wvd(ig,k)+wvd(ig,k+1))/masse(ig, k)
4270 END DO
4271 END DO
4272
4273 ✗ RETURN
4274 END SUBROUTINE dvthermcell
4275 ✗ SUBROUTINE dqthermcell2(ngrid, nlay, ptimestep, fm, entr, masse, frac, q, dq, &
4276 ✗ qa)
4277 USE dimphy
4278 IMPLICIT NONE
4279
4280 ! =======================================================================
4281
4282 ! Calcul du transport verticale dans la couche limite en presence
4283 ! de "thermiques" explicitement representes
4284 ! calcul du dq/dt une fois qu'on connait les ascendances
4285
4286 ! =======================================================================
4287
4288 INTEGER ngrid, nlay
4289
4290 REAL ptimestep
4291 REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4292 REAL entr(ngrid, nlay), frac(ngrid, nlay)
4293 REAL q(ngrid, nlay)
4294 REAL dq(ngrid, nlay)
4295
4296 ✗ REAL qa(klon, klev), detr(klon, klev), wqd(klon, klev+1)
4297 ✗ REAL qe(klon, klev), zf, zf2
4298
4299 INTEGER ig, k
4300
4301 ! calcul du detrainement
4302
4303 ✗ DO k = 1, nlay
4304 ✗ DO ig = 1, ngrid
4305 ✗ detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4306 END DO
4307 END DO
4308
4309 ! calcul de la valeur dans les ascendances
4310 ✗ DO ig = 1, ngrid
4311 ✗ qa(ig, 1) = q(ig, 1)
4312 ✗ qe(ig, 1) = q(ig, 1)
4313 END DO
4314
4315 ✗ DO k = 2, nlay
4316 ✗ DO ig = 1, ngrid
4317 ✗ IF ((fm(ig,k+1)+detr(ig,k))*ptimestep>1.E-5*masse(ig,k)) THEN
4318 ✗ zf = 0.5*(frac(ig,k)+frac(ig,k+1))
4319 ✗ zf2 = 1./(1.-zf)
4320 qa(ig, k) = (fm(ig,k)*qa(ig,k-1)+zf2*entr(ig,k)*q(ig,k))/ &
4321 ✗ (fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2)
4322 ✗ qe(ig, k) = (q(ig,k)-zf*qa(ig,k))*zf2
4323 ELSE
4324 ✗ qa(ig, k) = q(ig, k)
4325 ✗ qe(ig, k) = q(ig, k)
4326 END IF
4327 END DO
4328 END DO
4329
4330 ✗ DO k = 2, nlay
4331 ✗ DO ig = 1, ngrid
4332 ! wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
4333 ✗ wqd(ig, k) = fm(ig, k)*qe(ig, k)
4334 END DO
4335 END DO
4336 ✗ DO ig = 1, ngrid
4337 ✗ wqd(ig, 1) = 0.
4338 ✗ wqd(ig, nlay+1) = 0.
4339 END DO
4340
4341 ✗ DO k = 1, nlay
4342 ✗ DO ig = 1, ngrid
4343 dq(ig, k) = (detr(ig,k)*qa(ig,k)-entr(ig,k)*qe(ig,k)-wqd(ig,k)+wqd(ig,k &
4344 ✗ +1))/masse(ig, k)
4345 END DO
4346 END DO
4347
4348 ✗ RETURN
4349 END SUBROUTINE dqthermcell2
4350 ✗ SUBROUTINE dvthermcell2(ngrid, nlay, ptimestep, fm, entr, masse, fraca, &
4351 larga, u, v, du, dv, ua, va)
4352 USE dimphy
4353 IMPLICIT NONE
4354
4355 ! =======================================================================
4356
4357 ! Calcul du transport verticale dans la couche limite en presence
4358 ! de "thermiques" explicitement representes
4359 ! calcul du dq/dt une fois qu'on connait les ascendances
4360
4361 ! =======================================================================
4362
4363 INTEGER ngrid, nlay
4364
4365 REAL ptimestep
4366 REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4367 REAL fraca(ngrid, nlay+1)
4368 REAL larga(ngrid)
4369 REAL entr(ngrid, nlay)
4370 REAL u(ngrid, nlay)
4371 REAL ua(ngrid, nlay)
4372 REAL du(ngrid, nlay)
4373 REAL v(ngrid, nlay)
4374 REAL va(ngrid, nlay)
4375 REAL dv(ngrid, nlay)
4376
4377 ✗ REAL qa(klon, klev), detr(klon, klev), zf, zf2
4378 ✗ REAL wvd(klon, klev+1), wud(klon, klev+1)
4379 ✗ REAL gamma0, gamma(klon, klev+1)
4380 ✗ REAL ue(klon, klev), ve(klon, klev)
4381 REAL dua, dva
4382 INTEGER iter
4383
4384 INTEGER ig, k
4385
4386 ! calcul du detrainement
4387
4388 ✗ DO k = 1, nlay
4389 ✗ DO ig = 1, ngrid
4390 ✗ detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4391 END DO
4392 END DO
4393
4394 ! calcul de la valeur dans les ascendances
4395 ✗ DO ig = 1, ngrid
4396 ✗ ua(ig, 1) = u(ig, 1)
4397 ✗ va(ig, 1) = v(ig, 1)
4398 ✗ ue(ig, 1) = u(ig, 1)
4399 ✗ ve(ig, 1) = v(ig, 1)
4400 END DO
4401
4402 ✗ DO k = 2, nlay
4403 ✗ DO ig = 1, ngrid
4404 ✗ IF ((fm(ig,k+1)+detr(ig,k))*ptimestep>1.E-5*masse(ig,k)) THEN
4405 ! On it�re sur la valeur du coeff de freinage.
4406 ! gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
4407 gamma0 = masse(ig, k)*sqrt(0.5*(fraca(ig,k+1)+fraca(ig, &
4408 ✗ k)))*0.5/larga(ig)*1.
4409 ! s *0.5
4410 ! gamma0=0.
4411 zf = 0.5*(fraca(ig,k)+fraca(ig,k+1))
4412 zf = 0.
4413 zf2 = 1./(1.-zf)
4414 ! la premi�re fois on multiplie le coefficient de freinage
4415 ! par le module du vent dans la couche en dessous.
4416 ✗ dua = ua(ig, k-1) - u(ig, k-1)
4417 ✗ dva = va(ig, k-1) - v(ig, k-1)
4418 ✗ DO iter = 1, 5
4419 ! On choisit une relaxation lineaire.
4420 gamma(ig, k) = gamma0
4421 ! On choisit une relaxation quadratique.
4422 ✗ gamma(ig, k) = gamma0*sqrt(dua**2+dva**2)
4423 ua(ig, k) = (fm(ig,k)*ua(ig,k-1)+(zf2*entr(ig,k)+gamma(ig, &
4424 k))*u(ig,k))/(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2+gamma(ig,k) &
4425 ✗ )
4426 va(ig, k) = (fm(ig,k)*va(ig,k-1)+(zf2*entr(ig,k)+gamma(ig, &
4427 k))*v(ig,k))/(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2+gamma(ig,k) &
4428 ✗ )
4429 ! print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
4430 ✗ dua = ua(ig, k) - u(ig, k)
4431 ✗ dva = va(ig, k) - v(ig, k)
4432 ✗ ue(ig, k) = (u(ig,k)-zf*ua(ig,k))*zf2
4433 ✗ ve(ig, k) = (v(ig,k)-zf*va(ig,k))*zf2
4434 END DO
4435 ELSE
4436 ✗ ua(ig, k) = u(ig, k)
4437 ✗ va(ig, k) = v(ig, k)
4438 ✗ ue(ig, k) = u(ig, k)
4439 ✗ ve(ig, k) = v(ig, k)
4440 ✗ gamma(ig, k) = 0.
4441 END IF
4442 END DO
4443 END DO
4444
4445 ✗ DO k = 2, nlay
4446 ✗ DO ig = 1, ngrid
4447 ✗ wud(ig, k) = fm(ig, k)*ue(ig, k)
4448 ✗ wvd(ig, k) = fm(ig, k)*ve(ig, k)
4449 END DO
4450 END DO
4451 ✗ DO ig = 1, ngrid
4452 ✗ wud(ig, 1) = 0.
4453 ✗ wud(ig, nlay+1) = 0.
4454 ✗ wvd(ig, 1) = 0.
4455 ✗ wvd(ig, nlay+1) = 0.
4456 END DO
4457
4458 ✗ DO k = 1, nlay
4459 ✗ DO ig = 1, ngrid
4460 du(ig, k) = ((detr(ig,k)+gamma(ig,k))*ua(ig,k)-(entr(ig,k)+gamma(ig, &
4461 ✗ k))*ue(ig,k)-wud(ig,k)+wud(ig,k+1))/masse(ig, k)
4462 dv(ig, k) = ((detr(ig,k)+gamma(ig,k))*va(ig,k)-(entr(ig,k)+gamma(ig, &
4463 ✗ k))*ve(ig,k)-wvd(ig,k)+wvd(ig,k+1))/masse(ig, k)
4464 END DO
4465 END DO
4466
4467 ✗ RETURN
4468 END SUBROUTINE dvthermcell2
4469 ✗ SUBROUTINE thermcell_sec(ngrid, nlay, ptimestep, pplay, pplev, pphi, zlev, &
4470 ✗ pu, pv, pt, po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0 & ! s
4471 ! ,pu_therm,pv_therm
4472 , r_aspect, l_mix, w2di, tho)
4473
4474 USE dimphy
4475 IMPLICIT NONE
4476
4477 ! =======================================================================
4478
4479 ! Calcul du transport verticale dans la couche limite en presence
4480 ! de "thermiques" explicitement representes
4481
4482 ! R��criture � partir d'un listing papier � Habas, le 14/02/00
4483
4484 ! le thermique est suppos� homog�ne et dissip� par m�lange avec
4485 ! son environnement. la longueur l_mix contr�le l'efficacit� du
4486 ! m�lange
4487
4488 ! Le calcul du transport des diff�rentes esp�ces se fait en prenant
4489 ! en compte:
4490 ! 1. un flux de masse montant
4491 ! 2. un flux de masse descendant
4492 ! 3. un entrainement
4493 ! 4. un detrainement
4494
4495 ! =======================================================================
4496
4497 ! -----------------------------------------------------------------------
4498 ! declarations:
4499 ! -------------
4500
4501 include "YOMCST.h"
4502
4503 ! arguments:
4504 ! ----------
4505
4506 INTEGER ngrid, nlay, w2di
4507 REAL tho
4508 REAL ptimestep, l_mix, r_aspect
4509 REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
4510 REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
4511 REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
4512 REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
4513 REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
4514 REAL pphi(ngrid, nlay)
4515
4516 INTEGER idetr
4517 SAVE idetr
4518 DATA idetr/3/
4519 !$OMP THREADPRIVATE(idetr)
4520
4521 ! local:
4522 ! ------
4523
4524 ✗ INTEGER ig, k, l, lmaxa(klon), lmix(klon)
4525 REAL zsortie1d(klon)
4526 ! CR: on remplace lmax(klon,klev+1)
4527 ✗ INTEGER lmax(klon), lmin(klon), lentr(klon)
4528 ✗ REAL linter(klon)
4529 ✗ REAL zmix(klon), fracazmix(klon)
4530 ! RC
4531 ✗ REAL zmax(klon), zw, zz, zw2(klon, klev+1), ztva(klon, klev), zzz
4532
4533 ✗ REAL zlev(klon, klev+1), zlay(klon, klev)
4534 ✗ REAL zh(klon, klev), zdhadj(klon, klev)
4535 ✗ REAL ztv(klon, klev)
4536 ✗ REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
4537 REAL wh(klon, klev+1)
4538 REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
4539 REAL zla(klon, klev+1)
4540 REAL zwa(klon, klev+1)
4541 REAL zld(klon, klev+1)
4542 REAL zwd(klon, klev+1)
4543 REAL zsortie(klon, klev)
4544 ✗ REAL zva(klon, klev)
4545 ✗ REAL zua(klon, klev)
4546 ✗ REAL zoa(klon, klev)
4547
4548 ✗ REAL zha(klon, klev)
4549 ✗ REAL wa_moy(klon, klev+1)
4550 ✗ REAL fraca(klon, klev+1)
4551 ✗ REAL fracc(klon, klev+1)
4552 REAL zf, zf2
4553 ✗ REAL thetath2(klon, klev), wth2(klon, klev)
4554 ! common/comtherm/thetath2,wth2
4555
4556 REAL count_time
4557 INTEGER ialt
4558
4559 LOGICAL sorties
4560 ✗ REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
4561 ✗ REAL zpspsk(klon, klev)
4562
4563 ! real wmax(klon,klev),wmaxa(klon)
4564 ✗ REAL wmax(klon), wmaxa(klon)
4565 ✗ REAL wa(klon, klev, klev+1)
4566 ✗ REAL wd(klon, klev+1)
4567 REAL larg_part(klon, klev, klev+1)
4568 ✗ REAL fracd(klon, klev+1)
4569 ✗ REAL xxx(klon, klev+1)
4570 ✗ REAL larg_cons(klon, klev+1)
4571 ✗ REAL larg_detr(klon, klev+1)
4572 ✗ REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
4573 REAL pu_therm(klon, klev), pv_therm(klon, klev)
4574 ✗ REAL fm(klon, klev+1), entr(klon, klev)
4575 ✗ REAL fmc(klon, klev+1)
4576
4577 ! CR:nouvelles variables
4578 ✗ REAL f_star(klon, klev+1), entr_star(klon, klev)
4579 ✗ REAL entr_star_tot(klon), entr_star2(klon)
4580 ✗ REAL f(klon), f0(klon)
4581 ✗ REAL zlevinter(klon)
4582 LOGICAL first
4583 DATA first/.FALSE./
4584 SAVE first
4585 !$OMP THREADPRIVATE(first)
4586 ! RC
4587
4588 CHARACTER *2 str2
4589 CHARACTER *10 str10
4590
4591 CHARACTER (LEN=20) :: modname = 'thermcell_sec'
4592 CHARACTER (LEN=80) :: abort_message
4593
4594 LOGICAL vtest(klon), down
4595
4596 EXTERNAL scopy
4597
4598 INTEGER ncorrec, ll
4599 SAVE ncorrec
4600 DATA ncorrec/0/
4601 !$OMP THREADPRIVATE(ncorrec)
4602
4603
4604 ! -----------------------------------------------------------------------
4605 ! initialisation:
4606 ! ---------------
4607
4608 sorties = .TRUE.
4609 ✗ IF (ngrid/=klon) THEN
4610 ✗ PRINT *
4611 ✗ PRINT *, 'STOP dans convadj'
4612 ✗ PRINT *, 'ngrid =', ngrid
4613 ✗ PRINT *, 'klon =', klon
4614 END IF
4615
4616 ! -----------------------------------------------------------------------
4617 ! incrementation eventuelle de tendances precedentes:
4618 ! ---------------------------------------------------
4619
4620 ! print*,'0 OK convect8'
4621
4622 ✗ DO l = 1, nlay
4623 ✗ DO ig = 1, ngrid
4624 ✗ zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa
4625 ✗ zh(ig, l) = pt(ig, l)/zpspsk(ig, l)
4626 ✗ zu(ig, l) = pu(ig, l)
4627 ✗ zv(ig, l) = pv(ig, l)
4628 ✗ zo(ig, l) = po(ig, l)
4629 ✗ ztv(ig, l) = zh(ig, l)*(1.+0.61*zo(ig,l))
4630 END DO
4631 END DO
4632
4633 ! print*,'1 OK convect8'
4634 ! --------------------
4635
4636
4637 ! + + + + + + + + + + +
4638
4639
4640 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
4641 ! wh,wt,wo ...
4642
4643 ! + + + + + + + + + + + zh,zu,zv,zo,rho
4644
4645
4646 ! -------------------- zlev(1)
4647 ! \\\\\\\\\\\\\\\\\\\
4648
4649
4650 ! -----------------------------------------------------------------------
4651 ! Calcul des altitudes des couches
4652 ! -----------------------------------------------------------------------
4653
4654 ✗ DO l = 2, nlay
4655 ✗ DO ig = 1, ngrid
4656 ✗ zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
4657 END DO
4658 END DO
4659 ✗ DO ig = 1, ngrid
4660 ✗ zlev(ig, 1) = 0.
4661 ✗ zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
4662 END DO
4663 ✗ DO l = 1, nlay
4664 ✗ DO ig = 1, ngrid
4665 ✗ zlay(ig, l) = pphi(ig, l)/rg
4666 END DO
4667 END DO
4668
4669 ! print*,'2 OK convect8'
4670 ! -----------------------------------------------------------------------
4671 ! Calcul des densites
4672 ! -----------------------------------------------------------------------
4673
4674 ✗ DO l = 1, nlay
4675 ✗ DO ig = 1, ngrid
4676 ✗ rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*zh(ig,l))
4677 END DO
4678 END DO
4679
4680 ✗ DO l = 2, nlay
4681 ✗ DO ig = 1, ngrid
4682 ✗ rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
4683 END DO
4684 END DO
4685
4686 ✗ DO k = 1, nlay
4687 ✗ DO l = 1, nlay + 1
4688 ✗ DO ig = 1, ngrid
4689 ✗ wa(ig, k, l) = 0.
4690 END DO
4691 END DO
4692 END DO
4693
4694 ! print*,'3 OK convect8'
4695 ! ------------------------------------------------------------------
4696 ! Calcul de w2, quarre de w a partir de la cape
4697 ! a partir de w2, on calcule wa, vitesse de l'ascendance
4698
4699 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
4700 ! w2 est stoke dans wa
4701
4702 ! ATTENTION: dans convect8, on n'utilise le calcule des wa
4703 ! independants par couches que pour calculer l'entrainement
4704 ! a la base et la hauteur max de l'ascendance.
4705
4706 ! Indicages:
4707 ! l'ascendance provenant du niveau k traverse l'interface l avec
4708 ! une vitesse wa(k,l).
4709
4710 ! --------------------
4711
4712 ! + + + + + + + + + +
4713
4714 ! wa(k,l) ---- -------------------- l
4715 ! / ! /||\ + + + + + + + + + +
4716 ! ||
4717 ! || --------------------
4718 ! ||
4719 ! || + + + + + + + + + +
4720 ! ||
4721 ! || --------------------
4722 ! ||__
4723 ! |___ + + + + + + + + + + k
4724
4725 ! --------------------
4726
4727
4728
4729 ! ------------------------------------------------------------------
4730
4731 ! CR: ponderation entrainement des couches instables
4732 ! def des entr_star tels que entr=f*entr_star
4733 ✗ DO l = 1, klev
4734 ✗ DO ig = 1, ngrid
4735 ✗ entr_star(ig, l) = 0.
4736 END DO
4737 END DO
4738 ! determination de la longueur de la couche d entrainement
4739 ✗ DO ig = 1, ngrid
4740 ✗ lentr(ig) = 1
4741 END DO
4742
4743 ! on ne considere que les premieres couches instables
4744 ✗ DO k = nlay - 2, 1, -1
4745 ✗ DO ig = 1, ngrid
4746 ✗ IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN
4747 ✗ lentr(ig) = k
4748 END IF
4749 END DO
4750 END DO
4751
4752 ! determination du lmin: couche d ou provient le thermique
4753 ✗ DO ig = 1, ngrid
4754 ✗ lmin(ig) = 1
4755 END DO
4756 ✗ DO ig = 1, ngrid
4757 ✗ DO l = nlay, 2, -1
4758 ✗ IF (ztv(ig,l-1)>ztv(ig,l)) THEN
4759 ✗ lmin(ig) = l - 1
4760 END IF
4761 END DO
4762 END DO
4763
4764 ! definition de l'entrainement des couches
4765 ✗ DO l = 1, klev - 1
4766 ✗ DO ig = 1, ngrid
4767 ✗ IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<=lentr(ig)) THEN
4768 entr_star(ig, l) = (ztv(ig,l)-ztv(ig,l+1))** & ! s
4769 ! (zlev(ig,l+1)-zlev(ig,l))
4770 ✗ sqrt(zlev(ig,l+1))
4771 END IF
4772 END DO
4773 END DO
4774 ! pas de thermique si couche 1 stable
4775 ✗ DO ig = 1, ngrid
4776 ✗ IF (lmin(ig)>1) THEN
4777 ✗ DO l = 1, klev
4778 ✗ entr_star(ig, l) = 0.
4779 END DO
4780 END IF
4781 END DO
4782 ! calcul de l entrainement total
4783 ✗ DO ig = 1, ngrid
4784 ✗ entr_star_tot(ig) = 0.
4785 END DO
4786 ✗ DO ig = 1, ngrid
4787 ✗ DO k = 1, klev
4788 ✗ entr_star_tot(ig) = entr_star_tot(ig) + entr_star(ig, k)
4789 END DO
4790 END DO
4791
4792 ! print*,'fin calcul entr_star'
4793 ✗ DO k = 1, klev
4794 ✗ DO ig = 1, ngrid
4795 ✗ ztva(ig, k) = ztv(ig, k)
4796 END DO
4797 END DO
4798 ! RC
4799 ! print*,'7 OK convect8'
4800 ✗ DO k = 1, klev + 1
4801 ✗ DO ig = 1, ngrid
4802 ✗ zw2(ig, k) = 0.
4803 ✗ fmc(ig, k) = 0.
4804 ! CR
4805 ✗ f_star(ig, k) = 0.
4806 ! RC
4807 ✗ larg_cons(ig, k) = 0.
4808 ✗ larg_detr(ig, k) = 0.
4809 ✗ wa_moy(ig, k) = 0.
4810 END DO
4811 END DO
4812
4813 ! print*,'8 OK convect8'
4814 ✗ DO ig = 1, ngrid
4815 ✗ linter(ig) = 1.
4816 ✗ lmaxa(ig) = 1
4817 ✗ lmix(ig) = 1
4818 ✗ wmaxa(ig) = 0.
4819 END DO
4820
4821 ! CR:
4822 ✗ DO l = 1, nlay - 2
4823 ✗ DO ig = 1, ngrid
4824 ✗ IF (ztv(ig,l)>ztv(ig,l+1) .AND. entr_star(ig,l)>1.E-10 .AND. &
4825 zw2(ig,l)<1E-10) THEN
4826 ✗ f_star(ig, l+1) = entr_star(ig, l)
4827 ! test:calcul de dteta
4828 zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
4829 ✗ (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
4830 ✗ larg_detr(ig, l) = 0.
4831 ✗ ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+entr_star(ig, &
4832 l)>1.E-10)) THEN
4833 ✗ f_star(ig, l+1) = f_star(ig, l) + entr_star(ig, l)
4834 ztva(ig, l) = (f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)*ztv(ig,l))/ &
4835 ✗ f_star(ig, l+1)
4836 zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/f_star(ig,l+1))**2 + &
4837 ✗ 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
4838 END IF
4839 ! determination de zmax continu par interpolation lineaire
4840 ✗ IF (zw2(ig,l+1)<0.) THEN
4841 ! test
4842 IF (abs(zw2(ig,l+1)-zw2(ig,l))<1E-10) THEN
4843 ! print*,'pb linter'
4844 END IF
4845 linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
4846 ✗ ig,l))
4847 ✗ zw2(ig, l+1) = 0.
4848 ✗ lmaxa(ig) = l
4849 ELSE
4850 IF (zw2(ig,l+1)<0.) THEN
4851 ! print*,'pb1 zw2<0'
4852 END IF
4853 ✗ wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
4854 END IF
4855 ✗ IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
4856 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
4857 ✗ lmix(ig) = l + 1
4858 ✗ wmaxa(ig) = wa_moy(ig, l+1)
4859 END IF
4860 END DO
4861 END DO
4862 ! print*,'fin calcul zw2'
4863
4864 ! Calcul de la couche correspondant a la hauteur du thermique
4865 ✗ DO ig = 1, ngrid
4866 ✗ lmax(ig) = lentr(ig)
4867 END DO
4868 ✗ DO ig = 1, ngrid
4869 ✗ DO l = nlay, lentr(ig) + 1, -1
4870 ✗ IF (zw2(ig,l)<=1.E-10) THEN
4871 ✗ lmax(ig) = l - 1
4872 END IF
4873 END DO
4874 END DO
4875 ! pas de thermique si couche 1 stable
4876 ✗ DO ig = 1, ngrid
4877 ✗ IF (lmin(ig)>1) THEN
4878 ✗ lmax(ig) = 1
4879 ✗ lmin(ig) = 1
4880 END IF
4881 END DO
4882
4883 ! Determination de zw2 max
4884 ✗ DO ig = 1, ngrid
4885 ✗ wmax(ig) = 0.
4886 END DO
4887
4888 ✗ DO l = 1, nlay
4889 ✗ DO ig = 1, ngrid
4890 ✗ IF (l<=lmax(ig)) THEN
4891 IF (zw2(ig,l)<0.) THEN
4892 ! print*,'pb2 zw2<0'
4893 END IF
4894 ✗ zw2(ig, l) = sqrt(zw2(ig,l))
4895 ✗ wmax(ig) = max(wmax(ig), zw2(ig,l))
4896 ELSE
4897 ✗ zw2(ig, l) = 0.
4898 END IF
4899 END DO
4900 END DO
4901
4902 ! Longueur caracteristique correspondant a la hauteur des thermiques.
4903 ✗ DO ig = 1, ngrid
4904 ✗ zmax(ig) = 0.
4905 ✗ zlevinter(ig) = zlev(ig, 1)
4906 END DO
4907 ✗ DO ig = 1, ngrid
4908 ! calcul de zlevinter
4909 zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
4910 ✗ zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
4911 ✗ zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,lmin(ig)))
4912 END DO
4913
4914 ! print*,'avant fermeture'
4915 ! Fermeture,determination de f
4916 ✗ DO ig = 1, ngrid
4917 ✗ entr_star2(ig) = 0.
4918 END DO
4919 ✗ DO ig = 1, ngrid
4920 ✗ IF (entr_star_tot(ig)<1.E-10) THEN
4921 ✗ f(ig) = 0.
4922 ELSE
4923 ✗ DO k = lmin(ig), lentr(ig)
4924 entr_star2(ig) = entr_star2(ig) + entr_star(ig, k)**2/(rho(ig,k)*( &
4925 ✗ zlev(ig,k+1)-zlev(ig,k)))
4926 END DO
4927 ! Nouvelle fermeture
4928 f(ig) = wmax(ig)/(max(500.,zmax(ig))*r_aspect*entr_star2(ig))* &
4929 ✗ entr_star_tot(ig)
4930 ! test
4931 ! if (first) then
4932 ! f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
4933 ! s *wmax(ig))
4934 ! endif
4935 END IF
4936 ! f0(ig)=f(ig)
4937 ! first=.true.
4938 END DO
4939 ! print*,'apres fermeture'
4940
4941 ! Calcul de l'entrainement
4942 ✗ DO k = 1, klev
4943 ✗ DO ig = 1, ngrid
4944 ✗ entr(ig, k) = f(ig)*entr_star(ig, k)
4945 END DO
4946 END DO
4947 ! CR:test pour entrainer moins que la masse
4948 ✗ DO ig = 1, ngrid
4949 ✗ DO l = 1, lentr(ig)
4950 ✗ IF ((entr(ig,l)*ptimestep)>(0.9*masse(ig,l))) THEN
4951 entr(ig, l+1) = entr(ig, l+1) + entr(ig, l) - &
4952 ✗ 0.9*masse(ig, l)/ptimestep
4953 ✗ entr(ig, l) = 0.9*masse(ig, l)/ptimestep
4954 END IF
4955 END DO
4956 END DO
4957 ! CR: fin test
4958 ! Calcul des flux
4959 ✗ DO ig = 1, ngrid
4960 ✗ DO l = 1, lmax(ig) - 1
4961 ✗ fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
4962 END DO
4963 END DO
4964
4965 ! RC
4966
4967
4968 ! print*,'9 OK convect8'
4969 ! print*,'WA1 ',wa_moy
4970
4971 ! determination de l'indice du debut de la mixed layer ou w decroit
4972
4973 ! calcul de la largeur de chaque ascendance dans le cas conservatif.
4974 ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
4975 ! d'une couche est �gale � la hauteur de la couche alimentante.
4976 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
4977 ! de la vitesse d'entrainement horizontal dans la couche alimentante.
4978
4979 ✗ DO l = 2, nlay
4980 ✗ DO ig = 1, ngrid
4981 ✗ IF (l<=lmaxa(ig)) THEN
4982 ✗ zw = max(wa_moy(ig,l), 1.E-10)
4983 ✗ larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw)
4984 END IF
4985 END DO
4986 END DO
4987
4988 ✗ DO l = 2, nlay
4989 ✗ DO ig = 1, ngrid
4990 ✗ IF (l<=lmaxa(ig)) THEN
4991 ! if (idetr.eq.0) then
4992 ! cette option est finalement en dur.
4993 IF ((l_mix*zlev(ig,l))<0.) THEN
4994 ! print*,'pb l_mix*zlev<0'
4995 END IF
4996 ! CR: test: nouvelle def de lambda
4997 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
4998 ✗ IF (zw2(ig,l)>1.E-10) THEN
4999 ✗ larg_detr(ig, l) = sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
5000 ELSE
5001 ✗ larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
5002 END IF
5003 ! RC
5004 ! else if (idetr.eq.1) then
5005 ! larg_detr(ig,l)=larg_cons(ig,l)
5006 ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
5007 ! else if (idetr.eq.2) then
5008 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5009 ! s *sqrt(wa_moy(ig,l))
5010 ! else if (idetr.eq.4) then
5011 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5012 ! s *wa_moy(ig,l)
5013 ! endif
5014 END IF
5015 END DO
5016 END DO
5017
5018 ! print*,'10 OK convect8'
5019 ! print*,'WA2 ',wa_moy
5020 ! calcul de la fraction de la maille concern�e par l'ascendance en tenant
5021 ! compte de l'epluchage du thermique.
5022
5023 ! CR def de zmix continu (profil parabolique des vitesses)
5024 ✗ DO ig = 1, ngrid
5025 ✗ IF (lmix(ig)>1.) THEN
5026 ! test
5027 ✗ IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
5028 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
5029 zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- &
5030 (zlev(ig,lmix(ig)))))>1E-10) THEN
5031
5032 zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) &
5033 )**2-(zlev(ig,lmix(ig)+1))**2)-(zw2(ig,lmix(ig))-zw2(ig, &
5034 lmix(ig)+1))*((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))/ &
5035 (2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
5036 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
5037 ✗ zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
5038 ELSE
5039 ✗ zmix(ig) = zlev(ig, lmix(ig))
5040 ! print*,'pb zmix'
5041 END IF
5042 ELSE
5043 ✗ zmix(ig) = 0.
5044 END IF
5045 ! test
5046 ✗ IF ((zmax(ig)-zmix(ig))<0.) THEN
5047 ✗ zmix(ig) = 0.99*zmax(ig)
5048 ! print*,'pb zmix>zmax'
5049 END IF
5050 END DO
5051
5052 ! calcul du nouveau lmix correspondant
5053 ✗ DO ig = 1, ngrid
5054 ✗ DO l = 1, klev
5055 ✗ IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
5056 ✗ lmix(ig) = l
5057 END IF
5058 END DO
5059 END DO
5060
5061 ✗ DO l = 2, nlay
5062 ✗ DO ig = 1, ngrid
5063 ✗ IF (larg_cons(ig,l)>1.) THEN
5064 ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
5065 ✗ fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
5066 ! test
5067 ✗ fraca(ig, l) = max(fraca(ig,l), 0.)
5068 ✗ fraca(ig, l) = min(fraca(ig,l), 0.5)
5069 ✗ fracd(ig, l) = 1. - fraca(ig, l)
5070 ✗ fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
5071 ELSE
5072 ! wa_moy(ig,l)=0.
5073 ✗ fraca(ig, l) = 0.
5074 ✗ fracc(ig, l) = 0.
5075 ✗ fracd(ig, l) = 1.
5076 END IF
5077 END DO
5078 END DO
5079 ! CR: calcul de fracazmix
5080 ✗ DO ig = 1, ngrid
5081 fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
5082 (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
5083 fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca(ig &
5084 ✗ ,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
5085 END DO
5086
5087 ✗ DO l = 2, nlay
5088 ✗ DO ig = 1, ngrid
5089 ✗ IF (larg_cons(ig,l)>1.) THEN
5090 ✗ IF (l>lmix(ig)) THEN
5091 ! test
5092 ✗ IF (zmax(ig)-zmix(ig)<1.E-10) THEN
5093 ! print*,'pb xxx'
5094 ✗ xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
5095 ELSE
5096 ✗ xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
5097 END IF
5098 ✗ IF (idetr==0) THEN
5099 ✗ fraca(ig, l) = fracazmix(ig)
5100 ✗ ELSE IF (idetr==1) THEN
5101 ✗ fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
5102 ✗ ELSE IF (idetr==2) THEN
5103 ✗ fraca(ig, l) = fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
5104 ELSE
5105 ✗ fraca(ig, l) = fracazmix(ig)*xxx(ig, l)**2
5106 END IF
5107 ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
5108 ✗ fraca(ig, l) = max(fraca(ig,l), 0.)
5109 ✗ fraca(ig, l) = min(fraca(ig,l), 0.5)
5110 ✗ fracd(ig, l) = 1. - fraca(ig, l)
5111 ✗ fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
5112 END IF
5113 END IF
5114 END DO
5115 END DO
5116
5117 ! print*,'fin calcul fraca'
5118 ! print*,'11 OK convect8'
5119 ! print*,'Ea3 ',wa_moy
5120 ! ------------------------------------------------------------------
5121 ! Calcul de fracd, wd
5122 ! somme wa - wd = 0
5123 ! ------------------------------------------------------------------
5124
5125
5126 ✗ DO ig = 1, ngrid
5127 ✗ fm(ig, 1) = 0.
5128 ✗ fm(ig, nlay+1) = 0.
5129 END DO
5130
5131 ✗ DO l = 2, nlay
5132 ✗ DO ig = 1, ngrid
5133 ✗ fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
5134 ! CR:test
5135 ✗ IF (entr(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) THEN
5136 ✗ fm(ig, l) = fm(ig, l-1)
5137 ! write(1,*)'ajustement fm, l',l
5138 END IF
5139 ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
5140 ! RC
5141 END DO
5142 ✗ DO ig = 1, ngrid
5143 ✗ IF (fracd(ig,l)<0.1) THEN
5144 ✗ abort_message = 'fracd trop petit'
5145 ✗ CALL abort_physic(modname, abort_message, 1)
5146 ELSE
5147 ! vitesse descendante "diagnostique"
5148 ✗ wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
5149 END IF
5150 END DO
5151 END DO
5152
5153 ✗ DO l = 1, nlay
5154 ✗ DO ig = 1, ngrid
5155 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
5156 ✗ masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
5157 END DO
5158 END DO
5159
5160 ! print*,'12 OK convect8'
5161 ! print*,'WA4 ',wa_moy
5162 ! c------------------------------------------------------------------
5163 ! calcul du transport vertical
5164 ! ------------------------------------------------------------------
5165
5166 GO TO 4444
5167 ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
5168 DO l = 2, nlay - 1
5169 DO ig = 1, ngrid
5170 IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( &
5171 ig,l+1)) THEN
5172 ! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
5173 ! s ,fm(ig,l+1)*ptimestep
5174 ! s ,' M=',masse(ig,l),masse(ig,l+1)
5175 END IF
5176 END DO
5177 END DO
5178
5179 DO l = 1, nlay
5180 DO ig = 1, ngrid
5181 IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
5182 ! print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
5183 ! s ,entr(ig,l)*ptimestep
5184 ! s ,' M=',masse(ig,l)
5185 END IF
5186 END DO
5187 END DO
5188
5189 DO l = 1, nlay
5190 DO ig = 1, ngrid
5191 IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
5192 ! print*,'WARN!!! fm exagere ig=',ig,' l=',l
5193 ! s ,' FM=',fm(ig,l)
5194 END IF
5195 IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
5196 ! print*,'WARN!!! masse exagere ig=',ig,' l=',l
5197 ! s ,' M=',masse(ig,l)
5198 ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
5199 ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
5200 ! print*,'zlev(ig,l+1),zlev(ig,l)'
5201 ! s ,zlev(ig,l+1),zlev(ig,l)
5202 ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
5203 ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
5204 END IF
5205 IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
5206 ! print*,'WARN!!! entr exagere ig=',ig,' l=',l
5207 ! s ,' E=',entr(ig,l)
5208 END IF
5209 END DO
5210 END DO
5211
5212 4444 CONTINUE
5213
5214 ! CR:redefinition du entr
5215 ✗ DO l = 1, nlay
5216 ✗ DO ig = 1, ngrid
5217 ✗ detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1)
5218 ✗ IF (detr(ig,l)<0.) THEN
5219 ✗ entr(ig, l) = entr(ig, l) - detr(ig, l)
5220 ✗ detr(ig, l) = 0.
5221 ! print*,'WARNING !!! detrainement negatif ',ig,l
5222 END IF
5223 END DO
5224 END DO
5225 ! RC
5226 ✗ IF (w2di==1) THEN
5227 ✗ fm0 = fm0 + ptimestep*(fm-fm0)/tho
5228 ✗ entr0 = entr0 + ptimestep*(entr-entr0)/tho
5229 ELSE
5230 ✗ fm0 = fm
5231 ✗ entr0 = entr
5232 END IF
5233
5234 IF (1==1) THEN
5235 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, &
5236 ✗ zha)
5237 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, &
5238 ✗ zoa)
5239 ELSE
5240 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
5241 zdhadj, zha)
5242 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
5243 pdoadj, zoa)
5244 END IF
5245
5246 IF (1==0) THEN
5247 CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
5248 zu, zv, pduadj, pdvadj, zua, zva)
5249 ELSE
5250 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
5251 ✗ zua)
5252 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
5253 ✗ zva)
5254 END IF
5255
5256 ✗ DO l = 1, nlay
5257 ✗ DO ig = 1, ngrid
5258 ✗ zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
5259 ✗ zf2 = zf/(1.-zf)
5260 ✗ thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l))**2
5261 ✗ wth2(ig, l) = zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
5262 END DO
5263 END DO
5264
5265
5266
5267 ! print*,'13 OK convect8'
5268 ! print*,'WA5 ',wa_moy
5269 ✗ DO l = 1, nlay
5270 ✗ DO ig = 1, ngrid
5271 ✗ pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l)
5272 END DO
5273 END DO
5274
5275
5276 ! do l=1,nlay
5277 ! do ig=1,ngrid
5278 ! if(abs(pdtadj(ig,l))*86400..gt.500.) then
5279 ! print*,'WARN!!! ig=',ig,' l=',l
5280 ! s ,' pdtadj=',pdtadj(ig,l)
5281 ! endif
5282 ! if(abs(pdoadj(ig,l))*86400..gt.1.) then
5283 ! print*,'WARN!!! ig=',ig,' l=',l
5284 ! s ,' pdoadj=',pdoadj(ig,l)
5285 ! endif
5286 ! enddo
5287 ! enddo
5288
5289 ! print*,'14 OK convect8'
5290 ! ------------------------------------------------------------------
5291 ! Calculs pour les sorties
5292 ! ------------------------------------------------------------------
5293
5294 ✗ RETURN
5295 END SUBROUTINE thermcell_sec
5296
5297