source: LMDZ6/branches/LMDZ_DECOUPLE/libf/phylmd/cv3p1_closure.F90 @ 5321

Last change on this file since 5321 was 3671, checked in by jyg, 5 years ago

Bug fix in cv3p1_closure,F90 and cv3p2_closure.F90
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 22.0 KB
Line 
1
2! $Id: cv3p1_closure.F90 3671 2020-04-29 13:48:22Z 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  INTEGER, DIMENSION (nloc), INTENT (INOUT)          :: iflag
44  REAL, DIMENSION (nloc, nd), INTENT (INOUT)         :: sig, w0
45  REAL, DIMENSION (nloc), INTENT (INOUT)             :: ptop2
46
47  ! output:
48  REAL, DIMENSION (nloc), INTENT (OUT)               :: cape, cin
49  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: m
50  REAL, DIMENSION (nloc), INTENT (OUT)               :: plim1, plim2
51  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: asupmax
52  REAL, DIMENSION (nloc), INTENT (OUT)               :: supmax0
53  REAL, DIMENSION (nloc), INTENT (OUT)               :: asupmaxmin
54  REAL, DIMENSION (nloc), INTENT (OUT)               :: cbmf, plfc
55  REAL, DIMENSION (nloc), INTENT (OUT)               :: wbeff
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, nd          !jyg: initialization up to nd
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,nd           !jyg: initialization up to nd
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 (NINT(flag_wb)==0) THEN
540        wbeff(il) = wbmax
541      ELSE IF (NINT(flag_wb)==1) THEN
542        wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
543      ELSE IF (NINT(flag_wb)==2) THEN
544        wbeff(il) = wbmax*(0.01*(ph(il,1)-plfc(il)))**2
545      ELSE ! Option provisoire ou le iflag_wb/10 est considere comme une vitesse
546        wbeff(il) = flag_wb*0.01+wbmax/(1.+500./(ph(il,1)-plfc(il)))
547      END IF
548    END IF
549  END DO
550
551!CR:Compute k at plfc
552  DO il=1,ncum
553           klfc(il)=nl
554  ENDDO
555  DO k=1,nl
556     DO il=1,ncum
557        if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
558           klfc(il)=k
559        endif
560     ENDDO
561  ENDDO
562!RC
563
564  DO il = 1, ncum
565    ! jyg    Modification du coef de wb*wb pour conformite avec papier Wake
566    ! c       cbmf1(il) = alp2(il)/(0.5*wb*wb-Cin(il))
567    cbmf1(il) = alp2(il)/(2.*wbeff(il)*wbeff(il)-cin(il))
568!CR: Add large-scale component to the mass-flux
569!encore connu sous le nom "Experience du tube de dentifrice"
570    if ((coef_clos_ls.gt.0.).and.(plfc(il).gt.0.)) then
571       cbmf1(il) = cbmf1(il) - coef_clos_ls*min(0.,1./RG*omega(il,klfc(il)))
572    endif
573!RC
574    IF (cbmf1(il)==0 .AND. alp2(il)/=0.) THEN
575      WRITE (lunout, *) 'cv3p1_closure cbmf1=0 and alp NE 0 il alp2 alp cin ' &
576        , il, alp2(il), alp(il), cin(il)
577      abort_message = ''
578      CALL abort_physic(modname, abort_message, 1)
579    END IF
580    cbmfmax(il) = sigmax*wb2(il)*100.*p(il, icb(il))/(rrd*tv(il,icb(il)))
581  END DO
582
583  DO il = 1, ncum
584    IF (cbmflim(il)>1.E-6) THEN
585      ! ATTENTION TEST CR
586      ! if (cbmfmax(il).lt.1.e-12) then
587      cbmf(il) = min(cbmf1(il), cbmfmax(il))
588      ! else
589      ! cbmf(il) = cbmf1(il)
590      ! endif
591      ! print*,'cbmf',cbmf1(il),cbmfmax(il)
592    END IF
593  END DO
594  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim_testCR'
595
596  ! c 2. Compute coefficient and apply correction
597
598  DO il = 1, ncum
599    coef(il) = (cbmf(il)+1.E-10)/(cbmflim(il)+1.E-10)
600  END DO
601  IF (prt_level>=20) PRINT *, 'cv3p1_param apres coef_plantePLUS'
602
603  DO k = 1, nl
604    DO il = 1, ncum
605      IF (k>=icb(il)+1 .AND. k<=inb(il)) THEN
606        amu = beta*sig(il, k)*w0(il, k) + (1.-beta)*coef(il)*siglim(il, k)* &
607          wlim(il, k)
608        w0(il, k) = wlim(il, k)
609        w0(il, k) = max(w0(il,k), 1.E-10)
610        sig(il, k) = amu/w0(il, k)
611        sig(il, k) = min(sig(il,k), 1.)
612        ! c         amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k)
613        m(il, k) = amu*0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
614      END IF
615    END DO
616  END DO
617  ! jyg2
618  DO il = 1, ncum
619    w0(il, icb(il)) = 0.5*w0(il, icb(il)+1)
620    m(il, icb(il)) = 0.5*m(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb(il)+1))/ &
621      (ph(il,icb(il)+1)-ph(il,icb(il)+2))
622    sig(il, icb(il)) = sig(il, icb(il)+1)
623    sig(il, icb(il)-1) = sig(il, icb(il))
624  END DO
625  IF (prt_level>=20) PRINT *, 'cv3p1_param apres w0_sig_M'
626
627!CR: new erosion of adiabatic ascent: modification of m
628!computation of the sum of ascending fluxes
629  IF (iflag_mix_adiab.eq.1) THEN
630
631!Verification sum(me)=sum(m)
632  DO k = 1,nd                         !jyg: initialization up to nd
633    DO il = 1, ncum
634       md(il,k)=0.
635       med(il,k)=0.
636    ENDDO
637  ENDDO
638
639  DO k = nl,1,-1
640    DO il = 1, ncum
641           md(il,k)=md(il,k+1)+m(il,k+1)
642    ENDDO
643  ENDDO
644
645  DO k = nl,1,-1
646    DO il = 1, ncum
647        IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
648           mad(il,k)=mad(il,k+1)+m(il,k+1)
649        ENDIF
650!        print*,"mad",il,k,mad(il,k)
651    ENDDO
652  ENDDO
653
654!CR: erosion of each adiabatic ascent during its ascent
655
656!Computation of erosion coefficient beta_coef
657  DO k = 1, nl
658    DO il = 1, ncum
659       IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN     
660!          print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
661          beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
662       ELSE
663          beta_coef(il,k)=0.
664       ENDIF
665    ENDDO
666  ENDDO
667
668!  print*,"apres beta_coef"
669
670  DO k = 1, nl
671    DO il = 1, ncum
672
673      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
674
675!        print*,"dz",il,k,tv(il, k-1)
676        dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG)
677        betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz)
678!        betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz)
679!        print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
680        dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG)
681!        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))
682        me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k))
683!        print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz   
684     
685      END IF
686       
687!Modification of m
688      m(il,k)=me(il,k)
689    END DO
690  END DO
691 
692!  DO il = 1, ncum
693!     dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG)
694!     m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) &
695!                  /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il))
696!     print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
697!  ENDDO
698
699!Verification sum(me)=sum(m)
700  DO k = nl,1,-1
701    DO il = 1, ncum
702           med(il,k)=med(il,k+1)+m(il,k+1)
703!           print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
704    ENDDO
705  ENDDO
706
707
708  ENDIF !(iflag_mix_adiab)
709!RC
710
711
712
713  ! c 3. Compute final cloud base mass flux and set iflag to 3 if
714  ! c    cloud base mass flux is exceedingly small and is decreasing (i.e. if
715  ! c    the final mass flux (cbmflast) is greater than the target mass flux
716  ! c    (cbmf)).
717
718  DO il = 1, ncum
719    cbmflast(il) = 0.
720  END DO
721
722  DO k = 1, nl
723    DO il = 1, ncum
724      IF (k>=icb(il) .AND. k<=inb(il)) THEN
725          !IMpropo??      IF ((k.ge.(icb(il)+1)).and.(k.le.inb(il))) THEN
726        cbmflast(il) = cbmflast(il) + m(il, k)
727      END IF
728    END DO
729  END DO
730
731  DO il = 1, ncum
732    IF (cbmflast(il)<1.E-6 .AND. cbmflast(il)>=cbmf(il)) THEN
733      iflag(il) = 3
734    END IF
735  END DO
736
737  DO k = 1, nl
738    DO il = 1, ncum
739      IF (iflag(il)>=3) THEN
740        m(il, k) = 0.
741        sig(il, k) = 0.
742        w0(il, k) = 0.
743      END IF
744    END DO
745  END DO
746  IF (prt_level>=20) PRINT *, 'cv3p1_param apres iflag'
747
748  ! c 4. Introduce a correcting factor for coef, in order to obtain an
749  ! effective
750  ! c    sigdz larger in the present case (using cv3p1_closure) than in the
751  ! old
752  ! c    closure (using cv3_closure).
753  IF (1==0) THEN
754    DO il = 1, ncum
755      ! c      coef(il) = 2.*coef(il)
756      coef(il) = 5.*coef(il)
757    END DO
758    ! version CVS du ..2008
759  ELSE
760    IF (iflag_cvl_sigd==0) THEN
761      ! test pour verifier qu on fait la meme chose qu avant: sid constant
762      coef(1:ncum) = 1.
763    ELSE
764      coef(1:ncum) = min(2.*coef(1:ncum), 5.)
765      coef(1:ncum) = max(2.*coef(1:ncum), 0.2)
766    END IF
767  END IF
768
769  IF (prt_level>=20) PRINT *, 'cv3p1_param FIN'
770  RETURN
771END SUBROUTINE cv3p1_closure
772
773
Note: See TracBrowser for help on using the repository browser.