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

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

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