GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/cv3p2_closure.F90 Lines: 0 327 0.0 %
Date: 2023-06-30 12:51:15 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