source: LMDZ6/branches/Amaury_dev/libf/phylmd/cv3p1_closure.F90 @ 5157

Last change on this file since 5157 was 5144, checked in by abarral, 8 weeks ago

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