source: LMDZ6/trunk/libf/phylmd/cv3p2_closure.f90 @ 5283

Last change on this file since 5283 was 5283, checked in by abarral, 3 days ago

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