
! $Id: cv3p1_closure_mod.f90 6058 2026-01-28 23:52:39Z fhourdin $
MODULE cv3p1_closure_mod
PRIVATE
PUBLIC cv3p1_closure
CONTAINS

SUBROUTINE cv3p1_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
    tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, &
    iflag, coef, coeftrue, plim1, plim2, asupmax, supmax0, asupmaxmin, &
    cbmf, plfc, wbeff)


  ! **************************************************************
  ! *
  ! CV3P1_CLOSURE                                               *
  ! Ale & Alp Closure of Convect3              *
  ! *
  ! written by   :   Kerry Emanuel                              *
  ! vectorization:   S. Bony                                    *
  ! modified by :    Jean-Yves Grandpeix, 18/06/2003, 19.32.10  *
  ! Julie Frohwirth,     14/10/2005  17.44.22  *
  ! **************************************************************

! Replayisation USE yomcst2_mod_h
  USE lmdz_cv_ini, ONLY : iflag_mix_adiab,coef_clos_ls,supcrit1,supcrit2
   USE lmdz_cv_ini, ONLY : pbcrit,wbmax,rrd,minorig,flag_wb,alpha,alpha1,beta,coef_peel,nl
! Replayisation USE conema3_mod_h
  USE lmdz_cv_ini, ONLY : iflag_cvl_sigd
   USE lmdz_cv_ini, ONLY: prt_level, lunout
! Replayisation USE yomcst_mod_h
  USE lmdz_cv_ini, ONLY : g
  USE cv3_cine_mod, ONLY : cv3_cine
  USE cv3_buoy_mod, ONLY : cv3_buoy
IMPLICIT NONE



  ! input:
  INTEGER, INTENT (IN)                               :: ncum, nd, nloc
  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb
  REAL, DIMENSION (nloc), INTENT (IN)                :: pbase, plcl
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: p
  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: tv, tvp, buoy
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: supmax
  LOGICAL, INTENT (IN)                               :: ok_inhib ! enable convection inhibition by dryness
  REAL, DIMENSION (nloc), INTENT (IN)                :: ale, alp
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: omega

  ! input/output:
  INTEGER, DIMENSION (nloc), INTENT (INOUT)          :: iflag
  REAL, DIMENSION (nloc, nd), INTENT (INOUT)         :: sig, w0
  REAL, DIMENSION (nloc), INTENT (INOUT)             :: ptop2

  ! output:
  REAL, DIMENSION (nloc), INTENT (OUT)               :: cape, cin
  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: m
  REAL, DIMENSION (nloc), INTENT (OUT)               :: coef, coeftrue
  REAL, DIMENSION (nloc), INTENT (OUT)               :: plim1, plim2
  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: asupmax
  REAL, DIMENSION (nloc), INTENT (OUT)               :: supmax0
  REAL, DIMENSION (nloc), INTENT (OUT)               :: asupmaxmin
  REAL, DIMENSION (nloc), INTENT (OUT)               :: cbmf, plfc
  REAL, DIMENSION (nloc), INTENT (OUT)               :: wbeff

  ! local variables:
  INTEGER il, i, j, k, icbmax, i0(nloc), klfc(nloc)
  REAL deltap, fac, w, amu
  REAL rhodp, dz
  REAL pbmxup
  REAL dtmin(nloc, nd), sigold(nloc, nd)
  REAL coefmix(nloc, nd)
  REAL pzero(nloc), ptop2old(nloc)
  REAL cina(nloc), cinb(nloc)
  INTEGER ibeg(nloc)
  INTEGER nsupmax(nloc)
  REAL supcrit, temp(nloc, nd)
  REAL p1(nloc), pmin(nloc)
  REAL asupmax0(nloc)
  LOGICAL ok(nloc)
  REAL siglim(nloc, nd), wlim(nloc, nd), mlim(nloc, nd)
  REAL wb2(nloc)
  REAL cbmflim(nloc), cbmf1(nloc), cbmfmax(nloc)
  REAL cbmflast(nloc)
  REAL xp(nloc), xq(nloc), xr(nloc), discr(nloc), b3(nloc), b4(nloc)
  REAL theta(nloc), bb(nloc)
  REAL term1, term2, term3
  REAL alp2(nloc) ! Alp with offset
  !CR: variables for new erosion of adiabiatic ascent
  REAL mad(nloc, nd), me(nloc, nd), betalim(nloc, nd), beta_coef(nloc, nd)
  REAL med(nloc, nd), md(nloc,nd)
!jyg<
! coef_peel is now in the common cv3_param 
!!  REAL coef_peel
!!  PARAMETER (coef_peel=0.25)
!>jyg

  REAL,PARAMETER :: sigmax = 0.1

  CHARACTER (LEN=20), PARAMETER :: modname = 'cv3p1_closure'
  CHARACTER (LEN=80) :: abort_message

  ! print *,' -> cv3p1_closure, Ale ',ale(1)

  REAL :: rg

  rg=g ! FH 2026/01/27  Historique convergence constantes cv et yomcst a nettoyer

  ! -------------------------------------------------------
  ! -- Initialization
  ! -------------------------------------------------------


  DO il = 1, ncum
    alp2(il) = max(alp(il), 1.E-5)
    ! IM
    alp2(il) = max(alp(il), 1.E-12)
  END DO

  pbmxup = 50. ! PBMXUP+PBCRIT = cloud depth above which mixed updraughts
  ! exist (if any)

  IF (prt_level>=20) PRINT *, 'cv3p1_param nloc ncum nd icb inb nl', nloc, &
    ncum, nd, icb(nloc), inb(nloc), nl
  DO k = 1, nd          !jyg: initialization up to nd
    DO il = 1, ncum
      m(il, k) = 0.0
    END DO
  END DO

!CR: initializations for erosion of adiabatic ascent 
  DO k = 1,nd           !jyg: initialization up to nd
    DO il = 1, ncum
        mad(il,k)=0.
        me(il,k)=0.
        betalim(il,k)=1.
        wlim(il,k)=0.
    ENDDO
  ENDDO

  ! -------------------------------------------------------
  ! -- Reset sig(i) and w0(i) for i>inb and i<icb
  ! -------------------------------------------------------

  ! update sig and w0 above LNB:

  DO k = 1, nl - 1
    DO il = 1, ncum
      IF ((inb(il)<(nl-1)) .AND. (k>=(inb(il)+1))) THEN
        sig(il, k) = beta*sig(il, k) + 2.*alpha*buoy(il, inb(il))*abs(buoy(il &
          ,inb(il)))
        sig(il, k) = amax1(sig(il,k), 0.0)
        w0(il, k) = beta*w0(il, k)
      END IF
    END DO
  END DO

  ! if(prt.level.GE.20) print*,'cv3p1_param apres 100'
  ! compute icbmax:

!ym break the column independance
!ym  icbmax = 2
!ym  DO il = 1, ncum
!ym    icbmax = max(icbmax, icb(il))
!ym  END DO

  ! if(prt.level.GE.20) print*,'cv3p1_param apres 200'

  ! update sig and w0 below cloud base:
!ym column independance
!ym  DO k = 1, icbmax

! JYCF2026/01/20 m(k)=rho(k)*sig(k)*w(k)
! sig(k) est la fraction surfacique de l'ascendance arrivant au niveau k (en m^2/m^2)

  DO k = 1, nd
    DO il = 1, ncum
      IF (k<=MAX(2,icb(il))) THEN
        IF (k<=icb(il)) THEN
          sig(il, k) = beta*sig(il, k) - 2.*alpha*buoy(il, icb(il))*buoy(il, &
            icb(il))
          sig(il, k) = amax1(sig(il,k), 0.0)
          w0(il, k) = beta*w0(il, k)
        END IF
      ENDIF
    END DO
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 300'
  ! -------------------------------------------------------------
  ! -- Reset fractional areas of updrafts and w0 at initial time
  ! -- and after 10 time steps of no convection
  ! -------------------------------------------------------------

  DO k = 1, nl - 1
    DO il = 1, ncum
      IF (sig(il,nd)<1.5 .OR. sig(il,nd)>12.0) THEN
        sig(il, k) = 0.0
        w0(il, k) = 0.0
      END IF
    END DO
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 400'

  ! -------------------------------------------------------------
  ! jyg1
  ! --  Calculate adiabatic ascent top pressure (ptop)
  ! -------------------------------------------------------------


  ! c 1. Start at first level where precipitations form
  DO il = 1, ncum
    pzero(il) = plcl(il) - pbcrit
  END DO

  ! c 2. Add offset
  DO il = 1, ncum
    pzero(il) = pzero(il) - pbmxup
  END DO
  DO il = 1, ncum
    ptop2old(il) = ptop2(il)
  END DO

  DO il = 1, ncum
    ! CR:c est quoi ce 300??
    p1(il) = pzero(il) - 300.
  END DO

  ! compute asupmax=abs(supmax) up to lnm+1

  DO il = 1, ncum
    ok(il) = .TRUE.
    nsupmax(il) = inb(il)
  END DO

  DO i = 1, nl
    DO il = 1, ncum
      IF (i>icb(il) .AND. i<=inb(il)) THEN
        IF (p(il,i)<=pzero(il) .AND. supmax(il,i)<0 .AND. ok(il)) THEN
          nsupmax(il) = i
          ok(il) = .FALSE.
        END IF ! end IF (P(i) ...  )
      END IF ! end IF (icb+1 le i le inb)
    END DO
  END DO

  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 2.'
  DO i = 1, nl
    DO il = 1, ncum
      asupmax(il, i) = abs(supmax(il,i))
    END DO
  END DO


  DO il = 1, ncum
    asupmaxmin(il) = 10.
    pmin(il) = 100.
    ! IM ??
    asupmax0(il) = 0.
  END DO

  ! c 3.  Compute in which level is Pzero

  ! IM bug      i0 = 18
  DO il = 1, ncum
    i0(il) = nl
  END DO

  DO i = 1, nl
    DO il = 1, ncum
      IF (i>icb(il) .AND. i<=inb(il)) THEN
        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
          IF (pzero(il)>p(il,i) .AND. pzero(il)<p(il,i-1)) THEN
            i0(il) = i
          END IF
        END IF
      END IF
    END DO
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 3.'

  ! c 4.  Compute asupmax at Pzero

  DO i = 1, nl
    DO il = 1, ncum
      IF (i>icb(il) .AND. i<=inb(il)) THEN
        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
          asupmax0(il) = ((pzero(il)-p(il,i0(il)-1))*asupmax(il,i0(il))-( &
            pzero(il)-p(il,i0(il)))*asupmax(il,i0(il)-1))/(p(il,i0(il))-p(il, &
            i0(il)-1))
        END IF
      END IF
    END DO
  END DO


  DO i = 1, nl
    DO il = 1, ncum
      IF (p(il,i)==pzero(il)) THEN
        asupmax(i, il) = asupmax0(il)
      END IF
    END DO
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 4.'

  ! c 5. Compute asupmaxmin, minimum of asupmax

  DO i = 1, nl
    DO il = 1, ncum
      IF (i>icb(il) .AND. i<=inb(il)) THEN
        IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
          IF (asupmax(il,i)<asupmaxmin(il)) THEN
            asupmaxmin(il) = asupmax(il, i)
            pmin(il) = p(il, i)
          END IF
        END IF
      END IF
    END DO
  END DO

  DO il = 1, ncum
    ! IM
    IF (prt_level>=20) THEN
      PRINT *, 'cv3p1_closure il asupmax0 asupmaxmin', il, asupmax0(il), &
        asupmaxmin(il), pzero(il), pmin(il)
    END IF
    IF (asupmax0(il)<asupmaxmin(il)) THEN
      asupmaxmin(il) = asupmax0(il)
      pmin(il) = pzero(il)
    END IF
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 5.'


  ! Compute Supmax at Pzero

  DO i = 1, nl
    DO il = 1, ncum
      IF (i>icb(il) .AND. i<=inb(il)) THEN
        IF (p(il,i)<=pzero(il)) THEN
          supmax0(il) = ((p(il,i)-pzero(il))*asupmax(il,i-1)-(p(il, &
            i-1)-pzero(il))*asupmax(il,i))/(p(il,i)-p(il,i-1))
!ym WARNING : probably bad GOTO branching ===> to check !
!ym          GO TO 425
        END IF ! end IF (P(i) ... )
      END IF ! end IF (icb+1 le i le inb)
    END DO
  END DO
!ym bad branching
!ym 425 CONTINUE
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 425.'

  ! c 6. Calculate ptop2

  DO il = 1, ncum
    IF (asupmaxmin(il)<supcrit1) THEN
      ptop2(il) = pmin(il)
    END IF

    IF (asupmaxmin(il)>supcrit1 .AND. asupmaxmin(il)<supcrit2) THEN
      ptop2(il) = ptop2old(il)
    END IF

    IF (asupmaxmin(il)>supcrit2) THEN
      ptop2(il) = ph(il, inb(il))
    END IF
  END DO

  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 6.'

  ! c 7. Compute multiplying factor for adiabatic updraught mass flux


  IF (ok_inhib) THEN

! JYCF2026/01/20 m(k)=rho(k)*sig(k)*w(k)
! Possible que toutes les lignes au dessus ne  servent que pour ok_inhib=T
! Ce n'est pas le cas.
! ok_inhib est mis à iflag_mix==2 dans cva_driver.
! Il est donc faux par défaut, pour iflag_mix=1

    DO i = 1, nl
      DO il = 1, ncum
        IF (i<=nl) THEN
          coefmix(il, i) = (min(ptop2(il),ph(il,i))-ph(il,i))/(ph(il,i+1)-ph( &
            il,i))
          coefmix(il, i) = min(coefmix(il,i), 1.)
        END IF
      END DO
    END DO


  ELSE ! when inhibition is not taken into account, coefmix=1



    DO i = 1, nl
      DO il = 1, ncum
        IF (i<=nl) THEN
          coefmix(il, i) = 1.
        END IF
      END DO
    END DO

  END IF ! ok_inhib
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 7.'
  ! -------------------------------------------------------------------
  ! -------------------------------------------------------------------


  ! jyg2

  ! ==========================================================================


  ! -------------------------------------------------------------
  ! -- Calculate convective inhibition (CIN)
  ! -------------------------------------------------------------

  ! do i=1,nloc
  ! print*,'avant cine p',pbase(i),plcl(i)
  ! enddo
  ! do j=1,nd
  ! do i=1,nloc
  ! print*,'avant cine t',tv(i),tvp(i)
  ! enddo
  ! enddo
  CALL cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, &
    cinb, plfc)

! JYCF2026/01/20 somme de la partie en dessous et au dessus du LCL

  DO il = 1, ncum
    cin(il) = cina(il) + cinb(il)
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_cine'

  ! -------------------------------------------------------------
  ! --Update buoyancies to account for Ale
  ! -------------------------------------------------------------

! JYCF2026/01/20 Ajout de Ale dans le caclul de ALE pour le déclenchement
!  Est-ce que le calcul du déclenchement est fait dedans.
!  A verifier et documenter.
!  Notamment mettre les inten(in/out) dans cv3_buoy.f90
!  On ne sort buoy que si ALE+CIN>0 ?
!  Le "buoy" qui sort est un delta de température ?

  CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, &
    tvp, buoy)
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_buoy'

  ! -------------------------------------------------------------
  ! -- Calculate convective available potential energy (cape),
  ! -- vertical velocity (w), fractional area covered by
  ! -- undilute updraft (sig), and updraft mass flux (m)
  ! -------------------------------------------------------------

  DO il = 1, ncum
    cape(il) = 0.0
  END DO

  ! compute dtmin (minimum buoyancy between ICB and given level k):
  ! in fact a delta of temperature
  DO k = 1, nl
    DO il = 1, ncum
      dtmin(il, k) = 100.0
    END DO
  END DO

  DO k = 1, nl
    DO j = minorig, nl
      DO il = 1, ncum
        IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (j>=icb(il)) .AND. (j<= &
            (k-1))) THEN
          dtmin(il, k) = amin1(dtmin(il,k), buoy(il,j))
        END IF
      END DO
    END DO
  END DO

  ! the vertical interval on which cape is computed starts at pbase :

! JYCF2026/01/20 computing siglim(k), the target value for sig(k)
!   for the time relaxation
!   mlim=rho*dp*siglim*wlim avec wlim=sqrt(cape)


  DO k = 1, nl
    DO il = 1, ncum

      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
        IF (iflag_mix_adiab.eq.1) THEN
!CR:computation of cape from LCL: keep flag or to modify in all cases?
        deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k))
        ELSE
        deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
        ENDIF
        cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1)
        cape(il) = amax1(0.0, cape(il))
        sigold(il, k) = sig(il, k)


        ! jyg       Coefficient coefmix limits convection to levels where a
        ! sufficient
        ! fraction of mixed draughts are ascending.
        siglim(il, k) = coefmix(il, k)*alpha1*dtmin(il, k)*abs(dtmin(il,k))
        siglim(il, k) = amax1(siglim(il,k), 0.0)
        siglim(il, k) = amin1(siglim(il,k), 0.01)
        ! c         fac=AMIN1(((dtcrit-dtmin(il,k))/dtcrit),1.0)
        fac = 1.
        wlim(il, k) = fac*sqrt(cape(il))
        amu = siglim(il, k)*wlim(il, k)
        rhodp = 0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
        mlim(il, k) = amu*rhodp
        ! print*, 'siglim ', k,siglim(1,k)
      END IF

    END DO
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 600'

  DO il = 1, ncum
    ! IM beg
    IF (prt_level>=20) THEN
      PRINT *, 'cv3p1_closure il icb mlim ph ph+1 ph+2', il, icb(il), &
        mlim(il, icb(il)+1), ph(il, icb(il)), ph(il, icb(il)+1), &
        ph(il, icb(il)+2)
    END IF

! JYCF2026/01/20 computing mlim at CB
! Boucle à vérifier
! inb est le sommet

    IF (icb(il)+1<=inb(il)) THEN
      ! IM end
      mlim(il, icb(il)) = 0.5*mlim(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb( &
        il)+1))/(ph(il,icb(il)+1)-ph(il,icb(il)+2))
      ! IM beg
    END IF !(icb(il.le.inb(il))) then
    ! IM end
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres 700'

  ! jyg1
  ! ------------------------------------------------------------------------
  ! c     Correct mass fluxes so that power used to overcome CIN does not
  ! c     exceed Power Available for Lifting (PAL).
  ! ------------------------------------------------------------------------


! JYCF2026/01/20 cbmflim=sum(mlim)
  DO il = 1, ncum
    cbmflim(il) = 0.
    cbmf(il) = 0.
  END DO

  ! c 1. Compute cloud base mass flux of elementary system (Cbmf0=Cbmflim)

  DO k = 1, nl
    DO il = 1, ncum
      ! old       IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
      ! IM        IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
      IF (k>=icb(il) .AND. k<=inb(il) & !cor jyg
          .AND. icb(il)+1<=inb(il)) THEN !cor jyg
        cbmflim(il) = cbmflim(il) + mlim(il, k)
      END IF
    END DO
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim'

  ! 1.5 Compute cloud base mass flux given by Alp closure (Cbmf1), maximum
  !     allowed mass flux (Cbmfmax) and final target mass flux (Cbmf)
  !     Cbmf is set to zero if Cbmflim (the mass flux of elementary cloud)
  !     is exceedingly small.

  DO il = 1, ncum
    wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.))
  END DO

! JYCF2026/01/20 Various options to compute the vertical velocity
! in the adiabatic ascent at cloud base

  DO il = 1, ncum
    IF (plfc(il)<100.) THEN
      ! This is an irealistic value for plfc => no calculation of wbeff
      wbeff(il) = 100.1
    ELSE
      ! Calculate wbeff
      IF (NINT(flag_wb)==0) THEN
        wbeff(il) = wbmax
      ELSE IF (NINT(flag_wb)==1) THEN
        wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
      ELSE IF (NINT(flag_wb)==2) THEN
        wbeff(il) = wbmax*(0.01*(ph(il,1)-plfc(il)))**2
      ELSE ! Option provisoire ou le iflag_wb/10 est considere comme une vitesse
        wbeff(il) = flag_wb*0.01+wbmax/(1.+500./(ph(il,1)-plfc(il)))
      END IF
    END IF
  END DO

!CR:Compute k at plfc
  DO il=1,ncum
           klfc(il)=nl
  ENDDO
  DO k=1,nl
     DO il=1,ncum
        if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
           klfc(il)=k
        endif
     ENDDO
  ENDDO
!RC

! JYCF2026/01/20 computing cbmf1, mass flux at cloud base coming
! from the ALP closure

  DO il = 1, ncum
    ! jyg    Modification du coef de wb*wb pour conformite avec papier Wake
    ! c       cbmf1(il) = alp2(il)/(0.5*wb*wb-Cin(il))
    cbmf1(il) = alp2(il)/(2.*wbeff(il)*wbeff(il)-cin(il))
!CR: Add large-scale component to the mass-flux
!encore connu sous le nom "Experience du tube de dentifrice"
    if ((coef_clos_ls.gt.0.).and.(plfc(il).gt.0.)) then 
       cbmf1(il) = cbmf1(il) - coef_clos_ls*min(0.,1./RG*omega(il,klfc(il)))
    endif
!RC
    IF (cbmf1(il)==0 .AND. alp2(il)/=0.) THEN
      WRITE (lunout, *) 'cv3p1_closure cbmf1=0 and alp NE 0 il alp2 alp cin ' &
        , il, alp2(il), alp(il), cin(il)
      abort_message = ''
      CALL abort_physic(modname, abort_message, 1)
    END IF
    cbmfmax(il) = sigmax*wb2(il)*100.*p(il, icb(il))/(rrd*tv(il,icb(il)))
  END DO

! JYCF2026/01/20 cbmf=cbmf1, coming from ALP closure
  DO il = 1, ncum
    IF (cbmflim(il)>1.E-6) THEN
      ! ATTENTION TEST CR
      ! if (cbmfmax(il).lt.1.e-12) then
      cbmf(il) = min(cbmf1(il), cbmfmax(il))
      ! else
      ! cbmf(il) = cbmf1(il)
      ! endif
      ! print*,'cbmf',cbmf1(il),cbmfmax(il)
    END IF
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim_testCR'

  ! c 2. Compute coefficient and apply correction

! JYCF2026/01/20 coef : ratio of the ALP to the CAPE closure
! keeping the original ratio in coeftrue

  DO il = 1, ncum
    coef(il) = (cbmf(il)+1.E-10)/(cbmflim(il)+1.E-10)
    coeftrue(il) = coef(il)
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres coef_plantePLUS'

! JYCF2026/01/20 : resumé
! computing siglim(k), the target value for sig(k)
!   for the time relaxation
!   mlim=rho*dp*siglim*wlim avec wlim=sqrt(cape)
! cbmflim=sum(mlim)
! computing cbmf, mass flux at cloud base from the ALP closure
! time relaxation on sig*w (in fact on "m")
! w is kept unchanged and sig is computed as sig*w/w
!
! coef=cbmf(ALP)/cbmflim(CAPE) sans relaxation
! 


! JYCF2026/01/20 time relaxation on sig*w (in fact on "m")
! beta = 1.0 - delt / tau, in cv3_param
! w is kept unchanged and sig is computed as sig*w/w

  DO k = 1, nl
    DO il = 1, ncum
      IF (k>=icb(il)+1 .AND. k<=inb(il)) THEN
        amu = beta*sig(il, k)*w0(il, k) + (1.-beta)*coef(il)*siglim(il, k)* &
          wlim(il, k)
        w0(il, k) = wlim(il, k)
        w0(il, k) = max(w0(il,k), 1.E-10)
        sig(il, k) = amu/w0(il, k)
        sig(il, k) = min(sig(il,k), 1.)
        ! c         amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k)
! JYCF2026/01/20 0.007=2/R ?
!  m=sig w rho delta(p)
        m(il, k) = amu*0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
      END IF
    END DO
  END DO
  ! jyg2
  DO il = 1, ncum
    w0(il, icb(il)) = 0.5*w0(il, icb(il)+1)
    m(il, icb(il)) = 0.5*m(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb(il)+1))/ &
      (ph(il,icb(il)+1)-ph(il,icb(il)+2))
    sig(il, icb(il)) = sig(il, icb(il)+1)
    sig(il, icb(il)-1) = sig(il, icb(il))
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres w0_sig_M'

!CR: new erosion of adiabatic ascent: modification of m
!computation of the sum of ascending fluxes 
  IF (iflag_mix_adiab.eq.1) THEN
      !
      !Verification sum(me)=sum(m)
        DO k = 1,nd                         !jyg: initialization up to nd
          DO il = 1, ncum
             md(il,k)=0.
             med(il,k)=0.
          ENDDO
        ENDDO
      !
        DO k = nl,1,-1
          DO il = 1, ncum
                 md(il,k)=md(il,k+1)+m(il,k+1)
          ENDDO
        ENDDO
      !
        DO k = nl,1,-1
          DO il = 1, ncum
              IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
                 mad(il,k)=mad(il,k+1)+m(il,k+1)
              ENDIF
      !        print*,"mad",il,k,mad(il,k)
          ENDDO
        ENDDO
      !
      !CR: erosion of each adiabatic ascent during its ascent
      !
      !Computation of erosion coefficient beta_coef
        DO k = 1, nl
          DO il = 1, ncum
             IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN     
      !          print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
                beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
             ELSE
                beta_coef(il,k)=0.
             ENDIF
          ENDDO
        ENDDO
      !
      !  print*,"apres beta_coef"
      !
        DO k = 1, nl
          DO il = 1, ncum
      !
            IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
      !
      !        print*,"dz",il,k,tv(il, k-1)
              dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG)
              betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz)
      !        betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz)
      !        print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
              dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG)
      !        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))
              me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k))
      !        print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz    
      !
            END IF
      !
      !Modification of m
            m(il,k)=me(il,k) 
          END DO
        END DO
      !
      !  DO il = 1, ncum
      !     dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG)
      !     m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) &
      !                  /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il))
      !     print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
      !  ENDDO
      !
      !Verification sum(me)=sum(m)
        DO k = nl,1,-1
          DO il = 1, ncum
                 med(il,k)=med(il,k+1)+m(il,k+1)
      !           print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
          ENDDO
        ENDDO
      !
      !
  ENDIF !(iflag_mix_adiab)
!RC



  ! c 3. Compute final cloud base mass flux and set iflag to 3 if
  ! c    cloud base mass flux is exceedingly small and is decreasing (i.e. if
  ! c    the final mass flux (cbmflast) is greater than the target mass flux
  ! c    (cbmf)).

  DO il = 1, ncum
    cbmflast(il) = 0.
  END DO

  DO k = 1, nl
    DO il = 1, ncum
      IF (k>=icb(il) .AND. k<=inb(il)) THEN
          !IMpropo??      IF ((k.ge.(icb(il)+1)).and.(k.le.inb(il))) THEN
        cbmflast(il) = cbmflast(il) + m(il, k)
      END IF
    END DO
  END DO

  DO il = 1, ncum
    IF (cbmflast(il)<1.E-6 .AND. cbmflast(il)>=cbmf(il)) THEN
      iflag(il) = 3
    END IF
  END DO

  DO k = 1, nl
    DO il = 1, ncum
      IF (iflag(il)>=3) THEN
        m(il, k) = 0.
        sig(il, k) = 0.
        w0(il, k) = 0.
      END IF
    END DO
  END DO
  IF (prt_level>=20) PRINT *, 'cv3p1_param apres iflag'

  ! c 4. Introduce a correcting factor for coef, in order to obtain an
  ! effective sigdz larger in the present case (using cv3p1_closure)
  !  than in the old
  ! c    closure (using cv3_closure).
  IF (1==0) THEN
    DO il = 1, ncum
      ! c      coef(il) = 2.*coef(il)
      coef(il) = 5.*coef(il)
    END DO
    ! version CVS du ..2008
  ELSE
    IF (iflag_cvl_sigd==0) THEN
      ! test pour verifier qu on fait la meme chose qu avant: sid constant
      coef(1:ncum) = 1.
    ELSE
      coef(1:ncum) = min(2.*coef(1:ncum), 5.)
      coef(1:ncum) = max(2.*coef(1:ncum), 0.2)
    END IF
  END IF

  IF (prt_level>=20) PRINT *, 'cv3p1_param FIN'
  RETURN
END SUBROUTINE cv3p1_closure

END MODULE cv3p1_closure_mod

