module convection

  use vtype_module

  implicit none

  public :: ConvCalcRadTendency
  public :: ConvCalcRadConvTendency

contains
  !----------------------------------------------------------------------------

  subroutine ConvCalcRadTendencyExplicit( &
    & Grav, Cp, SurfHeatCap, DelTime, &
    & kmax, r_Press, r_Temp, &
    & r_RadUwFlux, r_RadDwFlux, &
    & r_DelTemp, DelSurfTemp &
    & )

    real(dp), intent(in ) :: Grav
    real(dp), intent(in ) :: Cp
    real(dp), intent(in ) :: SurfHeatCap
    real(dp), intent(in ) :: DelTime
    integer , intent(in ) :: kmax
    real(dp), intent(in ) :: r_Press        (0:kmax)
    real(dp), intent(in ) :: r_Temp         (0:kmax)
    real(dp), intent(in ) :: r_RadUwFlux    (0:kmax)
    real(dp), intent(in ) :: r_RadDwFlux    (0:kmax)
    real(dp), intent(out) :: r_DelTemp      (0:kmax)
    real(dp), intent(out) :: DelSurfTemp


    ! Local variables
    !
    real(dp) :: z_DelAtmMass (1:kmax)
    real(dp) :: r_RadFlux    (0:kmax)
    integer  :: k


    r_RadFlux     = r_RadUwFlux     - r_RadDwFlux


    do k = 1, kmax
      z_DelAtmMass(k) = ( r_Press(k-1) - r_Press(k) ) / Grav
    end do


    DelSurfTemp = - r_RadFlux(0) * DelTime / SurfHeatCap

    do k = 1, kmax
      r_DelTemp(k) = &
        & - ( r_RadFlux(k) - r_RadFlux(k-1) ) * DelTime &
        &   / ( Cp * z_DelAtmMass(k) )
    end do

    ! T0 = ( T2 - T1 ) / ( 2 - 1 ) * ( 0 - 1 ) + T1
    !    = - ( T2 - T1 ) + T1
    !    = 2 * T1 - T2
    r_DelTemp(0) = &
      &   2.0d0 * ( r_Temp(1) + r_DelTemp(1) ) - ( r_Temp(2) + r_DelTemp(2) ) &
      & - r_Temp(0)


  end subroutine ConvCalcRadTendencyExplicit

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

  subroutine ConvCalcRadTendency( &
    & Grav, Cp, SurfHeatCap, DelTime, &
    & kmax, r_Press, r_Temp, &
    & r_RadUwFlux, r_RadDwFlux, ra_DelRadUwFlux, ra_DelRadDwFlux, &
    & r_DelTemp, DelSurfTemp &
    & )

    real(dp), intent(in ) :: Grav
    real(dp), intent(in ) :: Cp
    real(dp), intent(in ) :: SurfHeatCap
    real(dp), intent(in ) :: DelTime
    integer , intent(in ) :: kmax
    real(dp), intent(in ) :: r_Press        (0:kmax)
    real(dp), intent(in ) :: r_Temp         (0:kmax)
    real(dp), intent(in ) :: r_RadUwFlux    (0:kmax)
    real(dp), intent(in ) :: r_RadDwFlux    (0:kmax)
    real(DP), intent(in ) :: ra_DelRadUwFlux(0:kmax,-1:1)
    real(DP), intent(in ) :: ra_DelRadDwFlux(0:kmax,-1:1)
    real(dp), intent(out) :: r_DelTemp      (0:kmax)
    real(dp), intent(out) :: DelSurfTemp


    ! Local variables
    !
    real(dp) :: z_DelAtmMass (1:kmax)
    real(dp) :: r_RadFlux    (0:kmax)
    real(DP) :: ra_DelRadFlux(0:kmax,-1:1)
    real(dp) :: rr_Mat(-1:kmax,-1:kmax)
    real(dp) :: r_Vec (-1:kmax)
    integer  :: k
    integer  :: kr
    integer  :: kc


    r_RadFlux     = r_RadUwFlux     - r_RadDwFlux
    ra_DelRadFlux = ra_DelRadUwFlux - ra_DelRadDwFlux


    do k = 1, kmax
      z_DelAtmMass(k) = ( r_Press(k-1) - r_Press(k) ) / Grav
    end do


    rr_Mat = 0.0d0
    r_Vec  = 0.0d0
    !
    ! Ts
    kr = -1
    !
    kc = -1
    rr_Mat(kr,kc) = SurfHeatCap / DelTime + ra_DelRadFlux(0,kc)
    kc = 0
    rr_Mat(kr,kc) = + ra_DelRadFlux(0,kc)
    kc = 1
    rr_Mat(kr,kc) = + ra_DelRadFlux(0,kc)
    !
    r_Vec(kr)     = - r_RadFlux(0)
    !
    ! T (k=0)
    !   T0 = ( T2 - T1 ) / ( 2 - 1 ) * ( 0 - 1 ) + T1
    !      = - ( T2 - T1 ) + T1
    !      = 2 * T1 - T2
    !   T0 - 2 * T1 + T2 = 0
    !   r_Temp(0) + r_DelTemp(0) - 2 * ( r_Temp(1) + r_DelTemp(1) )
    !     + r_Temp(2) + r_DelTemp(2) = 0
    !   r_DelTemp(0) - 2 * r_DelTemp(1) + r_DelTemp(2) 
    !     = - ( r_Temp(0) - 2 * r_Temp(1) + r_Temp(2) )
    kr = 0
    kc =   0
    rr_Mat(kr,kc) =   1.0d0
    kc =   1
    rr_Mat(kr,kc) = - 2.0d0
    kc =   2
    rr_Mat(kr,kc) =   1.0d0
    !
    r_Vec(kr)     = - ( r_Temp(0) - 2.0d0 * r_Temp(1) + r_Temp(2) )
    !
    ! T
    do kr = 1, kmax
      kc = -1
      rr_Mat(kr,kc) = - ( ra_DelRadFlux(kr-1,kc) - ra_DelRadFlux(kr,kc) )
      kc = 0
      rr_Mat(kr,kc) = - ( ra_DelRadFlux(kr-1,kc) - ra_DelRadFlux(kr,kc) )
      kc = 1
      rr_Mat(kr,kc) = - ( ra_DelRadFlux(kr-1,kc) - ra_DelRadFlux(kr,kc) )
      kc = kr
      rr_Mat(kr,kc) = rr_Mat(kr,kc) + z_DelAtmMass(kr) * Cp / DelTime
      !
      r_Vec(kr) = r_RadFlux(kr-1) - r_RadFlux(kr)
    end do

    call ludecomp_prep_simple ( kmax+2, rr_Mat )
    call ludecomp_solve_simple( kmax+2, rr_Mat, r_Vec )

    r_DelTemp   = r_Vec(0:kmax)
    DelSurfTemp = r_Vec(-1)


  end subroutine ConvCalcRadTendency

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

  subroutine ConvCalcRadConvTendency( &
    & FlagMoistConvection, &
    & Grav, Cp, SurfHeatCap, DelTime, &
    & kmax, r_Press, r_Temp, SurfTemp, r_VMRH2OVap, r_MeanMolWt, &
    & r_RadUwFlux, r_RadDwFlux, ra_DelRadUwFlux, ra_DelRadDwFlux, &
    & r_DelTemp, DelSurfTemp, r_DelVMRH2OVap, SurfSensFlux &
    & )

    logical , intent(in ) :: FlagMoistConvection

    real(dp), intent(in ) :: Grav
    real(dp), intent(in ) :: Cp
    real(dp), intent(in ) :: SurfHeatCap
    real(dp), intent(in ) :: DelTime
    integer , intent(in ) :: kmax
    real(dp), intent(in ) :: r_Press        (0:kmax)
    real(dp), intent(in ) :: r_Temp         (0:kmax)
    real(dp), intent(in ) :: SurfTemp
    real(dp), intent(in ) :: r_VMRH2OVap    (0:kmax)
    real(DP), intent(in ) :: r_MeanMolWt    (0:kmax)
    real(dp), intent(in ) :: r_RadUwFlux    (0:kmax)
    real(dp), intent(in ) :: r_RadDwFlux    (0:kmax)
    real(DP), intent(in ) :: ra_DelRadUwFlux(0:kmax,-1:1)
    real(DP), intent(in ) :: ra_DelRadDwFlux(0:kmax,-1:1)
    real(dp), intent(out) :: r_DelTemp      (0:kmax)
    real(dp), intent(out) :: DelSurfTemp
    real(DP), intent(out) :: r_DelVMRH2OVap (0:kmax)
    real(DP), intent(out) :: SurfSensFlux


    ! Local variables
    !
    real(DP) :: z_VertDelTemp  (1:kmax)
    real(DP) :: r_TempTentative(0:kmax)
    real(DP) :: z_RadFluxConv  (1:kmax)

    ! Full level index (mid-poing of layers)
    integer , parameter :: kConvS = 1
    integer             :: kConvE
    integer             :: kConvEMin

    real(dp) :: z_DelAtmMass (1:kmax)
    real(dp) :: TempTentative

    integer             :: k


    do k = 1, kmax
      z_RadFluxConv(k) = - ( r_RadUwFlux(k) - r_RadDwFlux(k) - ( r_RadUwFlux(k-1) - r_RadDwFlux(k-1) ) )
    end do

    ! Find top of convection
    do kConvEMin = 1, kmax
      if ( z_RadFluxConv(kConvEMin) < 0.0d0 ) exit
    end do
    if ( kConvEMin > kmax ) then
      write( 6, * ) 'Convection seems to reach top of model'
      stop
    end if


!!$    do kConvE = 1, kmax
    do kConvE = kConvEMin, kmax

      call ConvCalcRadConvTendencyCore( &
        & FlagMoistConvection, &
        & kConvS, kConvE, &
        & Grav, Cp, SurfHeatCap, DelTime, &
        & kmax, r_Press, r_Temp, SurfTemp, r_VMRH2OVap, r_MeanMolWt, &
        & r_RadUwFlux, r_RadDwFlux, ra_DelRadUwFlux, ra_DelRadDwFlux, &
        & r_DelTemp, DelSurfTemp, r_DelVMRH2OVap, SurfSensFlux &
        & )
!!$      call ConvCalcRadConvTendencyCoreTrial( &
!!$        & FlagMoistConvection, &
!!$        & kConvS, kConvE, &
!!$        & Grav, Cp, SurfHeatCap, DelTime, &
!!$        & kmax, r_Press, r_Temp, SurfTemp, r_VMRH2OVap, r_MeanMolWt, &
!!$        & r_RadUwFlux, r_RadDwFlux, ra_DelRadUwFlux, ra_DelRadDwFlux, &
!!$        & r_DelTemp, DelSurfTemp, r_DelVMRH2OVap, SurfSensFlux &
!!$        & )

      if ( kConvE < kmax ) then
        r_TempTentative = r_Temp + r_DelTemp
        call ConvCalcMoistConvDelTemp( &
          & FlagMoistConvection, &
          & kmax, r_Press, r_MeanMolWt, r_TempTentative, &
          & z_VertDelTemp &
          & )
        if (    r_TempTentative(kConvE)+z_VertDelTemp(kConvE+1) &
          &  <= r_TempTentative(kConvE+1) ) then
          exit
        end if
      end if
    end do

    write( 6, * ) "Convection : ", kConvE, r_Press(kConvE)


    do k = 1, kmax
      z_DelAtmMass(k) = ( r_Press(k-1) - r_Press(k) ) / Grav
    end do

    do k = kConvE+1, kmax-1
      r_TempTentative = r_Temp + r_DelTemp
      call ConvCalcMoistConvDelTemp( &
!!$        & FlagMoistConvection, &
!!$        & kmax, r_Press, r_MeanMolWt, r_TempTentative, &
!!$        & z_VertDelTemp &
        & .false., &
        & 1, r_Press(k:k+1), r_MeanMolWt(k:k+1), r_TempTentative(k:k+1), &
        & z_VertDelTemp(k+1:k+1) &
        & )
!!$      if (    r_TempTentative(kConvE)+z_VertDelTemp(kConvE+1) &
!!$        &  >  r_TempTentative(kConvE+1) ) then
      if (    r_TempTentative(k  )+z_VertDelTemp(k+1) &
        &  >  r_TempTentative(k+1) ) then
        write( 6, * ) "  Convection in free atmosphere: ", k

        ! Cp dM_k * T_k + Cp dM_{k+1} * T_{k+1}
        ! = Cp dM_k * T'_k + Cp dM_{k+1} * T'_{k+1}
        !
        ! T'_{k+1} = T'_{k} + VdT_{k+1}
        !
        ! Cp dM_k * T_k + Cp dM_{k+1} * T_{k+1}
        ! = Cp dM_k * T'_k + Cp dM_{k+1} * ( T'_{k} + VdT_{k+1} )
        ! = ( Cp dM_k + Cp dM_{k+1} ) * T'_k + Cp dM_{k+1} * VdT_{k+1}
        ! ( Cp dM_k + Cp dM_{k+1} ) * T'_k
        ! = Cp dM_k * T_k + Cp dM_{k+1} * T_{k+1} - Cp dM_{k+1} * VdT_{k+1}
        ! = Cp dM_k * T_k + Cp dM_{k+1} * ( T_{k+1} - VdT_{k+1} )
        ! T'_k
        ! = ( Cp dM_k * T_k + Cp dM_{k+1} * ( T_{k+1} - VdT_{k+1} ) )
        !   / ( Cp dM_k + Cp dM_{k+1} )
        !
        TempTentative = &
          & + (   Cp * z_DelAtmMass(k  ) * r_TempTentative(k) &
          &     + Cp * z_DelAtmMass(k+1) * ( r_TempTentative(k+1) - z_VertDelTemp(k+1) ) ) &
          & / ( Cp * z_DelAtmMass(k  ) + Cp * z_DelAtmMass(k+1) )
        r_DelTemp(k  ) = r_DelTemp(k  ) &
          & + TempTentative - r_TempTentative(k)
        r_DelTemp(k+1) = r_DelTemp(k+1) &
          & + TempTentative + z_VertDelTemp(k+1) - r_TempTentative(k+1)
      end if
    end do


  end subroutine ConvCalcRadConvTendency

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

  subroutine ConvCalcRadConvTendencyCore( &
    & FlagMoistConvection, &
    & kConvS, kConvE, &
    & Grav, Cp, SurfHeatCap, DelTime, &
    & kmax, r_Press, r_Temp, SurfTemp, r_VMRH2OVap, r_MeanMolWt, &
    & r_RadUwFlux, r_RadDwFlux, ra_DelRadUwFlux, ra_DelRadDwFlux, &
    & r_DelTemp, DelSurfTemp, r_DelVMRH2OVap, SurfSensFlux &
    & )

    logical , intent(in ) :: FlagMoistConvection
    ! Full level index (mid-poing of layers)
    integer , intent(in ) :: kConvS
    integer , intent(in ) :: kConvE
    real(dp), intent(in ) :: Grav
    real(dp), intent(in ) :: Cp
    real(dp), intent(in ) :: SurfHeatCap
    real(dp), intent(in ) :: DelTime
    integer , intent(in ) :: kmax
    real(dp), intent(in ) :: r_Press        (0:kmax)
    real(dp), intent(in ) :: r_Temp         (0:kmax)
    real(dp), intent(in ) :: SurfTemp
    real(dp), intent(in ) :: r_VMRH2OVap    (0:kmax)
    real(DP), intent(in ) :: r_MeanMolWt    (0:kmax)
    real(dp), intent(in ) :: r_RadUwFlux    (0:kmax)
    real(dp), intent(in ) :: r_RadDwFlux    (0:kmax)
    real(DP), intent(in ) :: ra_DelRadUwFlux(0:kmax,-1:1)
    real(DP), intent(in ) :: ra_DelRadDwFlux(0:kmax,-1:1)
    real(dp), intent(out) :: r_DelTemp      (0:kmax)
    real(dp), intent(out) :: DelSurfTemp
    real(DP), intent(out) :: r_DelVMRH2OVap (0:kmax)
    real(DP), intent(out) :: SurfSensFlux


    ! Function
    !
!!$    real(DP) :: LOCALPressH2OVapSat
!!$    real(DP) :: LOCALDPressH2OVapSatDTemp

    ! Local variables
    !
    real(DP), parameter :: LatentHeat = 2.5d6    ! J kg-1
    real(DP), parameter :: CpH2OVap   = 1810.0d0 ! J K-1 kg-1
    real(DP), parameter :: MolWtH2O   = 18.0d-3  ! kg mol-1
    real(dp)            :: LatentHeatPerMol


    real(dp) :: RelHumid = 1.0d0
    real(dp) :: z_DelAtmMass (1:kmax)
    real(dp) :: r_PressH2OVapSat      (0:kmax)
    real(dp) :: r_DPressH2OVapSatDTemp(0:kmax)
    real(dp) :: r_DMMRH2OVapSatDTemp  (0:kmax)
    real(dp) :: z_VertDelTemp(1:kmax)
    real(dp) :: r_RadFlux    (0:kmax)
    real(DP) :: ra_DelRadFlux(0:kmax,-1:1)
    real(dp) :: rr_Mat(-2:kmax,-2:kmax)
    real(dp) :: r_Vec (-2:kmax)

    real(dp) :: VMRH2OVapTentative

    integer  :: k
    integer  :: kr
    integer  :: kc


    LatentHeatPerMol = LatentHeat * MolWtH2O

    r_RadFlux     = r_RadUwFlux     - r_RadDwFlux
    ra_DelRadFlux = ra_DelRadUwFlux - ra_DelRadDwFlux


    do k = 1, kmax
      z_DelAtmMass(k) = ( r_Press(k-1) - r_Press(k) ) / Grav
    end do

    if ( FlagMoistConvection ) then
      do k = 0, kmax
        r_PressH2OVapSat(k) = &
          & LOCALPressH2OVapSat( LatentHeatPerMol, r_Temp(k) )
      end do
      do k = 0, kmax
        r_DPressH2OVapSatDTemp(k) = &
          & LOCALDPressH2OVapSatDTemp( LatentHeatPerMol, r_Temp(k) )
      end do
    else
      r_PressH2OVapSat = 0.0d0
      r_DPressH2OVapSatDTemp = 0.0d0
    end if
    r_DMMRH2OVapSatDTemp = r_DPressH2OVapSatDTemp / r_Press &
      & * MolWtH2O / r_MeanMolWt


    call ConvCalcMoistConvDelTemp( &
      & FlagMoistConvection, &
      & kmax, r_Press, r_MeanMolWt, r_Temp, &
      & z_VertDelTemp &
      & )


    rr_Mat = 0.0d0
    r_Vec  = 0.0d0
    !
    ! Energy conservation
    kr = -2
    kc = -2
    rr_Mat(kr,kc) = -1.0d0
    kc = -1
    rr_Mat(kr,kc) = + ( ra_DelRadFlux(kConvE,kc) - ra_DelRadFlux(0,kc) )
    kc =  0
    rr_Mat(kr,kc) = + ( ra_DelRadFlux(kConvE,kc) - ra_DelRadFlux(0,kc) )
    kc =  1
    rr_Mat(kr,kc) = + ( ra_DelRadFlux(kConvE,kc) - ra_DelRadFlux(0,kc) )
    do kc = 1, kConvE
!!$      rr_Mat(kr,kc) = rr_Mat(kr,kc) &
!!$        & + z_DelAtmMass(kc) / DelTime * ( Cp + RelHumid * LatentHeat * r_DMMRH2OVapSatDTemp(kc) )
      ! Latent heat is neglected to conserve energy
      rr_Mat(kr,kc) = rr_Mat(kr,kc) &
        & + z_DelAtmMass(kc) / DelTime * ( Cp + 0.0d0 )
    end do
    !
    r_Vec(kr) = - ( r_RadFlux(kConvE) - r_RadFlux(0) )

    !
    ! Ts
    kr = -1
    kc = -2
    rr_Mat(kr,kc) = 1.0d0
    kc = -1
    rr_Mat(kr,kc) = SurfHeatCap / DelTime + ra_DelRadFlux(0,kc)
    kc =  0
    rr_Mat(kr,kc) = + ra_DelRadFlux(0,kc)
    kc =  1
    rr_Mat(kr,kc) = + ra_DelRadFlux(0,kc)
    !
    r_Vec(kr) = - r_RadFlux(0)
    !
    ! T (k=0)
    !   T0 = Ts
    !   r_Temp(0) + r_DelTemp(0) = SurfTemp + DelSurfTemp
    !   - DelSurfTemp + r_DelTemp(0) = SurfTemp - r_Temp(0)
    kr = 0
    kc =  -1
    rr_Mat(kr,kc) = -1.0d0
    kc =   0
    rr_Mat(kr,kc) =  1.0d0
    !
    r_Vec(kr)     = SurfTemp - r_Temp(0)
    !
    ! T (1 <= k < kConvE)
    do kr = 1, kConvE
      kc = kr-1
      rr_Mat(kr,kc) = -1.0d0
      kc = kr
      rr_Mat(kr,kc) =  1.0d0
      !
!!$      if ( kConvE == 1 ) then
!!$        ! convective layer is assumed to be very thin
!!$        r_Vec(kr) = r_Temp(kr-1) &
!!$          & + ( r_Temp(2) - r_Temp(1) ) &
!!$          &   / log( r_Press(2) / r_Press(1) ) &
!!$          &   * log( r_Press(1) / r_Press(0) ) &
!!$          & - r_Temp(kr)
!!$      else
        ! convective layer is thicker than 1 layer
        r_Vec(kr) = r_Temp(kr-1) + z_VertDelTemp(kr) - r_Temp(kr)
!!$      end if
    end do



    ! test
!!$    do kr = 1, kConvE
!!$      ! general equation
!!$      ! - DT_{kr-1} + DT_{kr} = T_{kr-1} + vdelT_{kr} - T_{kr}
!!$      !  T_{kr} + DT_{kr} = T_{kr-1} + DT_{kr-1} + vdelT_{kr}
!!$      !
!!$      ! new formulation
!!$      ! vdelT_{kr} =
!!$      !    dT/dlogp_r * ( log(p_{kr}) - logp )
!!$      !  + dT/dlogp_c * ( logp - log(p_{kr-1}) )
!!$      !
!!$      !  T_{kr-1} + DT_{kr-1} + dT/dlogp_c * ( logp - log(p_{kr-1}) )
!!$      !    T_{kr} + DT_{kr} - dT/dlogp_r * ( log(p_{kr}) - logp )
!!$      !
!!$      kc = kr-1
!!$      rr_Mat(kr,kc) = -1.0d0
!!$      kc = kr
!!$      rr_Mat(kr,kc) =  1.0d0
!!$      !
!!$      r_Vec(kr) = r_Temp(kr-1) &
!!$        & + ( r_Temp(kr+1) - r_Temp(kr) ) &
!!$        &   / ( log( r_Press(kr+1) ) - log( r_Press(kr  ) ) ) &
!!$        &   * ( log( r_Press(kr  ) ) - log( r_Press(kr-1) ) ) &
!!$        & - r_Temp(kr)
!!$      if ( kConvE == 1 ) then
!!$        ! convective layer is assumed to be very thin
!!$        r_Vec(kr) = r_Temp(kr-1) &
!!$          & + ( r_Temp(2) - r_Temp(1) ) &
!!$          &   / log( r_Press(2) / r_Press(1) ) &
!!$          &   * log( r_Press(1) / r_Press(0) ) &
!!$          & - r_Temp(kr)
!!$      else
!!$        ! convective layer is thicker than 1 layer
!!$        r_Vec(kr) = r_Temp(kr-1) + z_VertDelTemp(kr) - r_Temp(kr)
!!$      end if
!!$    end do



    ! T (kConvE+1 <= k)
    do kr = kConvE+1, kmax
      kc = -1
      rr_Mat(kr,kc) = + ( ra_DelRadFlux(kr,kc) - ra_DelRadFlux(kr-1,kc) )
      kc = 0
      rr_Mat(kr,kc) = + ( ra_DelRadFlux(kr,kc) - ra_DelRadFlux(kr-1,kc) )
      kc = 1
      rr_Mat(kr,kc) = + ( ra_DelRadFlux(kr,kc) - ra_DelRadFlux(kr-1,kc) )
      kc = kr
      rr_Mat(kr,kc) = rr_Mat(kr,kc) + Cp * z_DelAtmMass(kc) / DelTime
      !
      r_Vec(kr) = - ( r_RadFlux(kr) - r_RadFlux(kr-1) )
    end do

    call ludecomp_prep_simple ( (kmax+3), rr_Mat )
    call ludecomp_solve_simple( (kmax+3), rr_Mat, r_Vec )
    !
    r_DelTemp    = r_Vec(0:kmax)
    DelSurfTemp  = r_Vec(-1)
    !
    SurfSensFlux = r_Vec(-2)


    if ( FlagMoistConvection ) then
      do k = 0, kmax
        r_PressH2OVapSat(k) = &
          & LOCALPressH2OVapSat( LatentHeatPerMol, r_Temp(k)+r_DelTemp(k) )
      end do
      do k = 0, kConvE
        r_DelVMRH2OVap(k) = RelHumid * r_PressH2OVapSat(k) / r_Press(k) - r_VMRH2OVap(k)
      end do
      do k = kConvE+1, kmax
        VMRH2OVapTentative = r_VMRH2OVap(k-1) + r_DelVMRH2OVap(k-1)
        VMRH2OVapTentative = min( VMRH2OVapTentative, r_PressH2OVapSat(k) / r_Press(k) )
        r_DelVMRH2OVap(k) = VMRH2OVapTentative - r_VMRH2OVap(k)
      end do
    else
      r_DelVMRH2OVap = 0.0d0
    end if


  end subroutine ConvCalcRadConvTendencyCore

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

  subroutine ConvCalcRadConvTendencyCoreTrial( &
    & FlagMoistConvection, &
    & kConvS, kConvE, &
    & Grav, Cp, SurfHeatCap, DelTime, &
    & kmax, r_Press, r_Temp, SurfTemp, r_VMRH2OVap, r_MeanMolWt, &
    & r_RadUwFlux, r_RadDwFlux, ra_DelRadUwFlux, ra_DelRadDwFlux, &
    & r_DelTemp, DelSurfTemp, r_DelVMRH2OVap, SurfSensFlux &
    & )

    logical , intent(in ) :: FlagMoistConvection
    ! Full level index (mid-poing of layers)
    integer , intent(in ) :: kConvS
    integer , intent(in ) :: kConvE
    real(dp), intent(in ) :: Grav
    real(dp), intent(in ) :: Cp
    real(dp), intent(in ) :: SurfHeatCap
    real(dp), intent(in ) :: DelTime
    integer , intent(in ) :: kmax
    real(dp), intent(in ) :: r_Press        (0:kmax)
    real(dp), intent(in ) :: r_Temp         (0:kmax)
    real(dp), intent(in ) :: SurfTemp
    real(dp), intent(in ) :: r_VMRH2OVap    (0:kmax)
    real(DP), intent(in ) :: r_MeanMolWt    (0:kmax)
    real(dp), intent(in ) :: r_RadUwFlux    (0:kmax)
    real(dp), intent(in ) :: r_RadDwFlux    (0:kmax)
    real(DP), intent(in ) :: ra_DelRadUwFlux(0:kmax,-1:1)
    real(DP), intent(in ) :: ra_DelRadDwFlux(0:kmax,-1:1)
    real(dp), intent(out) :: r_DelTemp      (0:kmax)
    real(dp), intent(out) :: DelSurfTemp
    real(DP), intent(out) :: r_DelVMRH2OVap (0:kmax)
    real(DP), intent(out) :: SurfSensFlux


    ! Function
    !
!!$    real(DP) :: LOCALPressH2OVapSat
!!$    real(DP) :: LOCALDPressH2OVapSatDTemp

    ! Local variables
    !
    real(DP), parameter :: LatentHeat = 2.5d6    ! J kg-1
    real(DP), parameter :: CpH2OVap   = 1810.0d0 ! J K-1 kg-1
    real(DP), parameter :: MolWtH2O   = 18.0d-3  ! kg mol-1
    real(dp)            :: LatentHeatPerMol


    real(dp) :: RelHumid = 1.0d0
    real(dp) :: z_DelAtmMass (1:kmax)
    real(dp) :: r_PressH2OVapSat      (0:kmax)
    real(dp) :: r_DPressH2OVapSatDTemp(0:kmax)
    real(dp) :: r_DMMRH2OVapSatDTemp  (0:kmax)
    real(dp) :: z_VertDelTemp(1:kmax)
    real(dp) :: r_RadFlux    (0:kmax)
    real(DP) :: ra_DelRadFlux(0:kmax,-1:1)
    real(dp) :: rr_Mat(-2:kmax,-2:kmax)
    real(dp) :: r_Vec (-2:kmax)
    integer  :: k
    integer  :: kr
    integer  :: kc


    LatentHeatPerMol = LatentHeat * MolWtH2O

    r_RadFlux     = r_RadUwFlux     - r_RadDwFlux
    ra_DelRadFlux = ra_DelRadUwFlux - ra_DelRadDwFlux


    do k = 1, kmax
      z_DelAtmMass(k) = ( r_Press(k-1) - r_Press(k) ) / Grav
    end do

    if ( FlagMoistConvection ) then
      do k = 0, kmax
        r_PressH2OVapSat(k) = &
          & LOCALPressH2OVapSat( LatentHeatPerMol, r_Temp(k) )
      end do
      do k = 0, kmax
        r_DPressH2OVapSatDTemp(k) = &
          & LOCALDPressH2OVapSatDTemp( LatentHeatPerMol, r_Temp(k) )
      end do
    else
      r_PressH2OVapSat = 0.0d0
      r_DPressH2OVapSatDTemp = 0.0d0
    end if
    r_DMMRH2OVapSatDTemp = r_DPressH2OVapSatDTemp / r_Press &
      & * MolWtH2O / r_MeanMolWt


    call ConvCalcMoistConvDelTemp( &
      & FlagMoistConvection, &
      & kmax, r_Press, r_MeanMolWt, r_Temp, &
      & z_VertDelTemp &
      & )


    rr_Mat = 0.0d0
    r_Vec  = 0.0d0
    !
    ! Energy conservation
    kr = -2
    kc = -2
    rr_Mat(kr,kc) = -1.0d0
    kc = -1
    rr_Mat(kr,kc) = &
      & + (   ( ra_DelRadUwFlux(kConvE,kc) - ra_DelRadDwFlux(kConvE+1,kc) ) &
      &     - ( ra_DelRadUwFlux(0     ,kc) - ra_DelRadDwFlux(0       ,kc) ) )
    kc = 0
    rr_Mat(kr,kc) = &
      & + (   ( ra_DelRadUwFlux(kConvE,kc) - ra_DelRadDwFlux(kConvE+1,kc) ) &
      &     - ( ra_DelRadUwFlux(0     ,kc) - ra_DelRadDwFlux(0       ,kc) ) )
    kc = 1
    rr_Mat(kr,kc) = &
      & + (   ( ra_DelRadUwFlux(kConvE,kc) - ra_DelRadDwFlux(kConvE+1,kc) ) &
      &     - ( ra_DelRadUwFlux(0     ,kc) - ra_DelRadDwFlux(0       ,kc) ) )
    !
    kc = 0
    rr_Mat(kr,kc) = rr_Mat(kr,kc) &
      & + Cp / DelTime * z_DelAtmMass(kc+1)/2.0d0
    do kc = 1, kConvE
!!$      rr_Mat(kr,kc) = rr_Mat(kr,kc) &
!!$        & + z_DelAtmMass(kc) / DelTime * ( Cp + RelHumid * LatentHeat * r_DMMRH2OVapSatDTemp(kc) )
      ! Latent heat is neglected to conserve energy
!!$      rr_Mat(kr,kc) = rr_Mat(kr,kc) &
!!$        & + z_DelAtmMass(kc) / DelTime * ( Cp + 0.0d0 )
      !
      rr_Mat(kr,kc) = rr_Mat(kr,kc) &
        & + Cp / DelTime * ( z_DelAtmMass(kc) + z_DelAtmMass(kc+1) ) / 2.0d0
    end do
    !
    r_Vec(kr)     = &
      & - (   ( r_RadUwFlux(kConvE) - r_RadDwFlux(kConvE+1) ) &
      &     - ( r_RadUwFlux(0     ) - r_RadDwFlux(0       ) ) )

    !
    ! Ts
    kr = -1
    !
    kc = -2
    rr_Mat(kr,kc) = 1.0d0
    kc = -1
    rr_Mat(kr,kc) = SurfHeatCap / DelTime + ra_DelRadUwFlux(0,kc)
    kc =  0
    rr_Mat(kr,kc) = + ra_DelRadFlux(0,kc)
    kc =  1
    rr_Mat(kr,kc) = + ra_DelRadFlux(0,kc)
    !
    r_Vec(kr)     = - r_RadFlux(0)
    !
    ! T (k=0)
    !   T0 = Ts
    !   r_Temp(0) + r_DelTemp(0) = SurfTemp + DelSurfTemp
    !   - DelSurfTemp + r_DelTemp(0) = SurfTemp - r_Temp(0)
    kr = 0
    kc =  -1
    rr_Mat(kr,kc) = -1.0d0
    kc =   0
    rr_Mat(kr,kc) =  1.0d0
    !
    r_Vec(kr)     = SurfTemp - r_Temp(0)
    !
    ! T
    do kr = 1, kConvE
      kc = kr-1
      rr_Mat(kr,kc) = -1.0d0
      kc = kr
      rr_Mat(kr,kc) =  1.0d0
      !
      r_Vec(kr) = r_Temp(kr-1) + z_VertDelTemp(kr) - r_Temp(kr)
    end do
    do kr = kConvE+1, kmax-1
      kc = -1
      rr_Mat(kr,kc) = &
        & + (   ( ra_DelRadUwFlux(kr  ,kc) - ra_DelRadDwFlux(kr+1,kc) ) &
        &     - ( ra_DelRadUwFlux(kr-1,kc) - ra_DelRadDwFlux(kr  ,kc) ) )
      kc =  0
      rr_Mat(kr,kc) = &
        & + (   ( ra_DelRadUwFlux(kr  ,kc) - ra_DelRadDwFlux(kr+1,kc) ) &
        &     - ( ra_DelRadUwFlux(kr-1,kc) - ra_DelRadDwFlux(kr  ,kc) ) )
      kc =  1
      rr_Mat(kr,kc) = &
        & + (   ( ra_DelRadUwFlux(kr  ,kc) - ra_DelRadDwFlux(kr+1,kc) ) &
        &     - ( ra_DelRadUwFlux(kr-1,kc) - ra_DelRadDwFlux(kr  ,kc) ) )
      kc = kr
      rr_Mat(kr,kc) = rr_Mat(kr,kc) &
        & + Cp / DelTime * ( z_DelAtmMass(kc) + z_DelAtmMass(kc+1) ) / 2.0d0
      !
      r_Vec(kr)     = &
        & - (   ( r_RadUwFlux(kr  ) - r_RadDwFlux(kr+1) ) &
        &     - ( r_RadUwFlux(kr-1) - r_RadDwFlux(kr  ) ) )
    end do
    kr = kmax
    kc = -1
    rr_Mat(kr,kc) = &
      & + (   ( ra_DelRadUwFlux(kr  ,kc) - 0.0d0 ) &
      &     - ( ra_DelRadUwFlux(kr-1,kc) - 0.0d0 ) )
    kc =  0
    rr_Mat(kr,kc) = &
      & + (   ( ra_DelRadUwFlux(kr  ,kc) - 0.0d0 ) &
      &     - ( ra_DelRadUwFlux(kr-1,kc) - 0.0d0 ) )
    kc =  1
    rr_Mat(kr,kc) = &
      & + (   ( ra_DelRadUwFlux(kr  ,kc) - 0.0d0 ) &
      &     - ( ra_DelRadUwFlux(kr-1,kc) - 0.0d0 ) )
    kc = kr
    rr_Mat(kr,kc) = rr_Mat(kr,kc) &
      & + Cp / DelTime * ( z_DelAtmMass(kc) + 0.0d0 ) / 2.0d0
    !
    r_Vec(kr)     = &
      & - (   ( r_RadUwFlux(kr  ) - 0.0d0 ) &
      &     - ( r_RadUwFlux(kr-1) - 0.0d0 ) )


    call ludecomp_prep_simple ( (kmax+3), rr_Mat(-2:kmax,-2:kmax) )
    call ludecomp_solve_simple( (kmax+3), rr_Mat(-2:kmax,-2:kmax), r_Vec(-2:kmax) )
    !
    DelSurfTemp       = r_Vec(-1)
    r_DelTemp(0:kmax) = r_Vec(0:kmax)
    !
    SurfSensFlux = r_Vec(-2)


    if ( FlagMoistConvection ) then
      do k = 0, kmax
        r_PressH2OVapSat(k) = &
          & LOCALPressH2OVapSat( LatentHeatPerMol, r_Temp(k)+r_DelTemp(k) )
      end do
      do k = 0, kConvE
        r_DelVMRH2OVap(k) = RelHumid * r_PressH2OVapSat(k) / r_Press(k) - r_VMRH2OVap(k)
      end do
      do k = kConvE+1, kmax
!!$        r_DelVMRH2OVap(k) = 0.0d0
        r_DelVMRH2OVap(k) = &
          & min( r_VMRH2OVap(k-1), r_PressH2OVapSat(k) ) - r_VMRH2OVap(k)
      end do
    else
      r_DelVMRH2OVap = 0.0d0
    end if


  end subroutine ConvCalcRadConvTendencyCoreTrial

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

  function LOCALPressH2OVapSat( LatentHeatPerMol, Temp ) result( PressH2OVapSat )

    use constants0, only : GasRUniv

    real(DP), intent(in ) :: LatentHeatPerMol
    real(DP), intent(in ) :: Temp

    real(DP) :: PressH2OVapSat


    PressH2OVapSat = 1.4d11 * exp( - LatentHeatPerMol / ( GasRUniv * Temp ) )

  end function LOCALPressH2OVapSat

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

  subroutine ConvCalcMoistConvDelTemp( &
    & FlagMoistConvection, &
    & kmax, r_Press, r_MeanMolWt, r_Temp, &
    & z_DelTemp &
    & )

    use constants0, only : GasRUniv
    use rad_planetconst, only : Cp

    implicit none

    logical , intent(in ) :: FlagMoistConvection
    integer , intent(in ) :: kmax
    real(DP), intent(in ) :: r_Press    (0:kmax)
    real(DP), intent(in ) :: r_MeanMolWt(0:kmax)
    real(DP), intent(in ) :: r_Temp     (0:kmax)
    real(DP), intent(out) :: z_DelTemp    (1:kmax)


    ! Local variables
    !
    real(DP), parameter :: LatentHeat = 2.5d6    ! J kg-1
    real(DP), parameter :: CpH2OVap   = 1810.0d0 ! J K-1 kg-1
    real(DP), parameter :: MolWtH2O   = 18.0d-3  ! kg mol-1

    real(DP) :: CPPerMol
    real(DP) :: CpH2OVapPerMol
    real(DP) :: CPTotPerMol
    real(DP) :: LatentHeatPerMol
    real(DP) :: PressH2OVapSat
    real(DP) :: DPressH2OVapSatDTemp
    real(DP) :: VMRSat

    real(DP) :: r_PressH2OVapSat      (0:kmax)
    real(DP) :: r_DPressH2OVapSatDTemp(0:kmax)

!!$    real(DP) :: LOCALPressH2OVapSat
!!$    real(DP) :: LOCALDPressH2OVapSatDTemp

    real(DP) :: r_DTempDLogPress(0:kmax)

    integer  :: k


    LatentHeatPerMol = LatentHeat * MolWtH2O
    CpH2OVapPerMol   = CpH2OVap * MolWtH2O

!!$    PressH2OVapSat   = LOCALPressH2OVapSat( LatentHeatPerMol, r_Temp(0) )
!!$    VMRSat           = PressH2OVapSat / r_Press(0)

    if ( FlagMoistConvection ) then
      do k = 0, kmax
        r_PressH2OVapSat      (k) = LOCALPressH2OVapSat      ( LatentHeatPerMol, r_Temp(k) )
        r_DPressH2OVapSatDTemp(k) = LOCALDPressH2OVapSatDTemp( LatentHeatPerMol, r_Temp(k) )
      end do
    else
      r_PressH2OVapSat       = 0.0d0
      r_DPressH2OVapSatDTemp = 0.0d0
    end if

    do k = 0, kmax

      ! moist convection

      CpPerMol = Cp * r_MeanMolWt(k)
      !
      PressH2OVapSat       = r_PressH2OVapSat      (k)
      DPressH2OVapSatDTemp = r_DPressH2OVapSatDTemp(k)
      !
      VMRSat = PressH2OVapSat / r_Press(k)
      CpTotPerMol = CpPerMol * ( 1.0d0 - VMRSat ) + CpH2OVapPerMol * VMRSat

      ! Is this wrong?
!!$      z_DelTemp(k) = ( GasRUniv * r_Temp(k-1) + LatentHeatPerMol * VMRSat ) &
!!$        & / ( CPTotPerMol + LatentHeatPerMol / r_Press(k-1) * DPressH2OVapSatDTemp ) &
!!$        & * log( r_Press(k) / r_Press(k-1) )
      ! Is this correct?
!!$      z_DelTemp(k) = &
!!$        &   ( GasRUniv * r_Temp(k-1) + LatentHeatPerMol * VMRSat / ( 1.0 - VMRSat ) ) &
!!$        & / ( CPPerMol + LatentHeatPerMol / ( r_Press(k-1) * ( 1.0 - VMRSat ) ) * DPressH2OVapSatDTemp ) &
!!$        & * log( r_Press(k) / r_Press(k-1) )
      r_DTempDLogPress(k) = &
        & - ( GasRUniv * r_Temp(k) + LatentHeatPerMol * VMRSat / ( 1.0 - VMRSat ) ) &
        & / ( CPPerMol + LatentHeatPerMol / ( r_Press(k) * ( 1.0 - VMRSat ) ) * DPressH2OVapSatDTemp )

    end do

    do k = 1, kmax
      z_DelTemp(k) = &
        &    ( r_DTempDLogPress(k-1) + r_DTempDLogPress(k) ) / 2.0d0 &
        & * log( r_Press(k-1) / r_Press(k) )
    end do


  end subroutine ConvCalcMoistConvDelTemp

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

  function LOCALDPressH2OVapSatDTemp( LatentHeatPerMol, Temp ) result( DPressH2OVapSatDTemp )

    use constants0, only : GasRUniv

    real(DP), intent(in ) :: LatentHeatPerMol
    real(DP), intent(in ) :: Temp

    real(DP) :: DPressH2OVapSatDTemp


    DPressH2OVapSatDTemp = &
      & 1.4d11 * exp( - LatentHeatPerMol / ( GasRUniv * Temp ) ) &
      & * LatentHeatPerMol / ( GasRUniv * Temp**2 )

  end function LOCALDPressH2OVapSatDTemp

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

end module convection
