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

Last change on this file was 5274, checked in by abarral, 55 minutes ago

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