source: LMDZ5/trunk/libf/phylmd/cv3p2_closure.F90 @ 3567

Last change on this file since 3567 was 2502, checked in by jyg, 9 years ago

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