source: LMDZ6/trunk/libf/phylmd/cv3p2_closure.F90 @ 4013

Last change on this file since 4013 was 3671, checked in by jyg, 5 years ago

Bug fix in cv3p1_closure,F90 and cv3p2_closure.F90
File size: 26.3 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
20  USE print_control_mod, ONLY: prt_level, lunout
21  IMPLICIT NONE
22
23  include "cvthermo.h"
24  include "cv3param.h"
[3496]25  include "cvflag.h"
[2374]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:
[3671]43  INTEGER, DIMENSION (nloc), INTENT (INOUT)          :: iflag
[2374]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
[2420]62  REAL                                               :: dz
[2374]63  REAL                                               :: pbmxup
64  REAL, DIMENSION (nloc, nd)                         :: dtmin, sigold
65  REAL, DIMENSION (nloc, nd)                         :: coefmix
[2398]66  REAL, DIMENSION (nloc)                             :: dtminmax
[2374]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
[2420]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
[2458]92!jyg<
93! coef_peel is now in the common cv3_param
94!!  REAL                                               :: coef_peel
95!!  PARAMETER (coef_peel=0.25)
96!>jyg
[2420]97
[2374]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
[2502]133!CR+jyg: initializations (up to nd) for erosion of adiabatic ascent and of m and wlim
134  DO k = 1,nd
[2420]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.
[2502]140        m(il, k) = 0.0
[2420]141    ENDDO
142  ENDDO
143
[2374]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
[2398]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
[2374]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
[2398]206  ENDIF  ! (ok_convstop)
207!>jyg
[2374]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
369425 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
[2398]464    dtminmax(il) = -100.
[2374]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
[2398]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
[2374]496!
[2398]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!
[2374]507  IF (prt_level >= 20) THEN
508    print *,'cv3p2_closure: dtmin ', (k, dtmin(igout,k), k=1,nl)
[2398]509    print *,'cv3p2_closure: dtminmax ', dtminmax(igout)
[2374]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
[2420]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
[2374]523        deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
[2420]524        ENDIF
[2374]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
[3571]611      IF (NINT(flag_wb)==0) THEN
[2374]612        wbeff(il) = wbmax
[3571]613      ELSE IF (NINT(flag_wb)==1) THEN
[2374]614        wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
[3571]615      ELSE IF (NINT(flag_wb)==2) THEN
[2374]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
[2398]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
[2374]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
[2398]674  ENDIF  !(OK_intermittent)
[2374]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
[2420]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)
[2502]714  DO k = 1,nd
[2420]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
[2374]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)).
[2398]797  ! c    If(ok_convstop): set iflag to 4 if no positive buoyancy has been met
[2374]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
[2398]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
[2374]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
[2398]837  ENDIF ! (ok_convstop)
[2374]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
867END SUBROUTINE cv3p2_closure
868
869
Note: See TracBrowser for help on using the repository browser.