source: LMDZ6/trunk/libf/phylmd/cv3p1_closure.f90 @ 5276

Last change on this file since 5276 was 5276, checked in by abarral, 7 hours ago

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