source: LMDZ5/trunk/libf/phylmd/cv3p1_closure.F90 @ 2253

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

1/ Introduction of two variables in the ".def" files: (i) cvl_sig2feed is
the top of the convective feeding layer in sigma coordinates (D=0.97);
(ii) cvl_comp_threshold is the threshold fraction of convective points
below which compression occurs (D=1.).
2/ Corrections of various bugs revealed by the changes in compression:

  • correct bugs in cv3a_uncompress.F90 for 3 fields used for convective

scavenging.

  • add a reset to zero of "sig" and "w0" for non-convective points

(cva_driver.F90).

  • in cv3_routines.F90, correct bounds of a few loops in cv3_undilute2,

correct the reset of the no-convection counter in cv3_yield.

  • in phys_output_write_mod.F90, correct output of wdtrainA and wdtrainM.

3/ Improve declarations in various subroutines.

Modified files:

conema3.h
cv3param.h
cv3p1_closure.F90
conf_phys_m.F90
cv3a_compress.F90
phys_output_write_mod.F90
cv3_routines.F90
concvl.F90
cva_driver.F90
cv3a_uncompress.F90
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.5 KB
Line 
1
2! $Id: cv3p1_closure.F90 2253 2015-03-30 09:08:45Z jyg $
3
4SUBROUTINE cv3p1_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
5    tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, &
6    iflag, coef, plim1, plim2, asupmax, supmax0, asupmaxmin, cbmf, plfc, &
7    wbeff)
8
9
10  ! **************************************************************
11  ! *
12  ! CV3P1_CLOSURE                                               *
13  ! Ale & Alp Closure of Convect3              *
14  ! *
15  ! written by   :   Kerry Emanuel                              *
16  ! vectorization:   S. Bony                                    *
17  ! modified by :    Jean-Yves Grandpeix, 18/06/2003, 19.32.10  *
18  ! Julie Frohwirth,     14/10/2005  17.44.22  *
19  ! **************************************************************
20
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  include "iniprint.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:
43  REAL, DIMENSION (nloc, nd), INTENT (INOUT)         :: sig, w0
44  REAL, DIMENSION (nloc), INTENT (INOUT)             :: ptop2
45
46  ! output:
47  REAL, DIMENSION (nloc), INTENT (OUT)               :: cape, cin
48  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: m
49  REAL, DIMENSION (nloc), INTENT (OUT)               :: plim1, plim2
50  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: asupmax
51  REAL, DIMENSION (nloc), INTENT (OUT)               :: supmax0
52  REAL, DIMENSION (nloc), INTENT (OUT)               :: asupmaxmin
53  REAL, DIMENSION (nloc), INTENT (OUT)               :: cbmf, plfc
54  REAL, DIMENSION (nloc), INTENT (OUT)               :: wbeff
55  INTEGER, DIMENSION (nloc), INTENT (OUT)            :: iflag
56
57  ! local variables:
58  INTEGER il, i, j, k, icbmax, i0(nloc), klfc(nloc)
59  REAL deltap, fac, w, amu
60  REAL rhodp
61  REAL pbmxup
62  REAL dtmin(nloc, nd), sigold(nloc, nd)
63  REAL coefmix(nloc, nd)
64  REAL pzero(nloc), ptop2old(nloc)
65  REAL cina(nloc), cinb(nloc)
66  INTEGER ibeg(nloc)
67  INTEGER nsupmax(nloc)
68  REAL supcrit, temp(nloc, nd)
69  REAL p1(nloc), pmin(nloc)
70  REAL asupmax0(nloc)
71  LOGICAL ok(nloc)
72  REAL siglim(nloc, nd), wlim(nloc, nd), mlim(nloc, nd)
73  REAL wb2(nloc)
74  REAL cbmflim(nloc), cbmf1(nloc), cbmfmax(nloc)
75  REAL cbmflast(nloc)
76  REAL coef(nloc)
77  REAL xp(nloc), xq(nloc), xr(nloc), discr(nloc), b3(nloc), b4(nloc)
78  REAL theta(nloc), bb(nloc)
79  REAL term1, term2, term3
80  REAL alp2(nloc) ! Alp with offset
81
82  REAL sigmax
83  PARAMETER (sigmax=0.1)
84
85  CHARACTER (LEN=20) :: modname = 'cv3p1_closure'
86  CHARACTER (LEN=80) :: abort_message
87
88  ! print *,' -> cv3p1_closure, Ale ',ale(1)
89
90
91  ! -------------------------------------------------------
92  ! -- Initialization
93  ! -------------------------------------------------------
94
95
96  DO il = 1, ncum
97    alp2(il) = max(alp(il), 1.E-5)
98    ! IM
99    alp2(il) = max(alp(il), 1.E-12)
100  END DO
101
102  pbmxup = 50. ! PBMXUP+PBCRIT = cloud depth above which mixed updraughts
103  ! exist (if any)
104
105  IF (prt_level>=20) PRINT *, 'cv3p1_param nloc ncum nd icb inb nl', nloc, &
106    ncum, nd, icb(nloc), inb(nloc), nl
107  DO k = 1, nl
108    DO il = 1, ncum
109      m(il, k) = 0.0
110    END DO
111  END DO
112
113  ! -------------------------------------------------------
114  ! -- Reset sig(i) and w0(i) for i>inb and i<icb
115  ! -------------------------------------------------------
116
117  ! update sig and w0 above LNB:
118
119  DO k = 1, nl - 1
120    DO il = 1, ncum
121      IF ((inb(il)<(nl-1)) .AND. (k>=(inb(il)+1))) THEN
122        sig(il, k) = beta*sig(il, k) + 2.*alpha*buoy(il, inb(il))*abs(buoy(il &
123          ,inb(il)))
124        sig(il, k) = amax1(sig(il,k), 0.0)
125        w0(il, k) = beta*w0(il, k)
126      END IF
127    END DO
128  END DO
129
130  ! if(prt.level.GE.20) print*,'cv3p1_param apres 100'
131  ! compute icbmax:
132
133  icbmax = 2
134  DO il = 1, ncum
135    icbmax = max(icbmax, icb(il))
136  END DO
137  ! if(prt.level.GE.20) print*,'cv3p1_param apres 200'
138
139  ! update sig and w0 below cloud base:
140
141  DO k = 1, icbmax
142    DO il = 1, ncum
143      IF (k<=icb(il)) THEN
144        sig(il, k) = beta*sig(il, k) - 2.*alpha*buoy(il, icb(il))*buoy(il, &
145          icb(il))
146        sig(il, k) = amax1(sig(il,k), 0.0)
147        w0(il, k) = beta*w0(il, k)
148      END IF
149    END DO
150  END DO
151  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 300'
152  ! -------------------------------------------------------------
153  ! -- Reset fractional areas of updrafts and w0 at initial time
154  ! -- and after 10 time steps of no convection
155  ! -------------------------------------------------------------
156
157  DO k = 1, nl - 1
158    DO il = 1, ncum
159      IF (sig(il,nd)<1.5 .OR. sig(il,nd)>12.0) THEN
160        sig(il, k) = 0.0
161        w0(il, k) = 0.0
162      END IF
163    END DO
164  END DO
165  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 400'
166
167  ! -------------------------------------------------------------
168  ! jyg1
169  ! --  Calculate adiabatic ascent top pressure (ptop)
170  ! -------------------------------------------------------------
171
172
173  ! c 1. Start at first level where precipitations form
174  DO il = 1, ncum
175    pzero(il) = plcl(il) - pbcrit
176  END DO
177
178  ! c 2. Add offset
179  DO il = 1, ncum
180    pzero(il) = pzero(il) - pbmxup
181  END DO
182  DO il = 1, ncum
183    ptop2old(il) = ptop2(il)
184  END DO
185
186  DO il = 1, ncum
187    ! CR:c est quoi ce 300??
188    p1(il) = pzero(il) - 300.
189  END DO
190
191  ! compute asupmax=abs(supmax) up to lnm+1
192
193  DO il = 1, ncum
194    ok(il) = .TRUE.
195    nsupmax(il) = inb(il)
196  END DO
197
198  DO i = 1, nl
199    DO il = 1, ncum
200      IF (i>icb(il) .AND. i<=inb(il)) THEN
201        IF (p(il,i)<=pzero(il) .AND. supmax(il,i)<0 .AND. ok(il)) THEN
202          nsupmax(il) = i
203          ok(il) = .FALSE.
204        END IF ! end IF (P(i) ...  )
205      END IF ! end IF (icb+1 le i le inb)
206    END DO
207  END DO
208
209  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 2.'
210  DO i = 1, nl
211    DO il = 1, ncum
212      asupmax(il, i) = abs(supmax(il,i))
213    END DO
214  END DO
215
216
217  DO il = 1, ncum
218    asupmaxmin(il) = 10.
219    pmin(il) = 100.
220    ! IM ??
221    asupmax0(il) = 0.
222  END DO
223
224  ! c 3.  Compute in which level is Pzero
225
226  ! IM bug      i0 = 18
227  DO il = 1, ncum
228    i0(il) = nl
229  END DO
230
231  DO i = 1, nl
232    DO il = 1, ncum
233      IF (i>icb(il) .AND. i<=inb(il)) THEN
234        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
235          IF (pzero(il)>p(il,i) .AND. pzero(il)<p(il,i-1)) THEN
236            i0(il) = i
237          END IF
238        END IF
239      END IF
240    END DO
241  END DO
242  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 3.'
243
244  ! c 4.  Compute asupmax at Pzero
245
246  DO i = 1, nl
247    DO il = 1, ncum
248      IF (i>icb(il) .AND. i<=inb(il)) THEN
249        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
250          asupmax0(il) = ((pzero(il)-p(il,i0(il)-1))*asupmax(il,i0(il))-( &
251            pzero(il)-p(il,i0(il)))*asupmax(il,i0(il)-1))/(p(il,i0(il))-p(il, &
252            i0(il)-1))
253        END IF
254      END IF
255    END DO
256  END DO
257
258
259  DO i = 1, nl
260    DO il = 1, ncum
261      IF (p(il,i)==pzero(il)) THEN
262        asupmax(i, il) = asupmax0(il)
263      END IF
264    END DO
265  END DO
266  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 4.'
267
268  ! c 5. Compute asupmaxmin, minimum of asupmax
269
270  DO i = 1, nl
271    DO il = 1, ncum
272      IF (i>icb(il) .AND. i<=inb(il)) THEN
273        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
274          IF (asupmax(il,i)<asupmaxmin(il)) THEN
275            asupmaxmin(il) = asupmax(il, i)
276            pmin(il) = p(il, i)
277          END IF
278        END IF
279      END IF
280    END DO
281  END DO
282
283  DO il = 1, ncum
284    ! IM
285    IF (prt_level>=20) THEN
286      PRINT *, 'cv3p1_closure il asupmax0 asupmaxmin', il, asupmax0(il), &
287        asupmaxmin(il), pzero(il), pmin(il)
288    END IF
289    IF (asupmax0(il)<asupmaxmin(il)) THEN
290      asupmaxmin(il) = asupmax0(il)
291      pmin(il) = pzero(il)
292    END IF
293  END DO
294  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 5.'
295
296
297  ! Compute Supmax at Pzero
298
299  DO i = 1, nl
300    DO il = 1, ncum
301      IF (i>icb(il) .AND. i<=inb(il)) THEN
302        IF (p(il,i)<=pzero(il)) THEN
303          supmax0(il) = ((p(il,i)-pzero(il))*asupmax(il,i-1)-(p(il, &
304            i-1)-pzero(il))*asupmax(il,i))/(p(il,i)-p(il,i-1))
305          GO TO 425
306        END IF ! end IF (P(i) ... )
307      END IF ! end IF (icb+1 le i le inb)
308    END DO
309  END DO
310
311425 CONTINUE
312  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 425.'
313
314  ! c 6. Calculate ptop2
315
316  DO il = 1, ncum
317    IF (asupmaxmin(il)<supcrit1) THEN
318      ptop2(il) = pmin(il)
319    END IF
320
321    IF (asupmaxmin(il)>supcrit1 .AND. asupmaxmin(il)<supcrit2) THEN
322      ptop2(il) = ptop2old(il)
323    END IF
324
325    IF (asupmaxmin(il)>supcrit2) THEN
326      ptop2(il) = ph(il, inb(il))
327    END IF
328  END DO
329
330  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 6.'
331
332  ! c 7. Compute multiplying factor for adiabatic updraught mass flux
333
334
335  IF (ok_inhib) THEN
336
337    DO i = 1, nl
338      DO il = 1, ncum
339        IF (i<=nl) THEN
340          coefmix(il, i) = (min(ptop2(il),ph(il,i))-ph(il,i))/(ph(il,i+1)-ph( &
341            il,i))
342          coefmix(il, i) = min(coefmix(il,i), 1.)
343        END IF
344      END DO
345    END DO
346
347
348  ELSE ! when inhibition is not taken into account, coefmix=1
349
350
351
352    DO i = 1, nl
353      DO il = 1, ncum
354        IF (i<=nl) THEN
355          coefmix(il, i) = 1.
356        END IF
357      END DO
358    END DO
359
360  END IF ! ok_inhib
361  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 7.'
362  ! -------------------------------------------------------------------
363  ! -------------------------------------------------------------------
364
365
366  ! jyg2
367
368  ! ==========================================================================
369
370
371  ! -------------------------------------------------------------
372  ! -- Calculate convective inhibition (CIN)
373  ! -------------------------------------------------------------
374
375  ! do i=1,nloc
376  ! print*,'avant cine p',pbase(i),plcl(i)
377  ! enddo
378  ! do j=1,nd
379  ! do i=1,nloc
380  ! print*,'avant cine t',tv(i),tvp(i)
381  ! enddo
382  ! enddo
383  CALL cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, &
384    cinb, plfc)
385
386  DO il = 1, ncum
387    cin(il) = cina(il) + cinb(il)
388  END DO
389  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_cine'
390  ! -------------------------------------------------------------
391  ! --Update buoyancies to account for Ale
392  ! -------------------------------------------------------------
393
394  CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, &
395    tvp, buoy)
396  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_buoy'
397
398  ! -------------------------------------------------------------
399  ! -- Calculate convective available potential energy (cape),
400  ! -- vertical velocity (w), fractional area covered by
401  ! -- undilute updraft (sig), and updraft mass flux (m)
402  ! -------------------------------------------------------------
403
404  DO il = 1, ncum
405    cape(il) = 0.0
406  END DO
407
408  ! compute dtmin (minimum buoyancy between ICB and given level k):
409
410  DO k = 1, nl
411    DO il = 1, ncum
412      dtmin(il, k) = 100.0
413    END DO
414  END DO
415
416  DO k = 1, nl
417    DO j = minorig, nl
418      DO il = 1, ncum
419        IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (j>=icb(il)) .AND. (j<= &
420            (k-1))) THEN
421          dtmin(il, k) = amin1(dtmin(il,k), buoy(il,j))
422        END IF
423      END DO
424    END DO
425  END DO
426
427  ! the interval on which cape is computed starts at pbase :
428
429  DO k = 1, nl
430    DO il = 1, ncum
431
432      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
433
434        deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
435        cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1)
436        cape(il) = amax1(0.0, cape(il))
437        sigold(il, k) = sig(il, k)
438
439
440        ! jyg       Coefficient coefmix limits convection to levels where a
441        ! sufficient
442        ! fraction of mixed draughts are ascending.
443        siglim(il, k) = coefmix(il, k)*alpha1*dtmin(il, k)*abs(dtmin(il,k))
444        siglim(il, k) = amax1(siglim(il,k), 0.0)
445        siglim(il, k) = amin1(siglim(il,k), 0.01)
446        ! c         fac=AMIN1(((dtcrit-dtmin(il,k))/dtcrit),1.0)
447        fac = 1.
448        wlim(il, k) = fac*sqrt(cape(il))
449        amu = siglim(il, k)*wlim(il, k)
450        rhodp = 0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
451        mlim(il, k) = amu*rhodp
452        ! print*, 'siglim ', k,siglim(1,k)
453      END IF
454
455    END DO
456  END DO
457  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 600'
458
459  DO il = 1, ncum
460    ! IM beg
461    IF (prt_level>=20) THEN
462      PRINT *, 'cv3p1_closure il icb mlim ph ph+1 ph+2', il, icb(il), &
463        mlim(il, icb(il)+1), ph(il, icb(il)), ph(il, icb(il)+1), &
464        ph(il, icb(il)+2)
465    END IF
466
467    IF (icb(il)+1<=inb(il)) THEN
468      ! IM end
469      mlim(il, icb(il)) = 0.5*mlim(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb( &
470        il)+1))/(ph(il,icb(il)+1)-ph(il,icb(il)+2))
471      ! IM beg
472    END IF !(icb(il.le.inb(il))) then
473    ! IM end
474  END DO
475  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 700'
476
477  ! jyg1
478  ! ------------------------------------------------------------------------
479  ! c     Correct mass fluxes so that power used to overcome CIN does not
480  ! c     exceed Power Available for Lifting (PAL).
481  ! ------------------------------------------------------------------------
482
483  DO il = 1, ncum
484    cbmflim(il) = 0.
485    cbmf(il) = 0.
486  END DO
487
488  ! c 1. Compute cloud base mass flux of elementary system (Cbmf0=Cbmflim)
489
490  DO k = 1, nl
491    DO il = 1, ncum
492      ! old       IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
493      ! IM        IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
494      IF (k>=icb(il) .AND. k<=inb(il) & !cor jyg
495          .AND. icb(il)+1<=inb(il)) THEN !cor jyg
496        cbmflim(il) = cbmflim(il) + mlim(il, k)
497      END IF
498    END DO
499  END DO
500  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim'
501
502  ! 1.5 Compute cloud base mass flux given by Alp closure (Cbmf1), maximum
503  !     allowed mass flux (Cbmfmax) and final target mass flux (Cbmf)
504  !     Cbmf is set to zero if Cbmflim (the mass flux of elementary cloud)
505  !     is exceedingly small.
506
507  DO il = 1, ncum
508    wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.))
509  END DO
510
511  DO il = 1, ncum
512    IF (plfc(il)<100.) THEN
513      ! This is an irealistic value for plfc => no calculation of wbeff
514      wbeff(il) = 100.1
515    ELSE
516      ! Calculate wbeff
517      IF (flag_wb==0) THEN
518        wbeff(il) = wbmax
519      ELSE IF (flag_wb==1) THEN
520        wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
521      ELSE IF (flag_wb==2) THEN
522        wbeff(il) = wbmax*(0.01*(ph(il,1)-plfc(il)))**2
523      END IF
524    END IF
525  END DO
526
527!CR:Compute k at plfc
528  DO il=1,ncum
529           klfc(il)=nl
530  ENDDO
531  DO k=1,nl
532     DO il=1,ncum
533        if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
534           klfc(il)=k
535        endif
536     ENDDO
537  ENDDO
538!RC
539
540  DO il = 1, ncum
541    ! jyg    Modification du coef de wb*wb pour conformite avec papier Wake
542    ! c       cbmf1(il) = alp2(il)/(0.5*wb*wb-Cin(il))
543    cbmf1(il) = alp2(il)/(2.*wbeff(il)*wbeff(il)-cin(il))
544!CR: Add large-scale component to the mass-flux
545!encore connu sous le nom "Experience du tube de dentifrice"
546    if ((coef_clos_ls.gt.0.).and.(plfc(il).gt.0.)) then
547       cbmf1(il) = cbmf1(il) - coef_clos_ls*min(0.,1./RG*omega(il,klfc(il)))
548    endif
549!RC
550    IF (cbmf1(il)==0 .AND. alp2(il)/=0.) THEN
551      WRITE (lunout, *) 'cv3p1_closure cbmf1=0 and alp NE 0 il alp2 alp cin ' &
552        , il, alp2(il), alp(il), cin(il)
553      abort_message = ''
554      CALL abort_gcm(modname, abort_message, 1)
555    END IF
556    cbmfmax(il) = sigmax*wb2(il)*100.*p(il, icb(il))/(rrd*tv(il,icb(il)))
557  END DO
558
559  DO il = 1, ncum
560    IF (cbmflim(il)>1.E-6) THEN
561      ! ATTENTION TEST CR
562      ! if (cbmfmax(il).lt.1.e-12) then
563      cbmf(il) = min(cbmf1(il), cbmfmax(il))
564      ! else
565      ! cbmf(il) = cbmf1(il)
566      ! endif
567      ! print*,'cbmf',cbmf1(il),cbmfmax(il)
568    END IF
569  END DO
570  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim_testCR'
571
572  ! c 2. Compute coefficient and apply correction
573
574  DO il = 1, ncum
575    coef(il) = (cbmf(il)+1.E-10)/(cbmflim(il)+1.E-10)
576  END DO
577  IF (prt_level>=20) PRINT *, 'cv3p1_param apres coef_plantePLUS'
578
579  DO k = 1, nl
580    DO il = 1, ncum
581      IF (k>=icb(il)+1 .AND. k<=inb(il)) THEN
582        amu = beta*sig(il, k)*w0(il, k) + (1.-beta)*coef(il)*siglim(il, k)* &
583          wlim(il, k)
584        w0(il, k) = wlim(il, k)
585        w0(il, k) = max(w0(il,k), 1.E-10)
586        sig(il, k) = amu/w0(il, k)
587        sig(il, k) = min(sig(il,k), 1.)
588        ! c         amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k)
589        m(il, k) = amu*0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
590      END IF
591    END DO
592  END DO
593  ! jyg2
594  DO il = 1, ncum
595    w0(il, icb(il)) = 0.5*w0(il, icb(il)+1)
596    m(il, icb(il)) = 0.5*m(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb(il)+1))/ &
597      (ph(il,icb(il)+1)-ph(il,icb(il)+2))
598    sig(il, icb(il)) = sig(il, icb(il)+1)
599    sig(il, icb(il)-1) = sig(il, icb(il))
600  END DO
601  IF (prt_level>=20) PRINT *, 'cv3p1_param apres w0_sig_M'
602
603  ! c 3. Compute final cloud base mass flux and set iflag to 3 if
604  ! c    cloud base mass flux is exceedingly small and is decreasing (i.e. if
605  ! c    the final mass flux (cbmflast) is greater than the target mass flux
606  ! c    (cbmf)).
607
608  DO il = 1, ncum
609    cbmflast(il) = 0.
610  END DO
611
612  DO k = 1, nl
613    DO il = 1, ncum
614      IF (k>=icb(il) .AND. k<=inb(il)) THEN
615          !IMpropo??      IF ((k.ge.(icb(il)+1)).and.(k.le.inb(il))) THEN
616        cbmflast(il) = cbmflast(il) + m(il, k)
617      END IF
618    END DO
619  END DO
620
621  DO il = 1, ncum
622    IF (cbmflast(il)<1.E-6 .AND. cbmflast(il)>=cbmf(il)) THEN
623      iflag(il) = 3
624    END IF
625  END DO
626
627  DO k = 1, nl
628    DO il = 1, ncum
629      IF (iflag(il)>=3) THEN
630        m(il, k) = 0.
631        sig(il, k) = 0.
632        w0(il, k) = 0.
633      END IF
634    END DO
635  END DO
636  IF (prt_level>=20) PRINT *, 'cv3p1_param apres iflag'
637
638  ! c 4. Introduce a correcting factor for coef, in order to obtain an
639  ! effective
640  ! c    sigdz larger in the present case (using cv3p1_closure) than in the
641  ! old
642  ! c    closure (using cv3_closure).
643  IF (1==0) THEN
644    DO il = 1, ncum
645      ! c      coef(il) = 2.*coef(il)
646      coef(il) = 5.*coef(il)
647    END DO
648    ! version CVS du ..2008
649  ELSE
650    IF (iflag_cvl_sigd==0) THEN
651      ! test pour verifier qu on fait la meme chose qu avant: sid constant
652      coef(1:ncum) = 1.
653    ELSE
654      coef(1:ncum) = min(2.*coef(1:ncum), 5.)
655      coef(1:ncum) = max(2.*coef(1:ncum), 0.2)
656    END IF
657  END IF
658
659  IF (prt_level>=20) PRINT *, 'cv3p1_param FIN'
660  RETURN
661END SUBROUTINE cv3p1_closure
662
663
Note: See TracBrowser for help on using the repository browser.