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

Last change on this file since 5833 was 5708, checked in by yann meurdesoif, 5 months ago

Convection GPU porting : suppress potential dependency between columns, may change results (cv3p1_closure)

YM
  • 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.5 KB
Line 
1
2! $Id: cv3p1_closure.f90 5708 2025-06-16 14:26:48Z rkazeroni $
3MODULE cv3p1_closure_mod
4PRIVATE
5PUBLIC cv3p1_closure
6CONTAINS
7
8SUBROUTINE cv3p1_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
9    tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, &
10    iflag, coef, coeftrue, plim1, plim2, asupmax, supmax0, asupmaxmin, &
11    cbmf, plfc, wbeff)
12
13
14  ! **************************************************************
15  ! *
16  ! CV3P1_CLOSURE                                               *
17  ! Ale & Alp Closure of Convect3              *
18  ! *
19  ! written by   :   Kerry Emanuel                              *
20  ! vectorization:   S. Bony                                    *
21  ! modified by :    Jean-Yves Grandpeix, 18/06/2003, 19.32.10  *
22  ! Julie Frohwirth,     14/10/2005  17.44.22  *
23  ! **************************************************************
24
25USE yomcst2_mod_h
26   USE lmdz_cv_ini, ONLY : pbcrit,wbmax,rrd,minorig,flag_wb,alpha,alpha1,beta,coef_peel,nl
27  USE conema3_mod_h
28  USE print_control_mod, ONLY: prt_level, lunout
29  USE yomcst_mod_h
30  USE cv3_cine_mod, ONLY : cv3_cine
31  USE cv3_buoy_mod, ONLY : cv3_buoy
32IMPLICIT NONE
33
34
35
36  ! input:
37  INTEGER, INTENT (IN)                               :: ncum, nd, nloc
38  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb
39  REAL, DIMENSION (nloc), INTENT (IN)                :: pbase, plcl
40  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: p
41  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
42  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: tv, tvp, buoy
43  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: supmax
44  LOGICAL, INTENT (IN)                               :: ok_inhib ! enable convection inhibition by dryness
45  REAL, DIMENSION (nloc), INTENT (IN)                :: ale, alp
46  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: omega
47
48  ! input/output:
49  INTEGER, DIMENSION (nloc), INTENT (INOUT)          :: iflag
50  REAL, DIMENSION (nloc, nd), INTENT (INOUT)         :: sig, w0
51  REAL, DIMENSION (nloc), INTENT (INOUT)             :: ptop2
52
53  ! output:
54  REAL, DIMENSION (nloc), INTENT (OUT)               :: cape, cin
55  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: m
56  REAL, DIMENSION (nloc), INTENT (OUT)               :: coef, coeftrue
57  REAL, DIMENSION (nloc), INTENT (OUT)               :: plim1, plim2
58  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: asupmax
59  REAL, DIMENSION (nloc), INTENT (OUT)               :: supmax0
60  REAL, DIMENSION (nloc), INTENT (OUT)               :: asupmaxmin
61  REAL, DIMENSION (nloc), INTENT (OUT)               :: cbmf, plfc
62  REAL, DIMENSION (nloc), INTENT (OUT)               :: wbeff
63
64  ! local variables:
65  INTEGER il, i, j, k, icbmax, i0(nloc), klfc(nloc)
66  REAL deltap, fac, w, amu
67  REAL rhodp, dz
68  REAL pbmxup
69  REAL dtmin(nloc, nd), sigold(nloc, nd)
70  REAL coefmix(nloc, nd)
71  REAL pzero(nloc), ptop2old(nloc)
72  REAL cina(nloc), cinb(nloc)
73  INTEGER ibeg(nloc)
74  INTEGER nsupmax(nloc)
75  REAL supcrit, temp(nloc, nd)
76  REAL p1(nloc), pmin(nloc)
77  REAL asupmax0(nloc)
78  LOGICAL ok(nloc)
79  REAL siglim(nloc, nd), wlim(nloc, nd), mlim(nloc, nd)
80  REAL wb2(nloc)
81  REAL cbmflim(nloc), cbmf1(nloc), cbmfmax(nloc)
82  REAL cbmflast(nloc)
83  REAL xp(nloc), xq(nloc), xr(nloc), discr(nloc), b3(nloc), b4(nloc)
84  REAL theta(nloc), bb(nloc)
85  REAL term1, term2, term3
86  REAL alp2(nloc) ! Alp with offset
87  !CR: variables for new erosion of adiabiatic ascent
88  REAL mad(nloc, nd), me(nloc, nd), betalim(nloc, nd), beta_coef(nloc, nd)
89  REAL med(nloc, nd), md(nloc,nd)
90!jyg<
91! coef_peel is now in the common cv3_param
92!!  REAL coef_peel
93!!  PARAMETER (coef_peel=0.25)
94!>jyg
95
96  REAL,PARAMETER :: sigmax = 0.1
97
98  CHARACTER (LEN=20), PARAMETER :: modname = 'cv3p1_closure'
99  CHARACTER (LEN=80) :: abort_message
100
101  ! print *,' -> cv3p1_closure, Ale ',ale(1)
102
103
104  ! -------------------------------------------------------
105  ! -- Initialization
106  ! -------------------------------------------------------
107
108
109  DO il = 1, ncum
110    alp2(il) = max(alp(il), 1.E-5)
111    ! IM
112    alp2(il) = max(alp(il), 1.E-12)
113  END DO
114
115  pbmxup = 50. ! PBMXUP+PBCRIT = cloud depth above which mixed updraughts
116  ! exist (if any)
117
118  IF (prt_level>=20) PRINT *, 'cv3p1_param nloc ncum nd icb inb nl', nloc, &
119    ncum, nd, icb(nloc), inb(nloc), nl
120  DO k = 1, nd          !jyg: initialization up to nd
121    DO il = 1, ncum
122      m(il, k) = 0.0
123    END DO
124  END DO
125
126!CR: initializations for erosion of adiabatic ascent
127  DO k = 1,nd           !jyg: initialization up to nd
128    DO il = 1, ncum
129        mad(il,k)=0.
130        me(il,k)=0.
131        betalim(il,k)=1.
132        wlim(il,k)=0.
133    ENDDO
134  ENDDO
135
136  ! -------------------------------------------------------
137  ! -- Reset sig(i) and w0(i) for i>inb and i<icb
138  ! -------------------------------------------------------
139
140  ! update sig and w0 above LNB:
141
142  DO k = 1, nl - 1
143    DO il = 1, ncum
144      IF ((inb(il)<(nl-1)) .AND. (k>=(inb(il)+1))) THEN
145        sig(il, k) = beta*sig(il, k) + 2.*alpha*buoy(il, inb(il))*abs(buoy(il &
146          ,inb(il)))
147        sig(il, k) = amax1(sig(il,k), 0.0)
148        w0(il, k) = beta*w0(il, k)
149      END IF
150    END DO
151  END DO
152
153  ! if(prt.level.GE.20) print*,'cv3p1_param apres 100'
154  ! compute icbmax:
155
156!ym break the column independance
157!ym  icbmax = 2
158!ym  DO il = 1, ncum
159!ym    icbmax = max(icbmax, icb(il))
160!ym  END DO
161
162  ! if(prt.level.GE.20) print*,'cv3p1_param apres 200'
163
164  ! update sig and w0 below cloud base:
165!ym column independance
166!ym  DO k = 1, icbmax
167  DO k = 1, nd
168    DO il = 1, ncum
169      IF (k<=MAX(2,icb(il))) THEN
170        IF (k<=icb(il)) THEN
171          sig(il, k) = beta*sig(il, k) - 2.*alpha*buoy(il, icb(il))*buoy(il, &
172            icb(il))
173          sig(il, k) = amax1(sig(il,k), 0.0)
174          w0(il, k) = beta*w0(il, k)
175        END IF
176      ENDIF
177    END DO
178  END DO
179  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 300'
180  ! -------------------------------------------------------------
181  ! -- Reset fractional areas of updrafts and w0 at initial time
182  ! -- and after 10 time steps of no convection
183  ! -------------------------------------------------------------
184
185  DO k = 1, nl - 1
186    DO il = 1, ncum
187      IF (sig(il,nd)<1.5 .OR. sig(il,nd)>12.0) THEN
188        sig(il, k) = 0.0
189        w0(il, k) = 0.0
190      END IF
191    END DO
192  END DO
193  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 400'
194
195  ! -------------------------------------------------------------
196  ! jyg1
197  ! --  Calculate adiabatic ascent top pressure (ptop)
198  ! -------------------------------------------------------------
199
200
201  ! c 1. Start at first level where precipitations form
202  DO il = 1, ncum
203    pzero(il) = plcl(il) - pbcrit
204  END DO
205
206  ! c 2. Add offset
207  DO il = 1, ncum
208    pzero(il) = pzero(il) - pbmxup
209  END DO
210  DO il = 1, ncum
211    ptop2old(il) = ptop2(il)
212  END DO
213
214  DO il = 1, ncum
215    ! CR:c est quoi ce 300??
216    p1(il) = pzero(il) - 300.
217  END DO
218
219  ! compute asupmax=abs(supmax) up to lnm+1
220
221  DO il = 1, ncum
222    ok(il) = .TRUE.
223    nsupmax(il) = inb(il)
224  END DO
225
226  DO i = 1, nl
227    DO il = 1, ncum
228      IF (i>icb(il) .AND. i<=inb(il)) THEN
229        IF (p(il,i)<=pzero(il) .AND. supmax(il,i)<0 .AND. ok(il)) THEN
230          nsupmax(il) = i
231          ok(il) = .FALSE.
232        END IF ! end IF (P(i) ...  )
233      END IF ! end IF (icb+1 le i le inb)
234    END DO
235  END DO
236
237  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 2.'
238  DO i = 1, nl
239    DO il = 1, ncum
240      asupmax(il, i) = abs(supmax(il,i))
241    END DO
242  END DO
243
244
245  DO il = 1, ncum
246    asupmaxmin(il) = 10.
247    pmin(il) = 100.
248    ! IM ??
249    asupmax0(il) = 0.
250  END DO
251
252  ! c 3.  Compute in which level is Pzero
253
254  ! IM bug      i0 = 18
255  DO il = 1, ncum
256    i0(il) = nl
257  END DO
258
259  DO i = 1, nl
260    DO il = 1, ncum
261      IF (i>icb(il) .AND. i<=inb(il)) THEN
262        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
263          IF (pzero(il)>p(il,i) .AND. pzero(il)<p(il,i-1)) THEN
264            i0(il) = i
265          END IF
266        END IF
267      END IF
268    END DO
269  END DO
270  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 3.'
271
272  ! c 4.  Compute asupmax at Pzero
273
274  DO i = 1, nl
275    DO il = 1, ncum
276      IF (i>icb(il) .AND. i<=inb(il)) THEN
277        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
278          asupmax0(il) = ((pzero(il)-p(il,i0(il)-1))*asupmax(il,i0(il))-( &
279            pzero(il)-p(il,i0(il)))*asupmax(il,i0(il)-1))/(p(il,i0(il))-p(il, &
280            i0(il)-1))
281        END IF
282      END IF
283    END DO
284  END DO
285
286
287  DO i = 1, nl
288    DO il = 1, ncum
289      IF (p(il,i)==pzero(il)) THEN
290        asupmax(i, il) = asupmax0(il)
291      END IF
292    END DO
293  END DO
294  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 4.'
295
296  ! c 5. Compute asupmaxmin, minimum of asupmax
297
298  DO i = 1, nl
299    DO il = 1, ncum
300      IF (i>icb(il) .AND. i<=inb(il)) THEN
301        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
302          IF (asupmax(il,i)<asupmaxmin(il)) THEN
303            asupmaxmin(il) = asupmax(il, i)
304            pmin(il) = p(il, i)
305          END IF
306        END IF
307      END IF
308    END DO
309  END DO
310
311  DO il = 1, ncum
312    ! IM
313    IF (prt_level>=20) THEN
314      PRINT *, 'cv3p1_closure il asupmax0 asupmaxmin', il, asupmax0(il), &
315        asupmaxmin(il), pzero(il), pmin(il)
316    END IF
317    IF (asupmax0(il)<asupmaxmin(il)) THEN
318      asupmaxmin(il) = asupmax0(il)
319      pmin(il) = pzero(il)
320    END IF
321  END DO
322  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 5.'
323
324
325  ! Compute Supmax at Pzero
326
327  DO i = 1, nl
328    DO il = 1, ncum
329      IF (i>icb(il) .AND. i<=inb(il)) THEN
330        IF (p(il,i)<=pzero(il)) THEN
331          supmax0(il) = ((p(il,i)-pzero(il))*asupmax(il,i-1)-(p(il, &
332            i-1)-pzero(il))*asupmax(il,i))/(p(il,i)-p(il,i-1))
333!ym WARNING : probably bad GOTO branching ===> to check !
334!ym          GO TO 425
335        END IF ! end IF (P(i) ... )
336      END IF ! end IF (icb+1 le i le inb)
337    END DO
338  END DO
339!ym bad branching
340!ym 425 CONTINUE
341  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 425.'
342
343  ! c 6. Calculate ptop2
344
345  DO il = 1, ncum
346    IF (asupmaxmin(il)<supcrit1) THEN
347      ptop2(il) = pmin(il)
348    END IF
349
350    IF (asupmaxmin(il)>supcrit1 .AND. asupmaxmin(il)<supcrit2) THEN
351      ptop2(il) = ptop2old(il)
352    END IF
353
354    IF (asupmaxmin(il)>supcrit2) THEN
355      ptop2(il) = ph(il, inb(il))
356    END IF
357  END DO
358
359  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 6.'
360
361  ! c 7. Compute multiplying factor for adiabatic updraught mass flux
362
363
364  IF (ok_inhib) THEN
365
366    DO i = 1, nl
367      DO il = 1, ncum
368        IF (i<=nl) THEN
369          coefmix(il, i) = (min(ptop2(il),ph(il,i))-ph(il,i))/(ph(il,i+1)-ph( &
370            il,i))
371          coefmix(il, i) = min(coefmix(il,i), 1.)
372        END IF
373      END DO
374    END DO
375
376
377  ELSE ! when inhibition is not taken into account, coefmix=1
378
379
380
381    DO i = 1, nl
382      DO il = 1, ncum
383        IF (i<=nl) THEN
384          coefmix(il, i) = 1.
385        END IF
386      END DO
387    END DO
388
389  END IF ! ok_inhib
390  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 7.'
391  ! -------------------------------------------------------------------
392  ! -------------------------------------------------------------------
393
394
395  ! jyg2
396
397  ! ==========================================================================
398
399
400  ! -------------------------------------------------------------
401  ! -- Calculate convective inhibition (CIN)
402  ! -------------------------------------------------------------
403
404  ! do i=1,nloc
405  ! print*,'avant cine p',pbase(i),plcl(i)
406  ! enddo
407  ! do j=1,nd
408  ! do i=1,nloc
409  ! print*,'avant cine t',tv(i),tvp(i)
410  ! enddo
411  ! enddo
412  CALL cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, &
413    cinb, plfc)
414
415  DO il = 1, ncum
416    cin(il) = cina(il) + cinb(il)
417  END DO
418  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_cine'
419  ! -------------------------------------------------------------
420  ! --Update buoyancies to account for Ale
421  ! -------------------------------------------------------------
422
423  CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, &
424    tvp, buoy)
425  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_buoy'
426
427  ! -------------------------------------------------------------
428  ! -- Calculate convective available potential energy (cape),
429  ! -- vertical velocity (w), fractional area covered by
430  ! -- undilute updraft (sig), and updraft mass flux (m)
431  ! -------------------------------------------------------------
432
433  DO il = 1, ncum
434    cape(il) = 0.0
435  END DO
436
437  ! compute dtmin (minimum buoyancy between ICB and given level k):
438
439  DO k = 1, nl
440    DO il = 1, ncum
441      dtmin(il, k) = 100.0
442    END DO
443  END DO
444
445  DO k = 1, nl
446    DO j = minorig, nl
447      DO il = 1, ncum
448        IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (j>=icb(il)) .AND. (j<= &
449            (k-1))) THEN
450          dtmin(il, k) = amin1(dtmin(il,k), buoy(il,j))
451        END IF
452      END DO
453    END DO
454  END DO
455
456  ! the interval on which cape is computed starts at pbase :
457
458  DO k = 1, nl
459    DO il = 1, ncum
460
461      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
462        IF (iflag_mix_adiab.eq.1) THEN
463!CR:computation of cape from LCL: keep flag or to modify in all cases?
464        deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k))
465        ELSE
466        deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
467        ENDIF
468        cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1)
469        cape(il) = amax1(0.0, cape(il))
470        sigold(il, k) = sig(il, k)
471
472
473        ! jyg       Coefficient coefmix limits convection to levels where a
474        ! sufficient
475        ! fraction of mixed draughts are ascending.
476        siglim(il, k) = coefmix(il, k)*alpha1*dtmin(il, k)*abs(dtmin(il,k))
477        siglim(il, k) = amax1(siglim(il,k), 0.0)
478        siglim(il, k) = amin1(siglim(il,k), 0.01)
479        ! c         fac=AMIN1(((dtcrit-dtmin(il,k))/dtcrit),1.0)
480        fac = 1.
481        wlim(il, k) = fac*sqrt(cape(il))
482        amu = siglim(il, k)*wlim(il, k)
483        rhodp = 0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
484        mlim(il, k) = amu*rhodp
485        ! print*, 'siglim ', k,siglim(1,k)
486      END IF
487
488    END DO
489  END DO
490  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 600'
491
492  DO il = 1, ncum
493    ! IM beg
494    IF (prt_level>=20) THEN
495      PRINT *, 'cv3p1_closure il icb mlim ph ph+1 ph+2', il, icb(il), &
496        mlim(il, icb(il)+1), ph(il, icb(il)), ph(il, icb(il)+1), &
497        ph(il, icb(il)+2)
498    END IF
499
500    IF (icb(il)+1<=inb(il)) THEN
501      ! IM end
502      mlim(il, icb(il)) = 0.5*mlim(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb( &
503        il)+1))/(ph(il,icb(il)+1)-ph(il,icb(il)+2))
504      ! IM beg
505    END IF !(icb(il.le.inb(il))) then
506    ! IM end
507  END DO
508  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 700'
509
510  ! jyg1
511  ! ------------------------------------------------------------------------
512  ! c     Correct mass fluxes so that power used to overcome CIN does not
513  ! c     exceed Power Available for Lifting (PAL).
514  ! ------------------------------------------------------------------------
515
516  DO il = 1, ncum
517    cbmflim(il) = 0.
518    cbmf(il) = 0.
519  END DO
520
521  ! c 1. Compute cloud base mass flux of elementary system (Cbmf0=Cbmflim)
522
523  DO k = 1, nl
524    DO il = 1, ncum
525      ! old       IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
526      ! IM        IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
527      IF (k>=icb(il) .AND. k<=inb(il) & !cor jyg
528          .AND. icb(il)+1<=inb(il)) THEN !cor jyg
529        cbmflim(il) = cbmflim(il) + mlim(il, k)
530      END IF
531    END DO
532  END DO
533  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim'
534
535  ! 1.5 Compute cloud base mass flux given by Alp closure (Cbmf1), maximum
536  !     allowed mass flux (Cbmfmax) and final target mass flux (Cbmf)
537  !     Cbmf is set to zero if Cbmflim (the mass flux of elementary cloud)
538  !     is exceedingly small.
539
540  DO il = 1, ncum
541    wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.))
542  END DO
543
544  DO il = 1, ncum
545    IF (plfc(il)<100.) THEN
546      ! This is an irealistic value for plfc => no calculation of wbeff
547      wbeff(il) = 100.1
548    ELSE
549      ! Calculate wbeff
550      IF (NINT(flag_wb)==0) THEN
551        wbeff(il) = wbmax
552      ELSE IF (NINT(flag_wb)==1) THEN
553        wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
554      ELSE IF (NINT(flag_wb)==2) THEN
555        wbeff(il) = wbmax*(0.01*(ph(il,1)-plfc(il)))**2
556      ELSE ! Option provisoire ou le iflag_wb/10 est considere comme une vitesse
557        wbeff(il) = flag_wb*0.01+wbmax/(1.+500./(ph(il,1)-plfc(il)))
558      END IF
559    END IF
560  END DO
561
562!CR:Compute k at plfc
563  DO il=1,ncum
564           klfc(il)=nl
565  ENDDO
566  DO k=1,nl
567     DO il=1,ncum
568        if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
569           klfc(il)=k
570        endif
571     ENDDO
572  ENDDO
573!RC
574
575  DO il = 1, ncum
576    ! jyg    Modification du coef de wb*wb pour conformite avec papier Wake
577    ! c       cbmf1(il) = alp2(il)/(0.5*wb*wb-Cin(il))
578    cbmf1(il) = alp2(il)/(2.*wbeff(il)*wbeff(il)-cin(il))
579!CR: Add large-scale component to the mass-flux
580!encore connu sous le nom "Experience du tube de dentifrice"
581    if ((coef_clos_ls.gt.0.).and.(plfc(il).gt.0.)) then
582       cbmf1(il) = cbmf1(il) - coef_clos_ls*min(0.,1./RG*omega(il,klfc(il)))
583    endif
584!RC
585    IF (cbmf1(il)==0 .AND. alp2(il)/=0.) THEN
586      WRITE (lunout, *) 'cv3p1_closure cbmf1=0 and alp NE 0 il alp2 alp cin ' &
587        , il, alp2(il), alp(il), cin(il)
588      abort_message = ''
589      CALL abort_physic(modname, abort_message, 1)
590    END IF
591    cbmfmax(il) = sigmax*wb2(il)*100.*p(il, icb(il))/(rrd*tv(il,icb(il)))
592  END DO
593
594  DO il = 1, ncum
595    IF (cbmflim(il)>1.E-6) THEN
596      ! ATTENTION TEST CR
597      ! if (cbmfmax(il).lt.1.e-12) then
598      cbmf(il) = min(cbmf1(il), cbmfmax(il))
599      ! else
600      ! cbmf(il) = cbmf1(il)
601      ! endif
602      ! print*,'cbmf',cbmf1(il),cbmfmax(il)
603    END IF
604  END DO
605  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim_testCR'
606
607  ! c 2. Compute coefficient and apply correction
608
609  DO il = 1, ncum
610    coef(il) = (cbmf(il)+1.E-10)/(cbmflim(il)+1.E-10)
611    coeftrue(il) = coef(il)
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
785END MODULE cv3p1_closure_mod
786
Note: See TracBrowser for help on using the repository browser.