GCC Code Coverage Report


Directory: ./
File: phys/cv3p2_closure.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 327 0.0%
Branches: 0 422 0.0%

Line Branch Exec Source
1
2
3 SUBROUTINE cv3p2_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
4 tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, &
5 iflag, coef, plim1, plim2, asupmax, supmax0, asupmaxmin, cbmflast, plfc, &
6 wbeff)
7
8
9 ! **************************************************************
10 ! *
11 ! CV3P2_CLOSURE *
12 ! Ale & Alp Closure of Convect3 *
13 ! *
14 ! written by : Kerry Emanuel *
15 ! vectorization: S. Bony *
16 ! modified by : Jean-Yves Grandpeix, 18/06/2003, 19.32.10 *
17 ! Julie Frohwirth, 14/10/2005 17.44.22 *
18 ! **************************************************************
19
20 USE print_control_mod, ONLY: prt_level, lunout
21 IMPLICIT NONE
22
23 include "cvthermo.h"
24 include "cv3param.h"
25 include "cvflag.h"
26 include "YOMCST2.h"
27 include "YOMCST.h"
28 include "conema3.h"
29
30 ! input:
31 INTEGER, INTENT (IN) :: ncum, nd, nloc
32 INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
33 REAL, DIMENSION (nloc), INTENT (IN) :: pbase, plcl
34 REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
35 REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
36 REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, buoy
37 REAL, DIMENSION (nloc, nd), INTENT (IN) :: supmax
38 LOGICAL, INTENT (IN) :: ok_inhib ! enable convection inhibition by dryness
39 REAL, DIMENSION (nloc), INTENT (IN) :: ale, alp
40 REAL, DIMENSION (nloc, nd), INTENT (IN) :: omega
41
42 ! input/output:
43 INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
44 REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig, w0
45 REAL, DIMENSION (nloc), INTENT (INOUT) :: ptop2
46
47 ! output:
48 REAL, DIMENSION (nloc), INTENT (OUT) :: cape, cin
49 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: m
50 REAL, DIMENSION (nloc), INTENT (OUT) :: plim1, plim2
51 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: asupmax
52 REAL, DIMENSION (nloc), INTENT (OUT) :: supmax0
53 REAL, DIMENSION (nloc), INTENT (OUT) :: asupmaxmin
54 REAL, DIMENSION (nloc), INTENT (OUT) :: cbmflast, plfc
55 REAL, DIMENSION (nloc), INTENT (OUT) :: wbeff
56
57 ! local variables:
58 INTEGER :: il, i, j, k, icbmax
59 INTEGER, DIMENSION (nloc) :: i0, klfc
60 REAL :: deltap, fac, w, amu
61 REAL, DIMENSION (nloc, nd) :: rhodp ! Factor such that m=rhodp*sig*w
62 REAL :: dz
63 REAL :: pbmxup
64 REAL, DIMENSION (nloc, nd) :: dtmin, sigold
65 REAL, DIMENSION (nloc, nd) :: coefmix
66 REAL, DIMENSION (nloc) :: dtminmax
67 REAL, DIMENSION (nloc) :: pzero, ptop2old
68 REAL, DIMENSION (nloc) :: cina, cinb
69 INTEGER, DIMENSION (nloc) :: ibeg
70 INTEGER, DIMENSION (nloc) :: nsupmax
71 REAL :: supcrit
72 REAL, DIMENSION (nloc, nd) :: temp
73 REAL, DIMENSION (nloc) :: p1, pmin
74 REAL, DIMENSION (nloc) :: asupmax0
75 LOGICAL, DIMENSION (nloc) :: ok
76 REAL, DIMENSION (nloc, nd) :: siglim, wlim, mlim
77 REAL, DIMENSION (nloc) :: wb2
78 REAL, DIMENSION (nloc) :: cbmf0 ! initial cloud base mass flux
79 REAL, DIMENSION (nloc) :: cbmflim ! cbmf given by Cape closure
80 REAL, DIMENSION (nloc) :: cbmfalp ! cbmf given by Alp closure
81 REAL, DIMENSION (nloc) :: cbmfalpb ! bounded cbmf given by Alp closure
82 REAL, DIMENSION (nloc) :: cbmfmax ! upper bound on cbmf
83 REAL, DIMENSION (nloc) :: coef
84 REAL, DIMENSION (nloc) :: xp, xq, xr, discr, b3, b4
85 REAL, DIMENSION (nloc) :: theta, bb
86 REAL :: term1, term2, term3
87 REAL, DIMENSION (nloc) :: alp2 ! Alp with offset
88
89 !CR: variables for new erosion of adiabiatic ascent
90 REAL, DIMENSION (nloc, nd) :: mad, me, betalim, beta_coef
91 REAL, DIMENSION (nloc, nd) :: med, md
92 !jyg<
93 ! coef_peel is now in the common cv3_param
94 !! REAL :: coef_peel
95 !! PARAMETER (coef_peel=0.25)
96 !>jyg
97
98 REAL :: sigmax
99 PARAMETER (sigmax=0.1)
100 !! PARAMETER (sigmax=10.)
101
102 CHARACTER (LEN=20) :: modname = 'cv3p2_closure'
103 CHARACTER (LEN=80) :: abort_message
104
105 INTEGER,SAVE :: igout=1
106 !$OMP THREADPRIVATE(igout)
107
108 IF (prt_level>=20) print *,' -> cv3p2_closure, Ale ',ale(igout)
109
110
111 ! -------------------------------------------------------
112 ! -- Initialization
113 ! -------------------------------------------------------
114
115
116 DO il = 1, ncum
117 alp2(il) = max(alp(il), 1.E-5)
118 ! IM
119 alp2(il) = max(alp(il), 1.E-12)
120 END DO
121
122 pbmxup = 50. ! PBMXUP+PBCRIT = cloud depth above which mixed updraughts
123 ! exist (if any)
124
125 IF (prt_level>=20) PRINT *, 'cv3p2_closure nloc ncum nd icb inb nl', nloc, &
126 ncum, nd, icb(nloc), inb(nloc), nl
127 DO k = 1, nl
128 DO il = 1, ncum
129 rhodp(il,k) = 0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
130 END DO
131 END DO
132
133 !CR+jyg: initializations (up to nd) for erosion of adiabatic ascent and of m and wlim
134 DO k = 1,nd
135 DO il = 1, ncum
136 mad(il,k)=0.
137 me(il,k)=0.
138 betalim(il,k)=1.
139 wlim(il,k)=0.
140 m(il, k) = 0.0
141 ENDDO
142 ENDDO
143
144 ! -------------------------------------------------------
145 ! -- Reset sig(i) and w0(i) for i>inb and i<icb
146 ! -------------------------------------------------------
147
148 ! update sig and w0 above LNB:
149
150 DO k = 1, nl - 1
151 DO il = 1, ncum
152 IF ((inb(il)<(nl-1)) .AND. (k>=(inb(il)+1))) THEN
153 sig(il, k) = beta*sig(il, k) + 2.*alpha*buoy(il, inb(il))*abs(buoy(il,inb(il)))
154 sig(il, k) = amax1(sig(il,k), 0.0)
155 w0(il, k) = beta*w0(il, k)
156 END IF
157 END DO
158 END DO
159
160 ! if(prt.level.GE.20) print*,'cv3p2_closure apres 100'
161 ! compute icbmax:
162
163 icbmax = 2
164 DO il = 1, ncum
165 icbmax = max(icbmax, icb(il))
166 END DO
167 ! if(prt.level.GE.20) print*,'cv3p2_closure apres 200'
168
169 ! update sig and w0 below cloud base:
170
171 DO k = 1, icbmax
172 DO il = 1, ncum
173 IF (k<=icb(il)) THEN
174 sig(il, k) = beta*sig(il, k) - 2.*alpha*buoy(il, icb(il))*buoy(il,icb(il))
175 sig(il, k) = amax1(sig(il,k), 0.0)
176 w0(il, k) = beta*w0(il, k)
177 END IF
178 END DO
179 END DO
180 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 300'
181
182 ! -------------------------------------------------------------
183 ! -- Reset fractional areas of updrafts and w0 at initial time
184 ! -- and after 10 time steps of no convection
185 ! -------------------------------------------------------------
186
187 !jyg<
188 IF (ok_convstop) THEN
189 DO k = 1, nl - 1
190 DO il = 1, ncum
191 IF (sig(il,nd)<1.5 .OR. sig(il,nd)>noconv_stop) THEN
192 sig(il, k) = 0.0
193 w0(il, k) = 0.0
194 END IF
195 END DO
196 END DO
197 ELSE
198 DO k = 1, nl - 1
199 DO il = 1, ncum
200 IF (sig(il,nd)<1.5 .OR. sig(il,nd)>12.0) THEN
201 sig(il, k) = 0.0
202 w0(il, k) = 0.0
203 END IF
204 END DO
205 END DO
206 ENDIF ! (ok_convstop)
207 !>jyg
208 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 400'
209
210 ! -------------------------------------------------------
211 ! -- Compute initial cloud base mass flux (Cbmf0)
212 ! -------------------------------------------------------
213 DO il = 1, ncum
214 cbmf0(il) = 0.0
215 END DO
216
217 DO k = 1, nl
218 DO il = 1, ncum
219 IF (k>=icb(il) .AND. k<=inb(il) &
220 .AND. icb(il)+1<=inb(il)) THEN
221 cbmf0(il) = cbmf0(il) + sig(il, k)*w0(il,k)*rhodp(il,k)
222 END IF
223 END DO
224 END DO
225
226 ! -------------------------------------------------------------
227 ! jyg1
228 ! -- Calculate adiabatic ascent top pressure (ptop)
229 ! -------------------------------------------------------------
230
231
232 ! c 1. Start at first level where precipitations form
233 DO il = 1, ncum
234 pzero(il) = plcl(il) - pbcrit
235 END DO
236
237 ! c 2. Add offset
238 DO il = 1, ncum
239 pzero(il) = pzero(il) - pbmxup
240 END DO
241 DO il = 1, ncum
242 ptop2old(il) = ptop2(il)
243 END DO
244
245 DO il = 1, ncum
246 ! CR:c est quoi ce 300??
247 p1(il) = pzero(il) - 300.
248 END DO
249
250 ! compute asupmax=abs(supmax) up to lnm+1
251
252 DO il = 1, ncum
253 ok(il) = .TRUE.
254 nsupmax(il) = inb(il)
255 END DO
256
257 DO i = 1, nl
258 DO il = 1, ncum
259 IF (i>icb(il) .AND. i<=inb(il)) THEN
260 IF (p(il,i)<=pzero(il) .AND. supmax(il,i)<0 .AND. ok(il)) THEN
261 nsupmax(il) = i
262 ok(il) = .FALSE.
263 END IF ! end IF (P(i) ... )
264 END IF ! end IF (icb+1 le i le inb)
265 END DO
266 END DO
267
268 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 2.'
269 DO i = 1, nl
270 DO il = 1, ncum
271 asupmax(il, i) = abs(supmax(il,i))
272 END DO
273 END DO
274
275
276 DO il = 1, ncum
277 asupmaxmin(il) = 10.
278 pmin(il) = 100.
279 ! IM ??
280 asupmax0(il) = 0.
281 END DO
282
283 ! c 3. Compute in which level is Pzero
284
285 ! IM bug i0 = 18
286 DO il = 1, ncum
287 i0(il) = nl
288 END DO
289
290 DO i = 1, nl
291 DO il = 1, ncum
292 IF (i>icb(il) .AND. i<=inb(il)) THEN
293 IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
294 IF (pzero(il)>p(il,i) .AND. pzero(il)<p(il,i-1)) THEN
295 i0(il) = i
296 END IF
297 END IF
298 END IF
299 END DO
300 END DO
301 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 3.'
302
303 ! c 4. Compute asupmax at Pzero
304
305 DO i = 1, nl
306 DO il = 1, ncum
307 IF (i>icb(il) .AND. i<=inb(il)) THEN
308 IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
309 asupmax0(il) = ((pzero(il)-p(il,i0(il)-1))*asupmax(il,i0(il))- &
310 (pzero(il)-p(il,i0(il)))*asupmax(il,i0(il)-1))/(p(il,i0(il))-p(il,i0(il)-1))
311 END IF
312 END IF
313 END DO
314 END DO
315
316
317 DO i = 1, nl
318 DO il = 1, ncum
319 IF (p(il,i)==pzero(il)) THEN
320 asupmax(i, il) = asupmax0(il)
321 END IF
322 END DO
323 END DO
324 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 4.'
325
326 ! c 5. Compute asupmaxmin, minimum of asupmax
327
328 DO i = 1, nl
329 DO il = 1, ncum
330 IF (i>icb(il) .AND. i<=inb(il)) THEN
331 IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
332 IF (asupmax(il,i)<asupmaxmin(il)) THEN
333 asupmaxmin(il) = asupmax(il, i)
334 pmin(il) = p(il, i)
335 END IF
336 END IF
337 END IF
338 END DO
339 END DO
340
341 DO il = 1, ncum
342 ! IM
343 IF (prt_level>=20) THEN
344 PRINT *, 'cv3p2_closure il asupmax0 asupmaxmin', il, asupmax0(il), &
345 asupmaxmin(il), pzero(il), pmin(il)
346 END IF
347 IF (asupmax0(il)<asupmaxmin(il)) THEN
348 asupmaxmin(il) = asupmax0(il)
349 pmin(il) = pzero(il)
350 END IF
351 END DO
352 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 5.'
353
354
355 ! Compute Supmax at Pzero
356
357 DO i = 1, nl
358 DO il = 1, ncum
359 IF (i>icb(il) .AND. i<=inb(il)) THEN
360 IF (p(il,i)<=pzero(il)) THEN
361 supmax0(il) = ((p(il,i)-pzero(il))*asupmax(il,i-1)- &
362 (p(il,i-1)-pzero(il))*asupmax(il,i))/(p(il,i)-p(il,i-1))
363 GO TO 425
364 END IF ! end IF (P(i) ... )
365 END IF ! end IF (icb+1 le i le inb)
366 END DO
367 END DO
368
369 425 CONTINUE
370 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 425.'
371
372 ! c 6. Calculate ptop2
373
374 DO il = 1, ncum
375 IF (asupmaxmin(il)<supcrit1) THEN
376 ptop2(il) = pmin(il)
377 END IF
378
379 IF (asupmaxmin(il)>supcrit1 .AND. asupmaxmin(il)<supcrit2) THEN
380 ptop2(il) = ptop2old(il)
381 END IF
382
383 IF (asupmaxmin(il)>supcrit2) THEN
384 ptop2(il) = ph(il, inb(il))
385 END IF
386 END DO
387
388 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 6.'
389
390 ! c 7. Compute multiplying factor for adiabatic updraught mass flux
391
392
393 IF (ok_inhib) THEN
394
395 DO i = 1, nl
396 DO il = 1, ncum
397 IF (i<=nl) THEN
398 coefmix(il, i) = (min(ptop2(il),ph(il,i))-ph(il,i))/(ph(il,i+1)-ph(il,i))
399 coefmix(il, i) = min(coefmix(il,i), 1.)
400 END IF
401 END DO
402 END DO
403
404
405 ELSE ! when inhibition is not taken into account, coefmix=1
406
407
408
409 DO i = 1, nl
410 DO il = 1, ncum
411 IF (i<=nl) THEN
412 coefmix(il, i) = 1.
413 END IF
414 END DO
415 END DO
416
417 END IF ! ok_inhib
418 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 7.'
419 ! -------------------------------------------------------------------
420 ! -------------------------------------------------------------------
421
422
423 ! jyg2
424
425 ! ==========================================================================
426
427
428 ! -------------------------------------------------------------
429 ! -- Calculate convective inhibition (CIN)
430 ! -------------------------------------------------------------
431
432 ! do i=1,nloc
433 ! print*,'avant cine p',pbase(i),plcl(i)
434 ! enddo
435 ! do j=1,nd
436 ! do i=1,nloc
437 ! print*,'avant cine t',tv(i),tvp(i)
438 ! enddo
439 ! enddo
440 CALL cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, &
441 cinb, plfc)
442
443 DO il = 1, ncum
444 cin(il) = cina(il) + cinb(il)
445 END DO
446 IF (prt_level>=20) PRINT *, 'cv3p2_closure after cv3_cine: cina, cinb, cin ', &
447 cina(igout), cinb(igout), cin(igout)
448 ! -------------------------------------------------------------
449 ! --Update buoyancies to account for Ale
450 ! -------------------------------------------------------------
451
452 CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, &
453 tvp, buoy)
454 IF (prt_level>=20) PRINT *, 'cv3p2_closure after cv3_buoy'
455
456 ! -------------------------------------------------------------
457 ! -- Calculate convective available potential energy (cape),
458 ! -- vertical velocity (w), fractional area covered by
459 ! -- undilute updraft (sig), and updraft mass flux (m)
460 ! -------------------------------------------------------------
461
462 DO il = 1, ncum
463 cape(il) = 0.0
464 dtminmax(il) = -100.
465 END DO
466
467 ! compute dtmin (minimum buoyancy between ICB and given level k):
468
469 DO k = 1, nl
470 DO il = 1, ncum
471 dtmin(il, k) = 100.0
472 END DO
473 END DO
474
475 DO k = 1, nl
476 DO j = minorig, nl
477 DO il = 1, ncum
478 IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (j>=icb(il)) &
479 .AND. (j<=(k-1))) THEN
480 dtmin(il, k) = amin1(dtmin(il,k), buoy(il,j))
481 END IF
482 END DO
483 END DO
484 END DO
485 !jyg<
486 ! Store maximum of dtmin
487 ! C est pas terrible d avoir ce test sur Ale+Cin encore une fois ici.
488 ! A REVOIR !
489 DO k = 1, nl
490 DO il = 1, ncum
491 IF (k>=(icb(il)+1) .AND. k<=inb(il) .AND. ale(il)+cin(il)>0.) THEN
492 dtminmax(il) = max(dtmin(il,k), dtminmax(il))
493 ENDIF
494 END DO
495 END DO
496 !
497 ! prevent convection when ale+cin <= 0
498 DO k = 1, nl
499 DO il = 1, ncum
500 IF (k>=(icb(il)+1) .AND. k<=inb(il)) THEN
501 dtmin(il,k) = min(dtmin(il,k), dtminmax(il))
502 ENDIF
503 END DO
504 END DO
505 !>jyg
506 !
507 IF (prt_level >= 20) THEN
508 print *,'cv3p2_closure: dtmin ', (k, dtmin(igout,k), k=1,nl)
509 print *,'cv3p2_closure: dtminmax ', dtminmax(igout)
510 ENDIF
511 !
512 ! the interval on which cape is computed starts at pbase :
513
514 DO k = 1, nl
515 DO il = 1, ncum
516
517 IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
518
519 IF (iflag_mix_adiab.eq.1) THEN
520 !CR:computation of cape from LCL: keep flag or to modify in all cases?
521 deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k))
522 ELSE
523 deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
524 ENDIF
525 cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1)
526 cape(il) = amax1(0.0, cape(il))
527 sigold(il, k) = sig(il, k)
528
529
530 ! jyg Coefficient coefmix limits convection to levels where a
531 ! sufficient
532 ! fraction of mixed draughts are ascending.
533 siglim(il, k) = coefmix(il, k)*alpha1*dtmin(il, k)*abs(dtmin(il,k))
534 siglim(il, k) = amax1(siglim(il,k), 0.0)
535 siglim(il, k) = amin1(siglim(il,k), 0.01)
536 ! c fac=AMIN1(((dtcrit-dtmin(il,k))/dtcrit),1.0)
537 fac = 1.
538 wlim(il, k) = fac*sqrt(cape(il))
539 amu = siglim(il, k)*wlim(il, k)
540 !! rhodp(il,k) = 0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k) !cor jyg : computed earlier
541 mlim(il, k) = amu*rhodp(il,k)
542 ! print*, 'siglim ', k,siglim(1,k)
543 END IF
544
545 END DO
546 END DO
547 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 600'
548
549 DO il = 1, ncum
550 ! IM beg
551 IF (prt_level>=20) THEN
552 PRINT *, 'cv3p2_closure il icb mlim ph ph+1 ph+2', il, icb(il), &
553 mlim(il, icb(il)+1), ph(il, icb(il)), ph(il, icb(il)+1), &
554 ph(il, icb(il)+2)
555 END IF
556
557 IF (icb(il)+1<=inb(il)) THEN
558 ! IM end
559 mlim(il, icb(il)) = 0.5*mlim(il,icb(il)+1)*(ph(il,icb(il))-ph(il,icb(il)+1))/ &
560 (ph(il,icb(il)+1)-ph(il,icb(il)+2))
561 ! IM beg
562 END IF !(icb(il.le.inb(il))) then
563 ! IM end
564 END DO
565 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 700'
566
567 !
568 ! ------------------------------------------------------------------------
569 ! c Compute Cloud base mass flux given by Cape closure (cbmflim = cbmf of
570 ! c elementary systems), cbmf given by Alp closure (cbmfalp), cbmf given by Alp
571 ! c closure with an upper bound imposed (cbmfalpb) and cbmf resulting from
572 ! c time integration (cbmflast).
573 ! ------------------------------------------------------------------------
574
575 DO il = 1, ncum
576 cbmflim(il) = 0.
577 cbmfalp(il) = 0.
578 cbmfalpb(il) = 0.
579 cbmflast(il) = 0.
580 END DO
581
582 ! c 1. Compute cloud base mass flux of elementary system (Cbmflim)
583
584 DO k = 1, nl
585 DO il = 1, ncum
586 ! old IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
587 ! IM IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
588 IF (k>=icb(il) .AND. k<=inb(il) & !cor jyg
589 .AND. icb(il)+1<=inb(il)) THEN !cor jyg
590 cbmflim(il) = cbmflim(il) + mlim(il, k)
591 END IF
592 END DO
593 END DO
594 IF (prt_level>=20) PRINT *, 'cv3p2_closure after cbmflim: cbmflim ', cbmflim(igout)
595
596 ! 1.5 Compute cloud base mass flux given by Alp closure (Cbmfalp), maximum
597 ! allowed mass flux (Cbmfmax) and bounded mass flux (Cbmfalpb)
598 ! Cbmfalpb is set to zero if Cbmflim (the mass flux of elementary cloud)
599 ! is exceedingly small.
600
601 DO il = 1, ncum
602 wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.))
603 END DO
604
605 DO il = 1, ncum
606 IF (plfc(il)<100.) THEN
607 ! This is an irealistic value for plfc => no calculation of wbeff
608 wbeff(il) = 100.1
609 ELSE
610 ! Calculate wbeff
611 IF (NINT(flag_wb)==0) THEN
612 wbeff(il) = wbmax
613 ELSE IF (NINT(flag_wb)==1) THEN
614 wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
615 ELSE IF (NINT(flag_wb)==2) THEN
616 wbeff(il) = wbmax*(0.01*(ph(il,1)-plfc(il)))**2
617 END IF
618 END IF
619 END DO
620
621 !CR:Compute k at plfc
622 DO il=1,ncum
623 klfc(il)=nl
624 ENDDO
625 DO k=1,nl
626 DO il=1,ncum
627 if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
628 klfc(il)=k
629 endif
630 ENDDO
631 ENDDO
632 !RC
633
634 DO il = 1, ncum
635 ! jyg Modification du coef de wb*wb pour conformite avec papier Wake
636 ! c cbmfalp(il) = alp2(il)/(0.5*wb*wb-Cin(il))
637 cbmfalp(il) = alp2(il)/(2.*wbeff(il)*wbeff(il)-cin(il))
638 !CR: Add large-scale component to the mass-flux
639 !encore connu sous le nom "Experience du tube de dentifrice"
640 if ((coef_clos_ls.gt.0.).and.(plfc(il).gt.0.)) then
641 cbmfalp(il) = cbmfalp(il) - coef_clos_ls*min(0.,1./RG*omega(il,klfc(il)))
642 endif
643 !RC
644 IF (cbmfalp(il)==0 .AND. alp2(il)/=0.) THEN
645 WRITE (lunout, *) 'cv3p2_closure cbmfalp=0 and alp NE 0 il alp2 alp cin ' , &
646 il, alp2(il), alp(il), cin(il)
647 abort_message = ''
648 CALL abort_physic(modname, abort_message, 1)
649 END IF
650 cbmfmax(il) = sigmax*wb2(il)*100.*p(il, icb(il))/(rrd*tv(il,icb(il)))
651 END DO
652
653 !jyg<
654 IF (OK_intermittent) THEN
655 DO il = 1, ncum
656 IF (cbmflim(il)>1.E-6) THEN
657 cbmfalpb(il) = min(cbmfalp(il), (cbmfmax(il)-beta*cbmf0(il))/(1.-beta))
658 ! print*,'cbmfalpb',cbmfalpb(il),cbmfmax(il)
659 END IF
660 END DO
661 ELSE
662 !>jyg
663 DO il = 1, ncum
664 IF (cbmflim(il)>1.E-6) THEN
665 ! ATTENTION TEST CR
666 ! if (cbmfmax(il).lt.1.e-12) then
667 cbmfalpb(il) = min(cbmfalp(il), cbmfmax(il))
668 ! else
669 ! cbmfalpb(il) = cbmfalp(il)
670 ! endif
671 ! print*,'cbmfalpb',cbmfalp(il),cbmfmax(il)
672 END IF
673 END DO
674 ENDIF !(OK_intermittent)
675 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres cbmfalpb: cbmfalpb ',cbmfalpb(igout)
676
677 ! c 2. Compute coefficient and apply correction
678
679 DO il = 1, ncum
680 coef(il) = (cbmfalpb(il)+1.E-10)/(cbmflim(il)+1.E-10)
681 END DO
682 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres coef_plantePLUS'
683
684 DO k = 1, nl
685 DO il = 1, ncum
686 IF (k>=icb(il)+1 .AND. k<=inb(il)) THEN
687 amu = beta*sig(il, k)*w0(il, k) + (1.-beta)*coef(il)*siglim(il, k)*wlim(il, k)
688 w0(il, k) = wlim(il, k)
689 w0(il, k) = max(w0(il,k), 1.E-10)
690 sig(il, k) = amu/w0(il, k)
691 sig(il, k) = min(sig(il,k), 1.)
692 ! c amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k)
693 !jyg m(il, k) = amu*0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
694 m(il, k) = amu*rhodp(il,k)
695 END IF
696 END DO
697 END DO
698 ! jyg2
699 DO il = 1, ncum
700 w0(il, icb(il)) = 0.5*w0(il, icb(il)+1)
701 m(il, icb(il)) = 0.5*m(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb(il)+1))/ &
702 (ph(il,icb(il)+1)-ph(il,icb(il)+2))
703 sig(il, icb(il)) = sig(il, icb(il)+1)
704 sig(il, icb(il)-1) = sig(il, icb(il))
705 END DO
706 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres w0_sig_M: w0, sig ', &
707 (k,w0(igout,k),sig(igout,k), k=icb(igout),inb(igout))
708
709 !CR: new erosion of adiabatic ascent: modification of m
710 !computation of the sum of ascending fluxes
711 IF (iflag_mix_adiab.eq.1) THEN
712
713 !Verification sum(me)=sum(m)
714 DO k = 1,nd
715 DO il = 1, ncum
716 md(il,k)=0.
717 med(il,k)=0.
718 ENDDO
719 ENDDO
720
721 DO k = nl,1,-1
722 DO il = 1, ncum
723 md(il,k)=md(il,k+1)+m(il,k+1)
724 ENDDO
725 ENDDO
726
727 DO k = nl,1,-1
728 DO il = 1, ncum
729 IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
730 mad(il,k)=mad(il,k+1)+m(il,k+1)
731 ENDIF
732 ! print*,"mad",il,k,mad(il,k)
733 ENDDO
734 ENDDO
735
736 !CR: erosion of each adiabatic ascent during its ascent
737
738 !Computation of erosion coefficient beta_coef
739 DO k = 1, nl
740 DO il = 1, ncum
741 IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN
742 ! print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
743 beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
744 ELSE
745 beta_coef(il,k)=0.
746 ENDIF
747 ENDDO
748 ENDDO
749
750 ! print*,"apres beta_coef"
751
752 DO k = 1, nl
753 DO il = 1, ncum
754
755 IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
756
757 ! print*,"dz",il,k,tv(il, k-1)
758 dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG)
759 betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz)
760 ! betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz)
761 ! print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
762 dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG)
763 ! me(il,k)=betalim(il,k)*(m(il,k)+RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz*mad(il,k))
764 me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k))
765 ! print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz
766
767 END IF
768
769 !Modification of m
770 m(il,k)=me(il,k)
771 END DO
772 END DO
773
774 ! DO il = 1, ncum
775 ! dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG)
776 ! m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) &
777 ! /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il))
778 ! print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
779 ! ENDDO
780
781 !Verification sum(me)=sum(m)
782 DO k = nl,1,-1
783 DO il = 1, ncum
784 med(il,k)=med(il,k+1)+m(il,k+1)
785 ! print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
786 ENDDO
787 ENDDO
788
789
790 ENDIF !(iflag_mix_adiab)
791 !RC
792
793 ! c 3. Compute final cloud base mass flux;
794 ! c set iflag to 3 if cloud base mass flux is exceedingly small and is
795 ! c decreasing (i.e. if the final mass flux (cbmflast) is greater than
796 ! c the target mass flux (cbmfalpb)).
797 ! c If(ok_convstop): set iflag to 4 if no positive buoyancy has been met
798
799 !jyg DO il = 1, ncum
800 !jyg cbmflast(il) = 0.
801 !jyg END DO
802
803 DO k = 1, nl
804 DO il = 1, ncum
805 IF (k>=icb(il) .AND. k<=inb(il)) THEN
806 !IMpropo?? IF ((k.ge.(icb(il)+1)).and.(k.le.inb(il))) THEN
807 cbmflast(il) = cbmflast(il) + m(il, k)
808 END IF
809 END DO
810 END DO
811 IF (prt_level>=20) PRINT *, 'cv3p2_closure apres cbmflast: cbmflast ',cbmflast(igout)
812
813 DO il = 1, ncum
814 IF (cbmflast(il)<1.E-6 .AND. cbmflast(il)>=cbmfalpb(il)) THEN
815 iflag(il) = 3
816 END IF
817 END DO
818
819 !jyg<
820 IF (ok_convstop) THEN
821 DO il = 1, ncum
822 IF (dtminmax(il) .LE. 0.) THEN
823 iflag(il) = 4
824 END IF
825 END DO
826 ELSE
827 !>jyg
828 DO k = 1, nl
829 DO il = 1, ncum
830 IF (iflag(il)>=3) THEN
831 m(il, k) = 0.
832 sig(il, k) = 0.
833 w0(il, k) = 0.
834 END IF
835 END DO
836 END DO
837 ENDIF ! (ok_convstop)
838 !
839 IF (prt_level >= 10) THEN
840 print *,'cv3p2_closure: iflag ',iflag(igout)
841 ENDIF
842 !
843
844 ! c 4. Introduce a correcting factor for coef, in order to obtain an
845 ! effective
846 ! c sigdz larger in the present case (using cv3p2_closure) than in the
847 ! old
848 ! c closure (using cv3_closure).
849 IF (1==0) THEN
850 DO il = 1, ncum
851 ! c coef(il) = 2.*coef(il)
852 coef(il) = 5.*coef(il)
853 END DO
854 ! version CVS du ..2008
855 ELSE
856 IF (iflag_cvl_sigd==0) THEN
857 ! test pour verifier qu on fait la meme chose qu avant: sid constant
858 coef(1:ncum) = 1.
859 ELSE
860 coef(1:ncum) = min(2.*coef(1:ncum), 5.)
861 coef(1:ncum) = max(2.*coef(1:ncum), 0.2)
862 END IF
863 END IF
864
865 IF (prt_level>=20) PRINT *, 'cv3p2_closure FIN'
866 RETURN
867 END SUBROUTINE cv3p2_closure
868
869
870