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

Last change on this file since 5304 was 5304, checked in by abarral, 22 hours ago

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