! $Id: cv3p1_closure.F90 2201 2015-02-10 14:23:19Z oboucher $

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, 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  *
  ! **************************************************************

  IMPLICIT NONE

  include "cvthermo.h"
  include "cv3param.h"
  include "YOMCST2.h"
  include "YOMCST.h"
  include "conema3.h"
  include "iniprint.h"

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

  ! input/output:
  REAL sig(nloc, nd), w0(nloc, nd), ptop2(nloc)

  ! output:
  REAL cape(nloc), cin(nloc)
  REAL m(nloc, nd)
  REAL plim1(nloc), plim2(nloc)
  REAL asupmax(nloc, nd), supmax0(nloc)
  REAL asupmaxmin(nloc)
  REAL cbmf(nloc), plfc(nloc)
  REAL wbeff(nloc)
  INTEGER iflag(nloc)

  ! local variables:
  INTEGER il, i, j, k, icbmax, i0(nloc), klfc
  REAL deltap, fac, w, amu
  REAL rhodp
  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 coef(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

  REAL sigmax
  PARAMETER (sigmax=0.1)

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

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


  ! -------------------------------------------------------
  ! -- 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, nl
    DO il = 1, ncum
      m(il, k) = 0.0
    END DO
  END DO

  ! -------------------------------------------------------
  ! -- 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:

  icbmax = 2
  DO il = 1, ncum
    icbmax = max(icbmax, icb(il))
  END DO
  ! if(prt.level.GE.20) print*,'cv3p1_param apres 200'

  ! update sig and w0 below cloud base:

  DO k = 1, icbmax
    DO il = 1, ncum
      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
    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))
          GO TO 425
        END IF ! end IF (P(i) ... )
      END IF ! end IF (icb+1 le i le inb)
    END DO
  END DO

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

    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)

  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
  ! -------------------------------------------------------------

  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):

  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 interval on which cape is computed starts at pbase :

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

      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN

        deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
        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

    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).
  ! ------------------------------------------------------------------------

  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'

  ! c 1.5 Compute cloud base mass flux given by Alp closure (Cbmf1), maximum
  ! c     allowed mass flux (Cbmfmax) and final target mass flux (Cbmf)
  ! c     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

  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 (flag_wb==0) THEN
        wbeff(il) = wbmax
      ELSE IF (flag_wb==1) THEN
        wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
      ELSE IF (flag_wb==2) THEN
        wbeff(il) = wbmax*(0.01*(ph(il,1)-plfc(il)))**2
      END IF
    END IF
  END DO

!CR:Compute k at plfc
  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=k
        endif
     ENDDO
  ENDDO
!RC

  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.) then 
       cbmf1(il) = cbmf1(il) - coef_clos_ls*min(0.,1./RG*omega(il,klfc))
    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_gcm(modname, abort_message, 1)
    END IF
    cbmfmax(il) = sigmax*wb2(il)*100.*p(il, icb(il))/(rrd*tv(il,icb(il)))
  END DO

  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

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

  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)
        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'

  ! 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
  ! c    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