  MODULE SHALLOW_CONV_MOD

!=======================================================================
! --- SHALLOW CONVECTION MODULE - GFDL SPECTRAL MODEL VERSION
!=======================================================================

 use  Sat_Vapor_Pres_Mod, ONLY: compute_qs, lookup_es_des
 use  mpp_mod,            only: input_nml_file
 use  Fms_Mod,            ONLY: FILE_EXIST, ERROR_MESG, FATAL,   &
                                CHECK_NML_ERROR, OPEN_NAMELIST_FILE,      &
                                CLOSE_FILE, mpp_pe, mpp_root_pe, &
                                write_version_number, stdlog

 use constants_mod, only: Hlv, Cp_Air, RDgas, RVgas, Kappa, grav

!---------------------------------------------------------------------
 implicit none
 private
!---------------------------------------------------------------------

 public  :: SHALLOW_CONV, SHALLOW_CONV_INIT
 public  :: MYLCL

!---------------------------------------------------------------------

 character(len=128) :: version = '$Id: shallow_conv.F90,v 19.0 2012/01/06 20:25:24 fms Exp $'
 character(len=128) :: tagname = '$Name: tikal $'

 logical :: module_is_initialized = .false.

!---------------------------------------------------------------------
! --- CONSTANTS
!---------------------------------------------------------------------

  real :: Hlv_by_Cp, Cp_by_RDgas, omkappa, dovcns, d622, d378
  real :: crtkons
  real, parameter :: p00 = 1000.0E2
  integer  :: kctopm1, kctopm2  

  real, allocatable, dimension(:) :: rhcrit, rhmax, rhmin, delrhc

!---------------------------------------------------------------------
! --- NAMELIST
!---------------------------------------------------------------------

  logical ::  lipps    = .false. 
  logical ::  ldetran  = .true. 
  real    ::  theqvcr  =    0.0
  real    ::  pshalow  =  750.0E2
  real    ::  akhsc0   =    5.0
  real    ::  crthum   =    0.85
  real    ::  hc       =    1.0
  integer ::  kctop    =    3

 
 NAMELIST / shallow_conv_nml /    &
         lipps, ldetran,  theqvcr, pshalow, akhsc0, kctop, crthum, hc

!---------------------------------------------------------------------

 contains

!#######################################################################
!#######################################################################

 SUBROUTINE SHALLOW_CONV_INIT( kx )

!=======================================================================
! ***** INITIALIZE SHALLOW CONVECTION
!=======================================================================
!---------------------------------------------------------------------
! Arguments (Intent in)
!     kx     - Number of levels in vertical
!---------------------------------------------------------------------
 integer, intent(in) :: kx
 
!---------------------------------------------------------------------
!  (Intent local)
!---------------------------------------------------------------------
 integer :: unit, io, ierr, logunit
 
!=====================================================================

  if (module_is_initialized) return

!---------------------------------------------------------------------
! --- Read namelist
!---------------------------------------------------------------------

#ifdef INTERNAL_FILE_NML
  read (input_nml_file, nml=shallow_conv_nml, iostat=io)
  ierr = check_nml_error(io,"shallow_conv_nml")
#else
  if( FILE_EXIST( 'input.nml' ) ) then
! -------------------------------------
   unit = OPEN_NAMELIST_FILE ( )
   ierr = 1
   do while( ierr .ne. 0 )
   READ ( unit,  nml = shallow_conv_nml, iostat = io, end = 10 ) 
   ierr = CHECK_NML_ERROR(io,'shallow_conv_nml')
   end do
10 continue
   CALL CLOSE_FILE ( unit )
! -------------------------------------
  end if
#endif

!------- write version number and namelist ---------

  if ( mpp_pe() == mpp_root_pe() ) then
       call write_version_number(version, tagname)
       logunit = stdlog()
       WRITE( logunit, nml = shallow_conv_nml ) 
  endif

!---------------------------------------------------------------------
! --- Initialize constants
!---------------------------------------------------------------------

  d622        = RDgas / RVgas
  d378        = 1.0 - d622
  Hlv_by_Cp   = Hlv / Cp_Air
  Cp_by_RDgas = Cp_Air / RDgas
  omkappa     = 1.0 - Kappa
  dovcns      = -0.5 * grav / RDgas

!----------------------------------------------------

              crtkons = -1.0 * theqvcr * RDgas / grav
  if( lipps ) crtkons = 0.0

!----------------------------------------------------

  kctopm1 = kctop - 1                                                                                           
  kctopm2 = kctop - 2   

!----------------------------------------------------

  allocate( rhcrit(kx) ) 
  allocate(  rhmax(kx) ) 
  allocate(  rhmin(kx) ) 
  allocate( delrhc(kx) ) 

  rhcrit = crthum
  rhmax  = hc                                                    
  rhmin  = 2.0 * rhcrit - hc - 1.0e-15                          
  delrhc = 1.0 / ( rhmax - rhmin )

!-------------------------------------------------------------------
  module_is_initialized = .true.
!---------------------------------------------------------------------

!=====================================================================
  end SUBROUTINE SHALLOW_CONV_INIT

!#######################################################################

  SUBROUTINE SHALLOW_CONV_END
!-------------------------------------------------------------------
  module_is_initialized = .false.
!---------------------------------------------------------------------

  end SUBROUTINE SHALLOW_CONV_END

!#######################################################################
  SUBROUTINE SHALLOW_CONV( Temp, qmix0, pfull, phalf, akhsc, kbot )

!=======================================================================
! --- SHALLOW CONVECTION
!=======================================================================
!----------------------------------------------------------------------
! Arguments (Intent in)
!       Temp    -  Temperature
!       qmix0   -  Specific humidity
!       pfull   -  Pressure at full levels
!       phalf   -  Pressure at half levels
!       kbot    -  OPTIONAL; lowest model level index (integer)
!----------------------------------------------------------------------
  real, intent(in), dimension(:,:,:) :: Temp, qmix0, pfull, phalf

  integer, intent(in), OPTIONAL, dimension(:,:) :: kbot

!----------------------------------------------------------------------
! Arguments (Intent out)
!       akhsc  -  mixing coefficient for heat and moisture
!                 due to shallow convection
!----------------------------------------------------------------------
  real, intent(out), dimension(:,:,:) :: akhsc

!----------------------------------------------------------------------
! --- local
!----------------------------------------------------------------------
  real,    dimension(SIZE(Temp,1),SIZE(Temp,2)) ::                &
           plcl, rhumjmp, rhumscl,  xy1, xy2, xy3

  integer, dimension(SIZE(Temp,1),SIZE(Temp,2)) ::                &
           ksiglcl

  real,    dimension(SIZE(Temp,1),SIZE(Temp,2),SIZE(Temp,3)) ::   &
           qmix,  dphalf, qsat,  qsat2, rhum, theta, thetav, buoy, xyz1 

  integer, dimension(SIZE(Temp,1),SIZE(Temp,2),SIZE(Temp,3)) ::   &
           kbuoy 

 integer :: k, kx, kxm, kxp, i, ix, j, jx

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

  ix  = SIZE(Temp,1)
  jx  = SIZE(Temp,2)
  kx  = SIZE(Temp,3)

  kxm = kx - 1
  kxp = kx + 1

!=======================================================================
! --- MOISTURE VARIABLES
!=======================================================================

  qmix = qmix0
  qmix = MAX( qmix, 1.0E-6 )
  qmix = MIN( qmix, 0.2    )

! --- saturation mixing ratio 
  call compute_qs (Temp, pfull, qsat)

! --- relative humidity
  call compute_qs (Temp, pfull, qsat2, q=qmix)
  rhum = qmix / qsat2

!=======================================================================
! --- POTENTIAL TEMPERATURE
!=======================================================================

  theta  = Temp  * ( ( p00 / pfull )**Kappa )

!=======================================================================
! --- CALCULATE THE LIFTING CONDENSATION LEVEL, IE CLOUB BASE
!=======================================================================

  if( PRESENT( kbot ) ) then
     do j = 1,jx
     do i = 1,ix
        k = kbot(i,j)
             xy1(i,j) =      Temp(i,j,k)
             xy2(i,j) = MIN( qmix(i,j,k), qsat(i,j,k) )
             xy3(i,j) =     pfull(i,j,k)
     end do
     end do
  else
             xy1(:,:) =      Temp(:,:,kx)
             xy2(:,:) = MIN( qmix(:,:,kx), qsat(:,:,kx) )
             xy3(:,:) =     pfull(:,:,kx)
  end if

  CALL MYLCL( xy1, xy2, xy3, phalf, plcl, ksiglcl )

!=======================================================================
! --- INITALIZE
!=======================================================================

  kbuoy   = kxp
  akhsc   = 0.0
 
!=======================================================================
! --- BUOYANCY
!=======================================================================

!---------------------------------------------------------------------
! --- DEFAULT:
! --- BASED ON EQUIVALENT POTENTIAL TEMPERATURE GRADIENT 
!---------------------------------------------------------------------

  if( .not. lipps ) then
! %%%%%%%%%%%%%%%%%%%%%%%

! --- Vertical differential of pressure 
  dphalf(:,:,1:kxm) = pfull(:,:,2:kx) - pfull(:,:,1:kxm)

  if( PRESENT( kbot ) ) then
  dphalf(:,:,1:kxm) = MAX( dphalf(:,:,1:kxm), 1.0e-5 )
  end if

! --- Equivalent potential temperature
  xyz1    = ( Hlv_by_Cp * qmix ) / Temp
! thetav = theta * (1.0 + xyz1 )
  thetav = theta * EXP( xyz1 )

! --- Equivalent potential temperature gradient
  xyz1(:,:,2:kx)  =     thetav(:,:,2:kx)  - thetav(:,:,1:kxm) 
  xyz1(:,:,2:kx)  =       xyz1(:,:,2:kx)  / dphalf(:,:,1:kxm)
  buoy(:,:,2:kxm) = 0.5*( xyz1(:,:,2:kxm) +   xyz1(:,:,3:kx) )

! %%%%%%%%%%%%%%%%%%%%%%%
  endif

!---------------------------------------------------------------------
! --- OPTION:
! --- BUOYANCY ALA FRANK LIPPS 
!---------------------------------------------------------------------

  if( lipps ) then
! %%%%%%%%%%%%%%%%%%%%%%%

! --- Virtual potential temperature gardient
  thetav = theta * ( 1.0 + 0.608 * qmix )

! --- Buoyancy
  do k = kctopm2,kxm
! -------------------
  xy1(:,:) = Hlv / ( Cp_Air*RVgas*Temp(:,:,k)*Temp(:,:,k) ) + 1.0 / qsat(:,:,k)

  xy2(:,:) = ( grav / ( Cp_Air*Temp(:,:,k) ) ) *                              &  
             ( Hlv / ( RVgas*Temp(:,:,k) ) - Cp_by_RDgas ) / xy1(:,:)

  xy3(:,:) = ( thetav(:,:,k+1) - thetav(:,:,k )  ) / phalf(:,:,k+1) +     &
             ( thetav(:,:,k  ) - thetav(:,:,k-1) ) / phalf(:,:,k  )

  buoy(:,:,k) = ( hlv / ( Cp_Air*Temp(:,:,k) ) - 1.608 ) * xy2(:,:)
  buoy(:,:,k) =    buoy(:,:,k) - dovcns*pfull(:,:,k) * xy3(:,:) / Temp(:,:,k)
! -------------------
  end do

! %%%%%%%%%%%%%%%%%%%%%%%
  endif

!=======================================================================
! --- COMPUTE THE LEVEL OF NO BUOYANCY 
! --- RETAIN ONLY THE LOWEST CONTIGUOUS BUOYANT SHALLOW CONVECTIVE LAYER.
!=======================================================================

  do k = kctopm1,kxm
! -------------------
  where ( ( pfull(:,:,k) >= pshalow   ) .and.    &
          ( pfull(:,:,k) <= plcl(:,:) ) .and.    &
          (  buoy(:,:,k) >= crtkons   ) )  
            kbuoy(:,:,k) =  k 
  endwhere
! -------------------
  end do

  do k = kctopm1,kxm
! -------------------
  where( ( pfull(:,:,k)   <  plcl(:,:) ) .and.   &
         ( kbuoy(:,:,k)   == kxp       ) .and.   &
         ( kbuoy(:,:,k-1) == k-1       ) )   
           kbuoy(:,:,k-1) =  kxp
  endwhere
! -------------------
  end do

!=======================================================================
! --- SHALLOW CONVECTION WILL OCCUR AT LEVELS WHERE KBUOY <= KSIGLCL
!=======================================================================

  do k = kctopm1,kxm
! -------------------
  where( kbuoy(:,:,k) <= ksiglcl(:,:) ) 
         akhsc(:,:,k+1) =  akhsc0
  endwhere
! -------------------
  end do

!=======================================================================
! --- DETRAINMENT THRU INVERSION LAYER
!=======================================================================

!---------------------------------------------------------------------
! --- DEFAULT:
! --- ENHANCED DETRAINMENT THRU INVERSION LAYER.
!---------------------------------------------------------------------

  if( ldetran) then
! %%%%%%%%%%%%%%%%%%%%%%%%

  do k = kctopm1,kxm
! -------------------
  where( ( kbuoy(:,:,k)   == k       ) .and.   &
         ( kbuoy(:,:,k-1) == kxp     ) .and.   &
         ( pfull(:,:,k)   >= pshalow ) ) 
           akhsc(:,:,k)   =  0.2 * akhsc0
           akhsc(:,:,k+1) =  0.6 * akhsc0
  endwhere
! -------------------
  end do

  do k = kctopm1,kxm
! -------------------
  where( ( pfull(:,:,k)   <= plcl(:,:) ) .and.  &
         ( pfull(:,:,k+1) >  plcl(:,:) ) .and.  &
         ( kbuoy(:,:,k)   == k         )  )
           akhsc(:,:,k+1) =  0.2 * akhsc0
  endwhere
! -------------------
  end do

! %%%%%%%%%%%%%%%%%%%%%%%%
  endif

!---------------------------------------------------------------------
! --- OPTION:
! --- NORMAL DETRAINMENT THRU INVERSION LAYER
!---------------------------------------------------------------------

  if( .not. ldetran ) then
! %%%%%%%%%%%%%%%%%%%%%%%%

  rhumscl = 0.0
  rhumjmp = 0.0

  do k = kctopm1,kxm
! -------------------
  where( ( kbuoy(:,:,k)   == k       ) .and.    &
         (  buoy(:,:,k-1) <  crtkons ) .and.    &
         (  buoy(:,:,k)   >= crtkons ) ) 
         rhumjmp(:,:) =   rhum(:,:,k) -   rhum(:,:,k-1)
         rhumscl(:,:) =     delrhc(k) * ( rhum(:,:,k) - rhmin(k) )
  endwhere
! -------------------
  end do

  rhumscl = MIN( 1.0, rhumscl )
  rhumscl = MAX( 0.0, rhumscl )
  rhumjmp = MIN( 1.0, rhumjmp )
  rhumjmp = MAX( 0.0, rhumjmp )

  do k = kctopm1,kxm
! -------------------
  where( ( kbuoy(:,:,k)   == k       ) .and.   &
         ( kbuoy(:,:,k-1) == kxp     ) .and.   &
         ( pfull(:,:,k)   >= pshalow ) ) 
           akhsc(:,:,k)   =  akhsc0 * rhumscl(:,:) * rhumjmp(:,:)
  endwhere
! -------------------
  end do

! %%%%%%%%%%%%%%%%%%%%%%%%
  endif

!=======================================================================
! --- CONFINE SHALLOW CONVECTION TO ( pshalow <= p <= plcl )
!=======================================================================

  do k = kctopm1,kxm
! -------------------
  where ( ( pfull(:,:,k) <= pshalow   ) .or.    &
          ( pfull(:,:,k) >= plcl(:,:) ) )  
            akhsc(:,:,k+1) =  0.0 
  endwhere
! -------------------
  end do

!=======================================================================
  end SUBROUTINE SHALLOW_CONV

!#######################################################################

  SUBROUTINE MYLCL ( tlparc, qlparc, plparc, phalf, plcl, kbase )

!=======================================================================
! ***** COMPUTE LCL ( CLOUD BASE )
!=======================================================================
!---------------------------------------------------------------------
! Arguments (Intent in)
!       tlparc   Initial parcel temperature
!       qlparc   Initial parcel mixing ratio
!      plparc   Initial parcel pressure
!       phalf    Pressure at half levels
! Arguments (Intent out)
!       plcl     Pressure at LCL
!       kbase    Index of LCL in column
!---------------------------------------------------------------------
  real,    intent(in),  dimension(:,:)   :: tlparc, qlparc, plparc
  real,    intent(in),  dimension(:,:,:) :: phalf
  real,    intent(out), dimension(:,:)   :: plcl
  integer, intent(out), dimension(:,:)   :: kbase

!---------------------------------------------------------------------
!  (Intent local)
!---------------------------------------------------------------------
  integer, parameter :: iter_max = 10      
  real,    parameter :: small    = 1.0E-2 

  real,    dimension(size(tlparc,1),size(tlparc,2)) ::    &
           tlcl, tlclo, clclo, esat, esato, desato, xy1, xy2

 logical, dimension(size(tlparc,1),size(tlparc,2)) ::  &
           non_cnvg

  integer :: k, kx, n, iter

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

! --- Index of lowest model level, etc
  kx  = size( phalf, 3 ) - 1

! --- Initial guess for temperature at LCL
  tlclo = tlparc

! --- Compute constant factor
  clclo = ( 1.0 + d622/qlparc ) / plparc
  clclo = kappa*LOG( clclo )
  clclo =       EXP( clclo )
  clclo =    tlclo * clclo

! --- Start with all points non-convergent
  non_cnvg = .true.

! $$$$$$$$$$$$$$$$$$$$
  do iter = 1,iter_max
! $$$$$$$$$$$$$$$$$$$$

! --- Compute saturation vapor pressure and derivative
  call lookup_es_des (tlclo, esato, desato)

! --- Compute new guess for temperature at LCL
  where (non_cnvg)
     xy1  = kappa * clclo * desato
     xy2  = omkappa*LOG( esato )
     xy2  = EXP(  xy2 )
     tlcl = ( xy1 * tlclo - clclo * esato ) / ( xy1 - xy2 )
     xy2  = abs( tlcl - tlclo )
  end where

! --- Test for convergence
  where (non_cnvg .and. xy2 <= small)
     esat = esato
     non_cnvg = .false.
  endwhere
  n = COUNT( non_cnvg )

  if( n .eq. 0 ) go to 1000

! --- Shift for next iteration
  tlclo = tlcl

! $$$$$$$$$$$$$$$$$$$$
  end do
! $$$$$$$$$$$$$$$$$$$$
       CALL ERROR_MESG ('MYLCL in SHALLOW_CONV_MOD',  &
                        'ITERATION LOOP FOR LCL FAILED', FATAL)
 1000 continue

! --- Compute pressure at LCL
  plcl = ( 1.0 + d622 / qlparc ) * esat

! --- Bound plcl 
  plcl(:,:) = MAX( plcl(:,:), pshalow      )
  plcl(:,:) = MIN( plcl(:,:),  plparc(:,:) )

! --- Find index of LCL
  do k = 2,kx
  where ( ( plcl(:,:) >= phalf(:,:,k  ) ) .and.   &
          ( plcl(:,:) <= phalf(:,:,k+1) ) )  
           kbase(:,:) = k
  end where
  end do

!=======================================================================
  end SUBROUTINE MYLCL

!#######################################################################
!#######################################################################
  end MODULE SHALLOW_CONV_MOD

