source: LMDZ6/branches/Amaury_dev/libf/phylmd/cv3p2_closure.F90 @ 5159

Last change on this file since 5159 was 5158, checked in by abarral, 7 weeks ago

Add missing klon on strataer_emiss_mod.F90
Correct various missing explicit declarations
Replace tabs by spaces (tabs are not part of the fortran charset)
Continue cleaning modules
Removed unused arguments and variables
File size: 26.3 KB
Line 
1
2
3SUBROUTINE 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 lmdz_print_control, ONLY: prt_level, lunout
21  USE lmdz_abort_physic, ONLY: abort_physic
22  USE lmdz_conema3
23  USE lmdz_cvflag
24  USE lmdz_cvthermo
25  USE lmdz_cv3param
26  USE lmdz_yomcst
27  USE lmdz_yomcst2
28
29  IMPLICIT NONE
30
31  ! input:
32  INTEGER, INTENT (IN)                               :: ncum, nd, nloc
33  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb
34  REAL, DIMENSION (nloc), INTENT (IN)                :: pbase, plcl
35  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: p
36  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
37  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: tv, tvp, buoy
38  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: supmax
39  LOGICAL, INTENT (IN)                               :: ok_inhib ! enable convection inhibition by dryness
40  REAL, DIMENSION (nloc), INTENT (IN)                :: ale, alp
41  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: omega
42
43  ! input/output:
44  INTEGER, DIMENSION (nloc), INTENT (INOUT)          :: iflag
45  REAL, DIMENSION (nloc, nd), INTENT (INOUT)         :: sig, w0
46  REAL, DIMENSION (nloc), INTENT (INOUT)             :: ptop2
47
48  ! output:
49  REAL, DIMENSION (nloc), INTENT (OUT)               :: cape, cin
50  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: m
51  REAL, DIMENSION (nloc), INTENT (OUT)               :: plim1, plim2
52  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: asupmax
53  REAL, DIMENSION (nloc), INTENT (OUT)               :: supmax0
54  REAL, DIMENSION (nloc), INTENT (OUT)               :: asupmaxmin
55  REAL, DIMENSION (nloc), INTENT (OUT)               :: cbmflast, plfc
56  REAL, DIMENSION (nloc), INTENT (OUT)               :: wbeff
57
58  ! local variables:
59  INTEGER                                            :: il, i, j, k, icbmax
60  INTEGER, DIMENSION (nloc)                          :: i0, klfc
61  REAL                                               :: deltap, fac, w, amu
62  REAL, DIMENSION (nloc, nd)                         :: rhodp               ! Factor such that m=rhodp*sig*w
63  REAL                                               :: dz
64  REAL                                               :: pbmxup
65  REAL, DIMENSION (nloc, nd)                         :: dtmin, sigold
66  REAL, DIMENSION (nloc, nd)                         :: coefmix
67  REAL, DIMENSION (nloc)                             :: dtminmax
68  REAL, DIMENSION (nloc)                             :: pzero, ptop2old
69  REAL, DIMENSION (nloc)                             :: cina, cinb
70  INTEGER, DIMENSION (nloc)                          :: ibeg
71  INTEGER, DIMENSION (nloc)                          :: nsupmax
72  REAL                                               :: supcrit
73  REAL, DIMENSION (nloc, nd)                         :: temp
74  REAL, DIMENSION (nloc)                             :: p1, pmin
75  REAL, DIMENSION (nloc)                             :: asupmax0
76  LOGICAL, DIMENSION (nloc)                          :: ok
77  REAL, DIMENSION (nloc, nd)                         :: siglim, wlim, mlim
78  REAL, DIMENSION (nloc)                             :: wb2
79  REAL, DIMENSION (nloc)                             :: cbmf0        ! initial cloud base mass flux
80  REAL, DIMENSION (nloc)                             :: cbmflim      ! cbmf given by Cape closure
81  REAL, DIMENSION (nloc)                             :: cbmfalp      ! cbmf given by Alp closure
82  REAL, DIMENSION (nloc)                             :: cbmfalpb     ! bounded cbmf given by Alp closure
83  REAL, DIMENSION (nloc)                             :: cbmfmax      ! upper bound on cbmf
84  REAL, DIMENSION (nloc)                             :: coef
85  REAL, DIMENSION (nloc)                             :: xp, xq, xr, discr, b3, b4
86  REAL, DIMENSION (nloc)                             :: theta, bb
87  REAL                                               :: term1, term2, term3
88  REAL, DIMENSION (nloc)                             :: alp2                  ! Alp with offset
89
90!CR: variables for new erosion of adiabiatic ascent
91  REAL, DIMENSION (nloc, nd)                         :: mad, me, betalim, beta_coef
92  REAL, DIMENSION (nloc, nd)                         :: med, md
93!jyg<
94! coef_peel is now in the common cv3_param
95!!  REAL                                               :: coef_peel
96!!  PARAMETER (coef_peel=0.25)
97!>jyg
98
99  REAL                                               :: sigmax
100  PARAMETER (sigmax=0.1)
101!!  PARAMETER (sigmax=10.)
102
103  CHARACTER (LEN=20)                                 :: modname = 'cv3p2_closure'
104  CHARACTER (LEN=80)                                 :: abort_message
105
106  INTEGER,SAVE                                       :: igout=1
107!$OMP THREADPRIVATE(igout)
108
109 IF (prt_level>=20) print *,' -> cv3p2_closure, Ale ',ale(igout)
110
111
112  ! -------------------------------------------------------
113  ! -- Initialization
114  ! -------------------------------------------------------
115
116
117  DO il = 1, ncum
118    alp2(il) = max(alp(il), 1.E-5)
119    ! IM
120    alp2(il) = max(alp(il), 1.E-12)
121  END DO
122
123  pbmxup = 50. ! PBMXUP+PBCRIT = cloud depth above which mixed updraughts
124  ! exist (if any)
125
126  IF (prt_level>=20) PRINT *, 'cv3p2_closure nloc ncum nd icb inb nl', nloc, &
127    ncum, nd, icb(nloc), inb(nloc), nl
128  DO k = 1, nl
129    DO il = 1, ncum
130      rhodp(il,k) = 0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
131    END DO
132  END DO
133
134!CR+jyg: initializations (up to nd) for erosion of adiabatic ascent and of m and wlim
135  DO k = 1,nd
136    DO il = 1, ncum
137        mad(il,k)=0.
138        me(il,k)=0.
139        betalim(il,k)=1.
140        wlim(il,k)=0.
141        m(il, k) = 0.0
142    ENDDO
143  ENDDO
144
145  ! -------------------------------------------------------
146  ! -- Reset sig(i) and w0(i) for i>inb and i<icb
147  ! -------------------------------------------------------
148
149  ! update sig and w0 above LNB:
150
151  DO k = 1, nl - 1
152    DO il = 1, ncum
153      IF ((inb(il)<(nl-1)) .AND. (k>=(inb(il)+1))) THEN
154        sig(il, k) = beta*sig(il, k) + 2.*alpha*buoy(il, inb(il))*abs(buoy(il,inb(il)))
155        sig(il, k) = amax1(sig(il,k), 0.0)
156        w0(il, k) = beta*w0(il, k)
157      END IF
158    END DO
159  END DO
160
161  ! IF(prt.level.GE.20) PRINT*,'cv3p2_closure apres 100'
162  ! compute icbmax:
163
164  icbmax = 2
165  DO il = 1, ncum
166    icbmax = max(icbmax, icb(il))
167  END DO
168  ! IF(prt.level.GE.20) PRINT*,'cv3p2_closure apres 200'
169
170  ! update sig and w0 below cloud base:
171
172  DO k = 1, icbmax
173    DO il = 1, ncum
174      IF (k<=icb(il)) THEN
175        sig(il, k) = beta*sig(il, k) - 2.*alpha*buoy(il, icb(il))*buoy(il,icb(il))
176        sig(il, k) = amax1(sig(il,k), 0.0)
177        w0(il, k) = beta*w0(il, k)
178      END IF
179    END DO
180  END DO
181  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 300'
182
183  ! -------------------------------------------------------------
184  ! -- Reset fractional areas of updrafts and w0 at initial time
185  ! -- and after 10 time steps of no convection
186  ! -------------------------------------------------------------
187
188!jyg<
189  IF (ok_convstop) THEN
190    DO k = 1, nl - 1
191      DO il = 1, ncum
192        IF (sig(il,nd)<1.5 .OR. sig(il,nd)>noconv_stop) THEN
193          sig(il, k) = 0.0
194          w0(il, k) = 0.0
195        END IF
196      END DO
197    END DO
198  ELSE
199  DO k = 1, nl - 1
200    DO il = 1, ncum
201      IF (sig(il,nd)<1.5 .OR. sig(il,nd)>12.0) THEN
202        sig(il, k) = 0.0
203        w0(il, k) = 0.0
204      END IF
205    END DO
206  END DO
207  ENDIF  ! (ok_convstop)
208!>jyg
209  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 400'
210
211  ! -------------------------------------------------------
212  ! -- Compute initial cloud base mass flux (Cbmf0)
213  ! -------------------------------------------------------
214  DO il = 1, ncum
215    cbmf0(il) = 0.0
216  END DO
217
218  DO k = 1, nl
219    DO il = 1, ncum
220      IF (k>=icb(il) .AND. k<=inb(il) &
221          .AND. icb(il)+1<=inb(il)) THEN
222        cbmf0(il) = cbmf0(il) + sig(il, k)*w0(il,k)*rhodp(il,k)
223      END IF
224    END DO
225  END DO
226
227  ! -------------------------------------------------------------
228  ! jyg1
229  ! --  Calculate adiabatic ascent top pressure (ptop)
230  ! -------------------------------------------------------------
231
232
233  ! c 1. Start at first level where precipitations form
234  DO il = 1, ncum
235    pzero(il) = plcl(il) - pbcrit
236  END DO
237
238  ! c 2. Add offset
239  DO il = 1, ncum
240    pzero(il) = pzero(il) - pbmxup
241  END DO
242  DO il = 1, ncum
243    ptop2old(il) = ptop2(il)
244  END DO
245
246  DO il = 1, ncum
247    ! CR:c est quoi ce 300??
248    p1(il) = pzero(il) - 300.
249  END DO
250
251  ! compute asupmax=abs(supmax) up to lnm+1
252
253  DO il = 1, ncum
254    ok(il) = .TRUE.
255    nsupmax(il) = inb(il)
256  END DO
257
258  DO i = 1, nl
259    DO il = 1, ncum
260      IF (i>icb(il) .AND. i<=inb(il)) THEN
261        IF (p(il,i)<=pzero(il) .AND. supmax(il,i)<0 .AND. ok(il)) THEN
262          nsupmax(il) = i
263          ok(il) = .FALSE.
264        END IF ! end IF (P(i) ...  )
265      END IF ! end IF (icb+1 le i le inb)
266    END DO
267  END DO
268
269  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 2.'
270  DO i = 1, nl
271    DO il = 1, ncum
272      asupmax(il, i) = abs(supmax(il,i))
273    END DO
274  END DO
275
276
277  DO il = 1, ncum
278    asupmaxmin(il) = 10.
279    pmin(il) = 100.
280    ! IM ??
281    asupmax0(il) = 0.
282  END DO
283
284  ! c 3.  Compute in which level is Pzero
285
286  ! IM bug      i0 = 18
287  DO il = 1, ncum
288    i0(il) = nl
289  END DO
290
291  DO i = 1, nl
292    DO il = 1, ncum
293      IF (i>icb(il) .AND. i<=inb(il)) THEN
294        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
295          IF (pzero(il)>p(il,i) .AND. pzero(il)<p(il,i-1)) THEN
296            i0(il) = i
297          END IF
298        END IF
299      END IF
300    END DO
301  END DO
302  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 3.'
303
304  ! c 4.  Compute asupmax at Pzero
305
306  DO i = 1, nl
307    DO il = 1, ncum
308      IF (i>icb(il) .AND. i<=inb(il)) THEN
309        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
310          asupmax0(il) = ((pzero(il)-p(il,i0(il)-1))*asupmax(il,i0(il))- &
311            (pzero(il)-p(il,i0(il)))*asupmax(il,i0(il)-1))/(p(il,i0(il))-p(il,i0(il)-1))
312        END IF
313      END IF
314    END DO
315  END DO
316
317
318  DO i = 1, nl
319    DO il = 1, ncum
320      IF (p(il,i)==pzero(il)) THEN
321        asupmax(i, il) = asupmax0(il)
322      END IF
323    END DO
324  END DO
325  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 4.'
326
327  ! c 5. Compute asupmaxmin, minimum of asupmax
328
329  DO i = 1, nl
330    DO il = 1, ncum
331      IF (i>icb(il) .AND. i<=inb(il)) THEN
332        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
333          IF (asupmax(il,i)<asupmaxmin(il)) THEN
334            asupmaxmin(il) = asupmax(il, i)
335            pmin(il) = p(il, i)
336          END IF
337        END IF
338      END IF
339    END DO
340  END DO
341
342  DO il = 1, ncum
343    ! IM
344    IF (prt_level>=20) THEN
345      PRINT *, 'cv3p2_closure il asupmax0 asupmaxmin', il, asupmax0(il), &
346        asupmaxmin(il), pzero(il), pmin(il)
347    END IF
348    IF (asupmax0(il)<asupmaxmin(il)) THEN
349      asupmaxmin(il) = asupmax0(il)
350      pmin(il) = pzero(il)
351    END IF
352  END DO
353  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 5.'
354
355
356  ! Compute Supmax at Pzero
357
358  DO i = 1, nl
359    DO il = 1, ncum
360      IF (i>icb(il) .AND. i<=inb(il)) THEN
361        IF (p(il,i)<=pzero(il)) THEN
362          supmax0(il) = ((p(il,i)-pzero(il))*asupmax(il,i-1)- &
363            (p(il,i-1)-pzero(il))*asupmax(il,i))/(p(il,i)-p(il,i-1))
364          GO TO 425
365        END IF ! end IF (P(i) ... )
366      END IF ! end IF (icb+1 le i le inb)
367    END DO
368  END DO
369
370425 CONTINUE
371  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 425.'
372
373  ! c 6. Calculate ptop2
374
375  DO il = 1, ncum
376    IF (asupmaxmin(il)<supcrit1) THEN
377      ptop2(il) = pmin(il)
378    END IF
379
380    IF (asupmaxmin(il)>supcrit1 .AND. asupmaxmin(il)<supcrit2) THEN
381      ptop2(il) = ptop2old(il)
382    END IF
383
384    IF (asupmaxmin(il)>supcrit2) THEN
385      ptop2(il) = ph(il, inb(il))
386    END IF
387  END DO
388
389  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 6.'
390
391  ! c 7. Compute multiplying factor for adiabatic updraught mass flux
392
393
394  IF (ok_inhib) THEN
395
396    DO i = 1, nl
397      DO il = 1, ncum
398        IF (i<=nl) THEN
399          coefmix(il, i) = (min(ptop2(il),ph(il,i))-ph(il,i))/(ph(il,i+1)-ph(il,i))
400          coefmix(il, i) = min(coefmix(il,i), 1.)
401        END IF
402      END DO
403    END DO
404
405
406  ELSE ! when inhibition is not taken into account, coefmix=1
407
408
409
410    DO i = 1, nl
411      DO il = 1, ncum
412        IF (i<=nl) THEN
413          coefmix(il, i) = 1.
414        END IF
415      END DO
416    END DO
417
418  END IF ! ok_inhib
419  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 7.'
420  ! -------------------------------------------------------------------
421  ! -------------------------------------------------------------------
422
423
424  ! jyg2
425
426  ! ==========================================================================
427
428
429  ! -------------------------------------------------------------
430  ! -- Calculate convective inhibition (CIN)
431  ! -------------------------------------------------------------
432
433  ! do i=1,nloc
434  ! PRINT*,'avant cine p',pbase(i),plcl(i)
435  ! enddo
436  ! do j=1,nd
437  ! do i=1,nloc
438  ! PRINT*,'avant cine t',tv(i),tvp(i)
439  ! enddo
440  ! enddo
441  CALL cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, &
442    cinb, plfc)
443
444  DO il = 1, ncum
445    cin(il) = cina(il) + cinb(il)
446  END DO
447  IF (prt_level>=20) PRINT *, 'cv3p2_closure after cv3_cine: cina, cinb, cin ', &
448                              cina(igout), cinb(igout), cin(igout)
449  ! -------------------------------------------------------------
450  ! --Update buoyancies to account for Ale
451  ! -------------------------------------------------------------
452
453  CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, &
454    tvp, buoy)
455  IF (prt_level>=20) PRINT *, 'cv3p2_closure after cv3_buoy'
456
457  ! -------------------------------------------------------------
458  ! -- Calculate convective available potential energy (cape),
459  ! -- vertical velocity (w), fractional area covered by
460  ! -- undilute updraft (sig), and updraft mass flux (m)
461  ! -------------------------------------------------------------
462
463  DO il = 1, ncum
464    cape(il) = 0.0
465    dtminmax(il) = -100.
466  END DO
467
468  ! compute dtmin (minimum buoyancy between ICB and given level k):
469
470  DO k = 1, nl
471    DO il = 1, ncum
472      dtmin(il, k) = 100.0
473    END DO
474  END DO
475
476  DO k = 1, nl
477    DO j = minorig, nl
478      DO il = 1, ncum
479        IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (j>=icb(il)) &
480                             .AND. (j<=(k-1))) THEN
481          dtmin(il, k) = amin1(dtmin(il,k), buoy(il,j))
482        END IF
483      END DO
484    END DO
485  END DO
486!jyg<
487!  Store maximum of dtmin
488!  C est pas terrible d avoir ce test sur Ale+Cin encore une fois ici.
489!                      A REVOIR !
490  DO k = 1, nl
491    DO il = 1, ncum
492      IF (k>=(icb(il)+1) .AND. k<=inb(il) .AND. ale(il)+cin(il)>0.) THEN
493        dtminmax(il) = max(dtmin(il,k), dtminmax(il))
494      ENDIF
495    END DO
496  END DO
497
498!    prevent convection when ale+cin <= 0
499  DO k = 1, nl
500    DO il = 1, ncum
501      IF (k>=(icb(il)+1) .AND. k<=inb(il)) THEN
502        dtmin(il,k) = min(dtmin(il,k), dtminmax(il))
503      ENDIF
504    END DO
505  END DO
506!>jyg
507
508  IF (prt_level >= 20) THEN
509    print *,'cv3p2_closure: dtmin ', (k, dtmin(igout,k), k=1,nl)
510    print *,'cv3p2_closure: dtminmax ', dtminmax(igout)
511  ENDIF
512
513  ! the interval on which cape is computed starts at pbase :
514
515  DO k = 1, nl
516    DO il = 1, ncum
517
518      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
519
520        IF (iflag_mix_adiab==1) THEN
521!CR:computation of cape from LCL: keep flag or to modify in all cases?
522        deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k))
523        ELSE
524        deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
525        ENDIF
526        cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1)
527        cape(il) = amax1(0.0, cape(il))
528        sigold(il, k) = sig(il, k)
529
530
531        ! jyg       Coefficient coefmix limits convection to levels where a
532        ! sufficient
533        ! fraction of mixed draughts are ascending.
534        siglim(il, k) = coefmix(il, k)*alpha1*dtmin(il, k)*abs(dtmin(il,k))
535        siglim(il, k) = amax1(siglim(il,k), 0.0)
536        siglim(il, k) = amin1(siglim(il,k), 0.01)
537        ! c         fac=AMIN1(((dtcrit-dtmin(il,k))/dtcrit),1.0)
538        fac = 1.
539        wlim(il, k) = fac*sqrt(cape(il))
540        amu = siglim(il, k)*wlim(il, k)
541!!        rhodp(il,k) = 0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k) !cor jyg : computed earlier
542        mlim(il, k) = amu*rhodp(il,k)
543        ! PRINT*, 'siglim ', k,siglim(1,k)
544      END IF
545
546    END DO
547  END DO
548  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 600'
549
550  DO il = 1, ncum
551    ! IM beg
552    IF (prt_level>=20) THEN
553      PRINT *, 'cv3p2_closure il icb mlim ph ph+1 ph+2', il, icb(il), &
554        mlim(il, icb(il)+1), ph(il, icb(il)), ph(il, icb(il)+1), &
555        ph(il, icb(il)+2)
556    END IF
557
558    IF (icb(il)+1<=inb(il)) THEN
559      ! IM end
560      mlim(il, icb(il)) = 0.5*mlim(il,icb(il)+1)*(ph(il,icb(il))-ph(il,icb(il)+1))/ &
561                                               (ph(il,icb(il)+1)-ph(il,icb(il)+2))
562      ! IM beg
563    END IF !(icb(il.le.inb(il))) THEN
564    ! IM end
565  END DO
566  IF (prt_level>=20) PRINT *, 'cv3p2_closure apres 700'
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)<ph(il,k)).AND.(plfc(il)>=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>0.).AND.(plfc(il)>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      ! END IF
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==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)>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) <= 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  ! c 4. Introduce a correcting factor for coef, in order to obtain an
844  ! effective
845  ! c    sigdz larger in the present case (using cv3p2_closure) than in the
846  ! old
847  ! c    closure (using cv3_closure).
848  IF (1==0) THEN
849    DO il = 1, ncum
850      ! c      coef(il) = 2.*coef(il)
851      coef(il) = 5.*coef(il)
852    END DO
853    ! version CVS du ..2008
854  ELSE
855    IF (iflag_cvl_sigd==0) THEN
856      ! test pour verifier qu on fait la meme chose qu avant: sid constant
857      coef(1:ncum) = 1.
858    ELSE
859      coef(1:ncum) = min(2.*coef(1:ncum), 5.)
860      coef(1:ncum) = max(2.*coef(1:ncum), 0.2)
861    END IF
862  END IF
863
864  IF (prt_level>=20) PRINT *, 'cv3p2_closure FIN'
865
866END SUBROUTINE cv3p2_closure
867
868
Note: See TracBrowser for help on using the repository browser.