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

Last change on this file since 2454 was 2420, checked in by crio, 9 years ago

Nouvelle option d'epluchage de l'ascendance adiabatique dans le schema d'Emanuel: epluchage fonction de B/w2 au lieu de w. S'active avec iflag_mix_adiab=1 (valeur par defaut iflag_mix_adiab=0). Fonctionne avec iflag_mix=0 et iflag_mix=1.
Correction de bugs dans le schema de convection pour le calcul de inb, cape et buoy (sous le meme flag pour l'instant).
New option for the erosion of the adiabatic ascent in the Emanuel scheme: erosion function of B/w2 instead of w. Activated by iflag_mix_adiab=1 (standard value iflag_mix_adiab=0). Should work with iflag_mix=0 and iflag_mix=1.
Various bug corrections in the convection scheme for the computation of inb, cape, buoy (protected by the same flag for now).

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