SUBROUTINE cv3_unsat_mf(nloc, ncum, inb, lwork,     &
                        ttk, qtk,                   &
                        cpnk,                       &
                        tvtk, tvtkm1, tvtkp1,       &
                        gzk,                        &
                        pk, pkm1, pkp1, phk, phkp1, &
                        lvk, lfk,                   &
                        sigd,                       &
                        fk,                         &
                        sevk, mltk,                 &
                        Mkp1_in,                    &
                        qkp1, hkp1,                 &
                        Mk_out,                     &
                        qk,  hk                     &
                       )        
                        
IMPLICIT NONE

 include "cvthermo.h"

!  Input variables

  INTEGER, INTENT (IN)                           :: nloc
  INTEGER, INTENT (IN)                           :: ncum
  INTEGER, DIMENSION (nloc), INTENT (IN)         :: inb    ! upmost level of convective column
  LOGICAL, DIMENSION (nloc), INTENT (IN)         :: lwork  ! flag of relevant levels
  REAL, DIMENSION (nloc), INTENT (IN)            :: ttk, qtk ! temperature and specific humidity
                                                                  ! of environment at level k
  REAL, DIMENSION (nloc), INTENT (IN)            :: cpnk        ! Cp of environment air at level k
  REAL, DIMENSION (nloc), INTENT (IN)            :: tvtk    ! virtual temperature of environment at level k
  REAL, DIMENSION (nloc), INTENT (IN)            :: tvtkm1  ! virtual temperature of environment at level k-1
  REAL, DIMENSION (nloc), INTENT (IN)            :: tvtkp1  ! virtual temperature of environment at level k+1
  REAL, DIMENSION (nloc), INTENT (IN)            :: gzk        ! geopotential
  REAL, DIMENSION (nloc), INTENT (IN)            :: pk, pkm1, pkp1   ! mean pressures of layers k, k-1 and k+1
  REAL, DIMENSION (nloc), INTENT (IN)            :: phk, phkp1 ! pressures at interfaces k and k+1
  REAL, DIMENSION (nloc), INTENT (IN)            :: lvk, lfk   ! latent heats of evaporation and fusion at level k
  REAL, DIMENSION (nloc), INTENT (IN)            :: sigd       ! fraction of space occupied by downdrafts

  REAL, DIMENSION (nloc), INTENT (IN)            :: fk         ! fraction of ice in precipitation
  REAL, DIMENSION (nloc), INTENT (IN)            :: sevk       ! product of evaporation at level k by sigs 
                                                               ! (= fraction of precipitation falling in clear air)
  REAL, DIMENSION (nloc), INTENT (IN)            :: mltk       ! melting at level k
  REAL, DIMENSION (nloc), INTENT (IN)            :: Mkp1_in   ! downdraft mass flux at level k+1, positive downward.
  REAL, DIMENSION (nloc), INTENT (IN)            :: qkp1, hkp1  ! specific humidity and 
                                                                      ! static energy of downdrafts at level k+1

! Output variables

  REAL, DIMENSION (nloc), INTENT (OUT)           :: Mk_out     ! downdraft mass flux at level k, positive downward.
  REAL, DIMENSION (nloc), INTENT (OUT)           :: qk, hk     ! specific humidity and 
                                                               ! static energy of downdrafts at level k

! Local variables

  INTEGER                                        :: i          ! horizontal grid cell number
  REAL                                           :: cprime
  REAL                                           :: rpi
  REAL                                           :: tvthk, tvthkp1
  REAL, DIMENSION (nloc)                         :: Mkp1       ! downdraft mass flux at level k+1, negative downward. 
  REAL, DIMENSION (nloc)                         :: Mk         ! downdraft mass flux at level k, negative downward. 
  REAL, DIMENSION (nloc)                         :: htk        !static energy of environment at level k
  REAL, DIMENSION (nloc)                         :: delta_pk   ! pressure thickness of layer k
  REAL, DIMENSION (nloc)                         :: rhohk      ! mass per unit volume at interface k
  REAL, DIMENSION (nloc)                         :: rhohkp1    ! mass per unit volume at interface k+1
  REAL, DIMENSION (nloc)                         :: Ak
  REAL, DIMENSION (nloc)                         :: Ak_crit
  REAL                                           :: gamma      ! rhohk/rhohkp1

  REAL                                           :: p, q, Delta, a, b
  REAL                                           :: alpha, beta
  REAL                                           :: r
  REAL                                           :: awk
  REAL                                           :: c, Mk2

  cprime = cpv - epsim1*cpd
  rpi = 4.*atan(1.)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                On Input, change the signs of amss fluxes
!
         Mkp1(:) = - Mkp1_in(:)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  DO i = 1,ncum
    IF (i<=inb(i) .AND. lwork(i)) THEN
     delta_pk(i) = phk(i) - phkp1(i)
!
! Density at interface k is computed by interpolating the virtual temperatures in layers k-1 and k
! Density at interface k+1 is computed by interpolating the virtual temperatures in layers k and k+1
     tvthk   = (tvtk(i)  *(pkm1(i)-phk(i)) + tvtkm1(i)*(phk(i)  -pk(i)  ))/(pkm1(i)-pk(i))
     rhohk(i) = phk(i)/(rrd*tvthk)
     tvthkp1 = (tvtkp1(i)*(pk(i)-phkp1(i)) + tvtk(i)  *(phkp1(i)-pkp1(i)))/(pk(i)-pkp1(i))
     rhohkp1(i) = phkp1(i)/(rrd*tvthkp1)
!
     htk(i) = cpnk(i)*ttk(i) + gzk(i)
     gamma = rhohk(i)/rhohkp1(i)
!
!     print *,'XXXmf0, tvthk, tvthkp1 ', tvthk, tvthkp1
!     print *,'XXXmf, delta_pk(i), rhohk(i), rhohkp1(i), htk(i) ', delta_pk(i), rhohk(i), rhohkp1(i), htk(i)
!    
     Ak(i) = (sigd(i)**2*delta_pk(i)*rhohk(i)/(cpd*tvtk(i))) * &
      (((cprime*ttk(i)+lvk(i)+fk(i)*lfk(i))*sevk(i) + lfk(i)*mltk(i))* &
       delta_pk(i)*ginv*sigd(i) + &   
       Mkp1(i)*(hkp1(i)-htk(i) - cprime*ttk(i)*(qkp1(i)-qtk(i))))
!
     Ak_crit(i) = Mkp1(i)*Mkp1(i)*Mkp1(i)*(gamma-1.)  
     print *,'XXXmf10, Ak(i), Ak_crit(i), Ak(i)-Ak_crit(i) ', Ak(i), Ak_crit(i), Ak(i)-Ak_crit(i)
!
     IF (Ak(i) .GT. Ak_crit(i)) THEN
       p = - gamma*Mkp1(i)**2
       q = Ak(i)
       Delta = 4.*p*p*p+27.*q*q
   print *,'XXXmf11 p, q, Delta ', p, q, Delta
       IF (Delta .GT.0.) THEN
        a=-q/2.-0.5*sqrt(Delta/27.)
        b=-q/2.+0.5*sqrt(Delta/27.)
!!        alpha=sign(1.,a-b)*(abs(a-b))**(1./3.)
        alpha=sign(1.,a)*(abs(a))**(1./3.)
!!        beta=sign(1.,a+b)*(abs(a+b))**(1./3.)
        beta=sign(1.,b)*(abs(b))**(1./3.)
        IF ((a-b)*(a+b) .LT. 0.) THEN
         Mk(i) = -q/(alpha*alpha-alpha*beta+beta*beta)  
        ELSE
         Mk(i) = alpha + beta  
        ENDIF ! ((a-b)*(a+b) .LT. 0.)
   print *,'XXXmf12 Delta>0, a, b, alpha, beta, Mk(i) ', a, b, alpha, beta, Mk(i)
       ELSE
!!        r = 0.5*sqrt(q*q+abs(Delta))
!!        alpha = Atan(sqrt(abs(Delta))/q)
!!        Mk(i) = 2.*r**(1./3.)*cos(alpha/3.) 
        alpha = Acos(-q/(2.*(-p/3.)**1.5))
        r = sqrt(-4.*p/3.)
        Mk(i) = r*cos(alpha/3.+2.*rpi/3.)

   print *,'XXXmf13 Delta<0, r, alpha, Mk(i) ', r, alpha, Mk(i)
       ENDIF ! (Delta .GT.0.)
   print*,'XXXmf15 verif : Mk(i)**3+p*Mk(i)+q = ', Mk(i)*Mk(i)*Mk(i)+p*Mk(i)+q
     ELSE 
      awk = 0.5
      a = awk*Mkp1(i) 
      b = (1.-awk*(1.+gamma))*Mkp1(i)*Mkp1(i) 
      c = Ak(i)-gamma*(1.-awk)*Mkp1(i)*Mkp1(i)*Mkp1(i)
      Delta = b*b-4.*a*c
   print *,'XXXmf20 a, b, c, Delta ', a, b, c, Delta
      IF (Delta .GT. 0.) THEN
        Mk(i) = (-b+sqrt(Delta))/(2.*a)
        Mk2   = (-b-sqrt(Delta))/(2.*a)
   print *,'XXXmf21 Delta>0 Mk(i), Mk2 ', Mk(i), Mk2
   print*,'XXXmf23 verif : a*Mk(i)*Mk(i)+b*Mk(i)+c = ', a*Mk(i)*Mk(i)+b*Mk(i)+c
      ELSE
        Mk(i) = 0.
   print *,'XXXmf22 Delta<0 Mk(i) ', Mk(i)
      ENDIF ! (Delta .GT. 0.)
    ENDIF ! (Ak(i) .GT. Ak_crit(i))
    ENDIF ! (i<=inb(i) .AND. lwork(i))
  ENDDO ! i = 1,ncum

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                On Output, change the signs of amss fluxes
!
         Mk_out(:) = - Mk(:)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

RETURN
END SUBROUTINE cv3_unsat_mf


