scale_atm_phy_cp_dgm_mconv_adjustment.F90

Go to the documentation of this file.

!> module FElib / Atmosphere / Physics cumulus parameterization
!!
!! @par Description
!! A module providing a moist convective adjustment scheme.
!!
!! Saturated convective instability is diagnosed from the vertical gradient of saturated moist static energy in nearly saturated layers.
!! The diagnosed layer is adjusted toward a moist-neutral temperature profile subject to column-integrated moist-static-energy and available-water constraints. 
!! Condensed water is immediately removed as precipitation.
!!
!! @author Yuta Kawai, Team SCALE
!!
!! @par Reference
!!      Manabe, S., J. Smagorinsky, and R. F. Strickler (1965):
!!      Simulated Climatology of a General Circulation Model with a
!!      Hydrologic Cycle. Monthly Weather Review, 93, 769-798.
!!
!! 
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_phy_cp_dgm_mconv_adjustment
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
  use scale_const, only: &
    grav => const_grav,   &
    cpdry => const_cpdry, &
    rdry => const_rdry,   &
    rvap => const_rvap,   &
    pres0 => const_pre00, &
    eps => const_eps
  use scale_atmos_hydrometeor, only: &
    cp_vapor, &
    lhv
  use scale_atmos_saturation, only: &
    atmos_saturation_psat_liq,     &
    atmos_saturation_pres2qsat_liq
 
  use scale_element_base, only: &
    elementbase1d, elementbase2d, elementbase3d
  use scale_localmesh_2d, only: localmesh2d  
  use scale_localmesh_3d, only: localmesh3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private  
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
  public :: atm_phy_cp_dgm_mconv_adjustment_setup
  public :: atm_phy_cp_dgm_mconv_adjustment_calc_tendency
  public :: atm_phy_cp_dgm_mconv_adjustment_finalize
 
  !-----------------------------------------------------------------------------
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  
  !* MCA paramters
 
  ! Iteration limits
  integer, parameter :: MCA_MAX_ADJUST_ITER = 200
  integer, parameter :: MCA_MAX_ENERGY_ITER = 50
  integer, parameter :: MCA_MAX_WATER_ITER  = 60
 
  ! Column-energy root solve
  real(RP), parameter :: MCA_COLUMN_ENERGY_RTOL   = 5.0e-7_rp
  real(RP), parameter :: MCA_COLUMN_ENERGY_ATOL   = 1.0_rp
  real(RP), parameter :: MCA_BOUNDARY_ENERGY_RTOL = 1.0e-6_rp
 
  ! Local saturated-MSE solve for moist-adiabat integration
  real(RP), parameter :: MCA_LOCAL_HMSE_RTOL = 1.0e-9_rp
  real(RP), parameter :: MCA_LOCAL_HMSE_ATOL = 1.0e-4_rp
  real(RP), parameter :: MCA_TEMP_TOL        = 1.0e-5_rp
 
  ! Water-feasibility constraint
  real(RP), parameter :: MCA_WATER_RTOL = 1.0e-9_rp
  real(RP), parameter :: MCA_WATER_ATOL = 1.0e-10_rp
 
  ! Saturation and instability thresholds
  real(RP), parameter :: MCA_RH_FORCED_SATURATION = 0.999_rp
  real(RP), parameter :: MCA_RH_TRIGGER           = 0.999_rp
 
  ! Thresholds for suppressing numerical-noise detection
  real(RP), parameter :: MCA_HMSE_GRAD_TOL      = 1.0e-3_rp
  real(RP), parameter :: MCA_HMSE_DIFF_TOL      = 1.0_rp
  real(RP), parameter :: MCA_MIN_UNSTABLE_DEPTH = 0.0_rp
  real(RP), parameter :: MCA_Z_TOL              = 100.0_rp * epsilon(1.0_rp)
 
  ! Moist-neutral-profile solver status
  integer, parameter :: MCA_PROFILE_SUCCESS                 = 0
  integer, parameter :: MCA_PROFILE_NO_ENERGY_ROOT          = 1
  integer, parameter :: MCA_PROFILE_ADIABAT_FAILURE         = 2
  integer, parameter :: MCA_PROFILE_ENERGY_MAXITER          = 3
  integer, parameter :: MCA_PROFILE_NO_ADJUSTED_NODE        = 4
  integer, parameter :: MCA_PROFILE_NO_WATER_FEASIBLE_STATE = 5
 
contains
  !> Setup a module for moist convective adjustment scheme
  subroutine atm_phy_cp_dgm_mconv_adjustment_setup()
    implicit none
    !----------------------------------------------------
    return
  end subroutine atm_phy_cp_dgm_mconv_adjustment_setup
 
  !> Calculate tendencies with moist convective adjustment scheme
!OCL SERIAL
  subroutine atm_phy_cp_dgm_mconv_adjustment_calc_tendency( &
    DENS_t, RHOT_t, RHOQV_t, SFLX_RAIN, SFLX_ENGI,     & ! (out)
    ddens, drhot, qv, pt, pres, dens_hyd, rtot, cptot, & ! (in)
    dtsec, lmesh, elem, elem1d                         ) ! (in)
 
    use scale_tracer, only: &
      tracer_cv    
    use scale_atmos_hydrometeor, only: &
      cv_water, i_qv    
    implicit none
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    class(elementbase1d), intent(in) :: elem1d
    real(rp), intent(out) :: dens_t(elem%np,lmesh%nea)
    real(rp), intent(out) :: rhot_t(elem%np,lmesh%nea)
    real(rp), intent(out) :: rhoqv_t(elem%np,lmesh%nea)
    real(rp), intent(out) :: sflx_rain(elem%nnode_h1d**2,lmesh%ne2da)
    real(rp), intent(out) :: sflx_engi(elem%nnode_h1d**2,lmesh%ne2da)
    real(rp), intent(in) :: ddens(elem%np,lmesh%nea)
    real(rp), intent(in) :: drhot(elem%np,lmesh%nea)
    real(rp), intent(in) :: qv(elem%np,lmesh%nea)
    real(rp), intent(in) :: pt(elem%np,lmesh%nea)
    real(rp), intent(in) :: pres(elem%np,lmesh%nea)
    real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
    real(rp), intent(in) :: rtot(elem%np,lmesh%nea)
    real(rp), intent(in) :: cptot(elem%np,lmesh%nea)
    real(rp), intent(in) :: dtsec
 
    integer :: ke, ke_xy, ke_z
    integer :: ph, pz, p
    real(rp) :: dens_z(elem%nnode_v,lmesh%nez)
    real(rp) :: pres_z(elem%nnode_v,lmesh%nez)
    real(rp) :: zlev_z(elem%nnode_v,lmesh%nez)
    real(rp) :: temp_z(elem%nnode_v,lmesh%nez)
    real(rp) :: pott_z(elem%nnode_v,lmesh%nez) ! work array
    real(rp) :: qvap_z(elem%nnode_v,lmesh%nez) ! work array
    real(rp) :: qsat
    real(rp) :: rh_z
    real(rp) :: rtot_, cptot_, qdry
 
    real(rp) :: elem_width_z
    real(rp) :: int_weight(elem%nnode_v,lmesh%nez)
 
    real(rp) :: pott_ini(elem%nnode_v,lmesh%nez)
    real(rp) :: qvap_ini(elem%nnode_v,lmesh%nez)
    real(rp) :: dens_ini(elem%nnode_v,lmesh%nez)
    real(rp) :: rhoqvap_ini(elem%nnode_v,lmesh%nez)
 
    integer :: iter_expand
    integer :: mca_max_mask_expand
    integer :: lbase, ltop, lbase_try, ltop_try
    real(rp) :: zbase_try, ztop_try
    logical :: mask_expanded
 
    logical :: is_unstable
    logical  :: unstable_core_mask(elem%nnode_v,lmesh%nez)    
    logical :: forced_saturation_mask(elem%nnode_v,lmesh%nez)
    logical :: adjustment_mask(elem%nnode_v,lmesh%nez)
 
    logical :: current_saturation_mask(elem%nnode_v,lmesh%nez)
    logical :: persistent_saturation_mask(elem%nnode_v,lmesh%nez)
 
    logical :: physical_active_mask(elem%nnode_v,lmesh%nez)
    integer :: iter_adj
    real(rp) :: temp_adj(elem%nnode_v,lmesh%nez)
    real(rp) :: qvap_adj(elem%nnode_v,lmesh%nez)
    real(rp) :: rhoqvap_adj(elem%nnode_v,lmesh%nez)
    real(rp) :: rhoprecip_adj(elem%nnode_v,lmesh%nez)    ! Water-vapor mass density diagnosed to be removed as precipitation during the current adjustment iteration [kg m-3]
    real(rp) :: rhoprecip_accum(elem%nnode_v,lmesh%nez)  ! Accumulated vapor mass density removed as precipitation [kg m-3]
    real(rp) :: precip_engi_accum 
 
    real(rp), allocatable :: zlev_diag(:)
    integer :: nlev_diag
 
    logical :: profile_converged
    integer :: profile_status
    logical :: do_adjustment
    logical :: adjustment_converged
 
    integer :: hslice_b(elem%nnode_h1d**2), hslice_t(elem%nnode_h1d**2)
 
    logical :: debug_flag
 
    integer  :: lactive_base
    integer  :: lactive_top
    logical  :: active_region_initialized    
    !----------------------------------------------------
 
    nlev_diag = lmesh%NeZ * (elem%Nnode_v-1) + 1
    allocate( zlev_diag(nlev_diag) )
    
    hslice_b(:) = elem%Hslice(:,1)
    hslice_t(:) = elem%Hslice(:,elem%Nnode_v)
 
    mca_max_mask_expand = nlev_diag
 
    !$omp parallel do collapse(2) private(ke,p, &
    !$omp dens_z, pres_z, temp_z, zlev_z, pott_z, qvap_z, qsat, rh_z, rtot_,cptot_,                 &
    !$omp dens_ini, pott_ini, qvap_ini, rhoqvap_ini,                                                &
    !$omp zlev_diag,                                                                                &
    !$omp lbase, ltop, lbase_try, ltop_try, zbase_try, ztop_try, iter_expand, mask_expanded,        &
    !$omp is_unstable, unstable_core_mask, forced_saturation_mask, adjustment_mask,                 &
    !$omp current_saturation_mask, persistent_saturation_mask, physical_active_mask,                &
    !$omp temp_adj, qdry, qvap_adj, rhoqvap_adj, rhoprecip_adj, rhoprecip_accum,                    &
    !$omp profile_converged, profile_status, do_adjustment, adjustment_converged,                   &
    !$omp int_weight,elem_width_z, debug_flag, lactive_base, lactive_top, active_region_initialized )
    do ke_xy=1, lmesh%Ne2D
    do ph=1, elem%Nnode_h1D**2
    
      ! if ( lmesh%PRC_myrank == 1 .and. ph == 32 .and. ke_xy == 13 ) then
      !   debug_flag = .true.
      ! else
        debug_flag = .false.
      ! end if
 
      !- Extract vertical 1D DG column
 
      do ke_z=1, lmesh%NeZ
        ke = ke_xy + (ke_z-1)*lmesh%Ne2D
 
        elem_width_z = lmesh%zlev(hslice_t(ph),ke) - lmesh%zlev(hslice_b(ph),ke)
        do pz=1, elem%Nnode_v
          p = ph + (pz-1)*elem%Nnode_h1D**2
          dens_z(pz,ke_z) = ddens(p,ke) + dens_hyd(p,ke)
          pres_z(pz,ke_z) = pres(p,ke)
          zlev_z(pz,ke_z) = lmesh%zlev(p,ke)
          int_weight(pz,ke_z) = elem1d%IntWeight_lgl(pz) * 0.5_rp * elem_width_z
 
          temp_z(pz,ke_z) = pres(p,ke)/ ( dens_z(pz,ke_z) * rtot(p,ke) )
          qvap_z(pz,ke_z) = qv(p,ke)
          pott_z(pz,ke_z) = pt(p,ke)
        end do
      end do
 
      !- Set initial state
 
      do ke_z=1, lmesh%NeZ
        dens_ini(:,ke_z) = dens_z(:,ke_z)
        pott_ini(:,ke_z) = pott_z(:,ke_z)
        qvap_ini(:,ke_z) = qvap_z(:,ke_z)
        rhoqvap_ini(:,ke_z) = dens_ini(:,ke_z) * qvap_ini(:,ke_z)
 
        temp_adj(:,ke_z) = temp_z(:,ke_z)
        rhoqvap_adj(:,ke_z) = rhoqvap_ini(:,ke_z)
        rhoprecip_accum(:,ke_z) = 0.0_rp
      end do
 
      !* Iteratively diagnose and remove remaining saturated convective instability.  ***************************************
      !  A single profile adjustment may generate a new unstable layer adjacent to or outside the previously adjusted region.
 
      adjustment_converged = .false.
      do_adjustment        = .false.
 
      active_region_initialized = .false.
      lactive_base = 0
      lactive_top  = 0
 
      forced_saturation_mask(:,:)     = .false.
      current_saturation_mask(:,:)    = .false.
      persistent_saturation_mask(:,:) = .false.
      physical_active_mask(:,:)       = .false.
 
      do iter_adj=1, mca_max_adjust_iter
 
        qvap_adj(:,:) = rhoqvap_adj(:,:) / dens_z(:,:)
 
        ! Reconstruct the saturation constraint from the current state.
        current_saturation_mask(:,:) = .false.
        forced_saturation_mask(:,:)  = .false.        
 
        !- Diagnose unstable layers 
 
        call diagnose_convective_layer( lbase, ltop, is_unstable, zlev_diag,     & ! (out)
          temp_adj, qvap_adj, pres_z, zlev_z, lmesh, elem, nlev_diag, debug_flag ) ! (in)
 
        if ( debug_flag ) then
          write(*,*) "ke_xy=", ke_xy, "ph=", ph, "iter_adj=", iter_adj
          write(*,*) "   lbase, ltop, is_unstable = ", lbase, ltop, is_unstable
        end if
 
        if ( .not. is_unstable ) then
          adjustment_converged = .true.
          exit
        end if
        do_adjustment = .true.
 
        ! Maintain the smallest diagnostic-level interval containing
        ! all physically diagnosed unstable regions belonging to the
        ! current connected convective-adjustment event.        
        if ( .not. active_region_initialized ) then
          lactive_base = lbase
          lactive_top  = ltop
          active_region_initialized = .true.
        else if ( lbase <= lactive_top + 1 .and. &
                  ltop  >= lactive_base - 1 ) then
          ! Overlapping or directly touching region:
          ! continue the same connected adjustment.
          lactive_base = min(lactive_base, lbase)
          lactive_top  = max(lactive_top,  ltop)
        else
          ! Disconnected instability:
          ! start a new independent adjustment region.
          lactive_base = lbase
          lactive_top  = ltop
          persistent_saturation_mask(:,:) = .false.
        end if
 
        if ( debug_flag ) then
          write(*,*) "   persistent active region=", lactive_base, lactive_top
        end if        
 
        ! Construct the physically connected active region.
        ! This does not include numerical expansion used only for water/energy feasibility.
        call make_dg_vertical_range_mask( &
          physical_active_mask,                                    & ! (out)
          zlev_z, zlev_diag(lactive_base), zlev_diag(lactive_top), & ! (in)
          elem%Nnode_v, lmesh%NeZ                                  ) ! (in)
        
        ! Mark nearly saturated DG nodes inside the diagnosed unstable core. 
        ! These nodes will be constrained to saturation in each trial moist-neutral profile.
        call make_dg_vertical_range_mask( &
          unstable_core_mask,                        & ! (out)
          zlev_z, zlev_diag(lbase), zlev_diag(ltop), & ! (in)
          elem%Nnode_v, lmesh%NeZ                    ) ! (in)        
 
        do ke_z = 1, lmesh%NeZ
        do pz = 1, elem%Nnode_v
          if ( unstable_core_mask(pz,ke_z) ) then
            call atmos_saturation_pres2qsat_liq( temp_adj(pz,ke_z), pres_z(pz,ke_z),  &
              qsat ) ! (out)
 
            rh_z = qvap_adj(pz,ke_z) / max(qsat, eps)
            current_saturation_mask(pz,ke_z) = rh_z >= mca_rh_forced_saturation
          end if
        end do
        end do
 
        ! Keep saturation constraints inherited from previous adjustment iterations within this connected convective event.
 
        forced_saturation_mask(:,:) = persistent_saturation_mask(:,:) .or. current_saturation_mask(:,:)
 
        !- Construct a moist-neutral profile.
        !  If no water-feasible and energy-conserving solution exists, expand the adjustment layer and retry.        
 
        lbase_try = lactive_base
        ltop_try  = lactive_top
 
        profile_converged     = .false.
        profile_status        = mca_profile_energy_maxiter
 
        do iter_expand=0, mca_max_mask_expand
          zbase_try = zlev_diag(lbase_try)
          ztop_try  = zlev_diag(ltop_try)
          
          !- Construct a mask for the adjustment layer
          call make_dg_vertical_range_mask( adjustment_mask,     & ! (out)
            zlev_z, zbase_try, ztop_try, elem%Nnode_v, lmesh%NeZ ) ! (in)
 
          if ( debug_flag ) then
            write(*,*) "* Adjustment-profile attempt: iter_expand=", iter_expand
            write(*,*) "   lactive_base, lactive_top=", lactive_base, lactive_top
            write(*,*) "   lbase_try, ltop_try=", lbase_try, ltop_try
            write(*,*) "   zbase_try, ztop_try=", zbase_try, ztop_try
            write(*,*) "   saturation, mixing nodes=", count(forced_saturation_mask), count(adjustment_mask)
          end if
 
          !- Solve for a water-feasible and energy-conserving trial moist-neutral profile
 
          call solve_moist_neutral_adjustment( &
            temp_adj, qvap_adj, rhoqvap_adj, rhoprecip_adj, profile_converged, profile_status,       & ! (out)
            qvap_z, pres_z, zlev_z, dens_z, temp_z,                                                  & ! (in)
            int_weight, adjustment_mask, forced_saturation_mask, elem%Nnode_v, lmesh%NeZ, debug_flag ) ! (in)
 
          if ( profile_converged ) then
            ! Update the persistent saturation history using the post-adjustment state in the physically connected region.
            do ke_z = 1, lmesh%NeZ
            do pz = 1, elem%Nnode_v
              if ( physical_active_mask(pz,ke_z) ) then
              
                call atmos_saturation_pres2qsat_liq( temp_adj(pz,ke_z), pres_z(pz,ke_z), &
                  qsat )
            
                rh_z = qvap_adj(pz,ke_z) / qsat
                if ( rh_z >= mca_rh_forced_saturation ) then
                  persistent_saturation_mask(pz,ke_z) = .true.
                end if
              end if
            end do
            end do
 
            if ( debug_flag ) then
              write(*,*) "Moist-neutral profile converged: iter_expand=", iter_expand
              write(*,*) "final solve region: lbase, ltop=", lbase_try, ltop_try
              write(*,*) "physical active region=", lactive_base, lactive_top
              write(*,*) "persistent saturated nodes=", count(persistent_saturation_mask)
            end if
            exit
          end if
 
          select case ( profile_status )
          case ( mca_profile_no_energy_root )
            call expand_adjustment_layer( lbase_try, ltop_try, mask_expanded, nlev_diag )
            if ( .not. mask_expanded ) then
              exit
            end if
          case default
            exit
          end select
 
        end do ! end loop for iter_expand
 
        if ( .not. profile_converged ) then
          log_info("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "iter_adj=", iter_adj, "ph=", ph, "ke_xy=", ke_xy, "iter_expand=", iter_expand
          select case ( profile_status )
          case ( mca_profile_no_energy_root )
            log_info("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "No energy-conserving root exists in the water-feasible temperature interval, even after adjustment-layer expansion."
            log_info("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "original lbase,ltop=", lbase, ltop, "final lbase,ltop=", lbase_try, ltop_try
            do_adjustment = .false.
            exit
            ! call PRC_abort
          case ( mca_profile_no_water_feasible_state )
            log_info("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "No water-feasible trial moist-neutral profile exists."
            call prc_abort
          case ( mca_profile_adiabat_failure )
            log_info("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "Moist-adiabat profile integration failed."
            call prc_abort
          case ( mca_profile_energy_maxiter )
            log_info("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "Energy root was bracketed, but bisection did not converge."
            call prc_abort
          case ( mca_profile_no_adjusted_node )
            log_info("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "No adjusted DG node was found."
            call prc_abort
          case default
            log_info("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "Unknown moist-neutral-profile construction error."
            call prc_abort
          end select
        end if
      
        ! Accumulate the precipitation density from the current iteration
        rhoprecip_accum(:,:) = rhoprecip_accum(:,:) + rhoprecip_adj(:,:)
 
        ! Update the state for the next iteration
        temp_z(:,:) = temp_adj(:,:)
        qvap_z(:,:) = qvap_adj(:,:)
      
        do ke_z = 1, lmesh%NeZ
        do pz   = 1, elem%Nnode_v
          qdry = 1.0_rp - qvap_z(pz,ke_z)
          rtot_ = rdry * qdry &
                + rvap * qvap_z(pz,ke_z)
          pres_z(pz,ke_z) = dens_z(pz,ke_z) * rtot_  * temp_z(pz,ke_z)
        enddo
        enddo        
      end do ! End loop for iteration
  
      if ( do_adjustment .and. ( .not. adjustment_converged ) ) then
        log_info("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "Moist convective adjustment did not converge: ph=", ph, "ke_xy=", ke_xy
        call prc_abort
      end if
      
      
      !- Convert the adjusted state to DG tendencies
 
      sflx_rain(ph,ke_xy) = 0.0_rp
      sflx_engi(ph,ke_xy) = 0.0_rp
 
      if ( do_adjustment ) then
 
        do ke_z=1, lmesh%NeZ
          ke = ke_xy + (ke_z-1)*lmesh%Ne2D
 
          dens_z(:,ke_z) = dens_ini(:,ke_z) - rhoprecip_accum(:,ke_z)
          qvap_z(:,ke_z) = rhoqvap_adj(:,ke_z) / dens_z(:,ke_z)
 
          do pz=1, elem%Nnode_v
            qdry = 1.0_rp - qvap_z(pz,ke_z)
            rtot_ =  rdry  * qdry +     rvap * qvap_z(pz,ke_z)
            cptot_ = cpdry * qdry + cp_vapor * qvap_z(pz,ke_z)
 
            pres_z(pz,ke_z) = dens_z(pz,ke_z) * rtot_ * temp_z(pz,ke_z)
            pott_z(pz,ke_z) = temp_adj(pz,ke_z) * ( pres0 / pres_z(pz,ke_z) )**( rtot_ / cptot_ )
 
            sflx_rain(ph,ke_xy) = sflx_rain(ph,ke_xy) &
                                - rhoprecip_accum(pz,ke_z) * int_weight(pz,ke_z) / dtsec
            sflx_engi(ph,ke_xy) = sflx_engi(ph,ke_xy) &
                                - temp_adj(pz,ke_z) * rhoprecip_accum(pz,ke_z) * int_weight(pz,ke_z) * cptot_ / dtsec
          end do
 
          do pz=1, elem%Nnode_v
            p = ph + (pz-1)*elem%Nnode_h1D**2
            dens_t(p,ke) = ( dens_z(pz,ke_z) - dens_ini(pz,ke_z) ) / dtsec
            rhot_t(p,ke) = ( dens_z(pz,ke_z) * pott_z(pz,ke_z) - dens_ini(pz,ke_z) * pott_ini(pz,ke_z) ) / dtsec
            rhoqv_t(p,ke) = ( dens_z(pz,ke_z) * qvap_z(pz,ke_z) - rhoqvap_ini(pz,ke_z) ) / dtsec
          end do
        end do
 
        if ( debug_flag ) then
          write(*,*) " MCA summary:"
          write(*,*) "   del_DENS: ", dens_z(:,:) - dens_ini(:,:)
          write(*,*) "   del_RHOT: ", dens_z(:,:) * pott_z(:,:) - dens_ini(:,:) * pott_ini(:,:)
          write(*,*) "   del_RHOQV: ", dens_z(:,:) * qvap_z(:,:) - rhoqvap_ini(:,:)
          write(*,*) "----------------------------------------------------------------"
        end if
      else
        do ke_z=1, lmesh%NeZ
          ke = ke_xy + (ke_z-1)*lmesh%Ne2D
          do pz=1, elem%Nnode_v
            p = ph + (pz-1)*elem%Nnode_h1D**2
            dens_t(p,ke) = 0.0_rp
            rhot_t(p,ke) = 0.0_rp
            rhoqv_t(p,ke) = 0.0_rp
          end do
        end do
      end if
 
    end do ! end loop for ph
    end do ! end loop for ke_xy
 
    return
  end subroutine atm_phy_cp_dgm_mconv_adjustment_calc_tendency
 
  !> Finalize a module for moist convective adjustment scheme
  subroutine atm_phy_cp_dgm_mconv_adjustment_finalize()
    implicit none
    !----------------------------------------------------
    return
  end subroutine atm_phy_cp_dgm_mconv_adjustment_finalize
 
 
!- private subroutines --------------------------
 
!OCL SERIAL
  subroutine build_diagnostic_profile_from_dg( &
    hmse_diag, rh_diag, zlev_diag,             & ! (out)
    hmse_z, rh_z, zlev_z, lmesh, elem, nlev    ) ! (in)
    implicit none
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: nlev
    real(rp), intent(out) :: hmse_diag(nlev)
    real(rp), intent(out) :: rh_diag(nlev)
    real(rp), intent(out) :: zlev_diag(nlev)
    real(rp), intent(in) :: hmse_z(elem%nnode_v,lmesh%nez)
    real(rp), intent(in) :: rh_z(elem%nnode_v,lmesh%nez)
    real(rp), intent(in) :: zlev_z(elem%nnode_v,lmesh%nez)
 
    integer :: pz, ke_z
    integer :: l
    integer :: nnode_v
    !----------------------------------------------------
 
    nnode_v = elem%Nnode_v
    l = 1
 
    hmse_diag(l) = hmse_z(1,1)
    rh_diag(l) = rh_z(1,1)
    zlev_diag(l) = zlev_z(1,1)
    do ke_z=1, lmesh%NeZ
      do pz=2, nnode_v-1
        l = l + 1
        hmse_diag(l) = hmse_z(pz,ke_z)
        rh_diag(l) = rh_z(pz,ke_z)
        zlev_diag(l) = zlev_z(pz,ke_z)
      end do
      l = l + 1
      if ( ke_z < lmesh%NeZ ) then
        hmse_diag(l) = 0.5_rp * ( hmse_z(nnode_v,ke_z) + hmse_z(1,ke_z+1) )
        rh_diag(l) = 0.5_rp * ( rh_z(nnode_v,ke_z) + rh_z(1,ke_z+1) )
        zlev_diag(l) = 0.5_rp * ( zlev_z(nnode_v,ke_z) + zlev_z(1,ke_z+1) )
      else
        hmse_diag(l) = hmse_z(nnode_v,ke_z)
        rh_diag(l) = rh_z(nnode_v,ke_z)
        zlev_diag(l) = zlev_z(nnode_v,ke_z)
      end if
    end do
    return
  end subroutine build_diagnostic_profile_from_dg
 
!OCL SERIAL
  subroutine diagnose_convective_layer( lbase, ltop, is_unstable, zlev_diag, & ! (out)
    temp_z, qvap_z, pres_z, zlev_z, lmesh, elem, nlev,  debug_flag )           ! (in)
    implicit none
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: nlev
    integer, intent(out) :: lbase
    integer, intent(out) :: ltop
    logical, intent(out) :: is_unstable
    real(rp), intent(out) :: zlev_diag(nlev)
    real(rp), intent(in) :: temp_z(elem%nnode_v,lmesh%nez)
    real(rp), intent(in) :: qvap_z(elem%nnode_v,lmesh%nez)
    real(rp), intent(in) :: pres_z(elem%nnode_v,lmesh%nez)
    real(rp), intent(in) :: zlev_z(elem%nnode_v,lmesh%nez)
    logical, intent(in) :: debug_flag
    
    integer :: pz, ke_z
    integer :: l
 
    real(rp) :: hmse_sat_z(elem%nnode_v,lmesh%nez)
    real(rp) :: rh_z(elem%nnode_v,lmesh%nez)
    real(rp) :: qsat, cptot
 
    real(rp) :: rh(nlev)
    real(rp) :: hmse_sat(nlev)
 
    real(rp) :: dz
    real(rp) :: dhmse_sat, dhmse_sat_dz
    real(rp) :: rh_lyr
    logical :: unstable_pair
    logical :: is_inside
 
    integer :: nnode_v
    !------------------------------------
 
    nnode_v = elem%Nnode_v
 
    do ke_z=1, lmesh%NeZ
    do pz=1, nnode_v
 
      call atmos_saturation_pres2qsat_liq( temp_z(pz,ke_z), pres_z(pz,ke_z), &
        qsat ) ! (out)
 
      rh_z(pz,ke_z) = qvap_z(pz,ke_z) / max(qsat, eps)
 
      cptot = cpdry * (1.0_rp - qsat) + cp_vapor * qsat
      hmse_sat_z(pz,ke_z) = cptot * temp_z(pz,ke_z) + grav * zlev_z(pz,ke_z) + lhv * qsat
    end do
    end do
 
    call build_diagnostic_profile_from_dg( hmse_sat, rh, zlev_diag, &
      hmse_sat_z, rh_z, zlev_z, &
      lmesh, elem, nlev )
 
    !- Diagnose unstable layers based on the saturated MSE profile and relative humidity
 
    is_unstable = .false.
    is_inside   = .false.
    lbase = 0; ltop = 0
 
    do l=1, nlev-1
      dz = zlev_diag(l+1) - zlev_diag(l)
      dhmse_sat = hmse_sat(l+1) - hmse_sat(l)
      dhmse_sat_dz = dhmse_sat / dz
      rh_lyr = 0.5_rp * ( rh(l) + rh(l+1) )
      ! write(*,*) "l=", l, "zlev(l)=", zlev(l), "zlev(l+1)=", zlev(l+1), &
      !            "rh_lyr=", rh_lyr, "dhmse_sat_dz=", dhmse_sat_dz
      unstable_pair =                                           &
            ( min(rh(l), rh(l+1)) >= mca_rh_trigger )           &
        .and. ( dhmse_sat < -mca_hmse_diff_tol )                &
        .and. ( dhmse_sat_dz < -mca_hmse_grad_tol )
 
 
      if ( unstable_pair ) then
        if ( .not. is_inside ) then
          lbase = l
          is_inside = .true.
        end if
        ltop = l + 1
      else if ( is_inside ) then
       exit
      end if
    end do
 
    if ( is_inside ) then
      if ( zlev_diag(ltop) - zlev_diag(lbase) >= mca_min_unstable_depth ) is_unstable = .true.
    end if
 
    if ( debug_flag .and. is_unstable ) then
      write(*,*) "Diagnose convective layer:"
      write(*,*) "hmse_sat=", hmse_sat(lbase:ltop)
      write(*,*) "rh=", rh(lbase:ltop)
      write(*,*) "hmse_sat_dg=", hmse_sat_z(:,:)
      write(*,*) "rh_dg=", rh_z(:,:)
    end if
 
    return
  end subroutine diagnose_convective_layer
 
!OCL SERIAL
  subroutine make_dg_vertical_range_mask( adjustment_mask, & ! (out)
    zlev, zbase, ztop, npz, nez                            ) ! (in)
    implicit none
    integer, intent(in) :: npz
    integer, intent(in) :: nez
    logical, intent(out) :: adjustment_mask(npz,nez)
    real(rp), intent(in) :: zlev(npz,nez)
    real(rp), intent(in) :: zbase
    real(rp), intent(in) :: ztop
 
    integer :: ke_z, pz
    !----------------------------------------------------
 
    adjustment_mask(:,:) = .false.
    do ke_z = 1, nez
    do pz = 1, npz
      if ( zlev(pz,ke_z) >= zbase - mca_z_tol .and. &
           zlev(pz,ke_z) <= ztop  + mca_z_tol ) then
        adjustment_mask(pz,ke_z) = .true.
      end if
    end do
    end do    
    return
  end subroutine make_dg_vertical_range_mask
 
  !> Expand the adjustment layer by one diagnostic level above and below, if possible
  !! Numerical fallback: enlarge the adjustment region by one diagnostic level on each available side. This is not a physically based entrainment closure.  
!OCL SERIAL
  subroutine expand_adjustment_layer( &
    lbase, ltop, expanded,            & ! (inout,out)
    nlev       ) ! (in)
    implicit none
 
    integer, intent(inout) :: lbase
    integer, intent(inout) :: ltop
    logical, intent(out) :: expanded
    integer, intent(in) :: nlev
 
    logical :: can_expand_below
    logical :: can_expand_above
    !------------------------------------------------------------
 
    can_expand_below = lbase > 1
    can_expand_above = ltop  < nlev
 
    expanded = .false.
 
    ! Expand by one diagnostic level on each side where possible.
 
    if ( can_expand_below ) then
      lbase = lbase - 1
      expanded = .true.
    end if
 
    if ( can_expand_above ) then
      ltop = ltop + 1
      expanded = .true.
    end if
    return
  end subroutine expand_adjustment_layer
 
  !> Solve for a trial moist-neutral profile satisfying: 1. available-water constraint, 2. column-integrated moist-static-energy conservation.
  !!
  !! The pressure and density profiles are held fixed during the solve.
  !! Nodes selected by forced_saturation_mask are set to saturation.
  !! At other adjusted nodes, the original vapor mixing ratio is retained unless it exceeds saturation.  
  !!
!OCL SERIAL
  subroutine solve_moist_neutral_adjustment( &
    temp_adj, qvap_adj, rhoqvap_adj, rhoprecip_adj, converged, profile_status, & ! (out)
    qvap, pres, zlev, dens, temp,                                              & ! (in)
    int_weight, mix_mask, saturation_mask, npz, nez, debug_flag                ) ! (in)
 
    implicit none
    integer, intent(in) :: npz
    integer, intent(in) :: nez
    real(rp), intent(out) :: temp_adj(npz,nez)
    real(rp), intent(out) :: qvap_adj(npz,nez)
    real(rp), intent(out) :: rhoqvap_adj(npz,nez)
    real(rp), intent(out) :: rhoprecip_adj(npz,nez)
    logical, intent(out) :: converged
    integer, intent(out) :: profile_status
    real(rp), intent(in) :: qvap(npz,nez)
    real(rp), intent(in) :: pres(npz,nez)
    real(rp), intent(in) :: zlev(npz,nez)
    real(rp), intent(in) :: dens(npz,nez)
    real(rp), intent(in) :: temp(npz,nez)
    real(rp), intent(in) :: int_weight(npz,nez)
    logical, intent(in) :: mix_mask(npz,nez)
    logical, intent(in) :: saturation_mask(npz,nez)
    logical, intent(in) :: debug_flag
 
    ! Initial column-integrated quantities
    real(rp) :: energy_target
    real(rp) :: water_mass_ini
 
    ! Trial-profile quantities
    real(rp) :: energy_trial
    real(rp) :: water_trial
    real(rp) :: residual
 
    ! Energy root bracket
    real(rp) :: tbase_lo, tbase_mid, tbase_hi
    real(rp) :: residual_lo, residual_hi
 
    real(rp) :: water_tmin, energy_tmin    
    real(rp) :: water_wlo, water_whi
    real(rp) :: energy_wlo, energy_whi
 
    real(rp) :: energy_lo, energy_hi
    real(rp) :: water_lo, water_hi
 
    ! Water-feasible upper-bound search
    real(rp) :: twater_lo, twater_mid, twater_hi
    real(rp) :: water_mid, energy_mid
 
    ! Work arrays
    real(rp) :: temp_trial(npz,nez)
    real(rp) :: pres_trial(npz,nez)
    real(rp) :: qvap_trial(npz,nez)
 
    ! Precipitation diagnostics
    real(rp) :: water_mass_adj
    real(rp) :: precip_mass_local
    real(rp) :: precip_mass_total
    real(rp) :: positive_cond_mass
    real(rp) :: precip_scale
 
    ! Tolerances
    real(rp) :: energy_tol
    real(rp) :: boundary_energy_tol    
    real(rp) :: water_tol
 
    integer :: ke_z, pz
    integer :: iter
    integer :: iter_water
 
    logical :: found_base
    logical :: adiabat_profile_ok
 
    logical :: water_lo_feasible
    logical :: water_hi_feasible
 
    real(rp), parameter :: mca_tbase_range = 40.0_rp
    !------------------------------------------------------------
 
    converged = .false.
    profile_status = mca_profile_energy_maxiter
 
    ! Initialize the adjusted profile with the original state
 
    temp_adj(:,:) = temp(:,:)
    qvap_adj(:,:) = qvap(:,:)
    rhoqvap_adj(:,:) = dens(:,:) * qvap(:,:)
    rhoprecip_adj(:,:) = 0.0_rp
 
    !- Calculate initial column-integrated MSE
 
    call integ_masked_column_mse( energy_target, & ! (out)
      temp, qvap, zlev, dens,           & ! (in)
      int_weight, mix_mask, npz, nez    ) ! (in)
 
    energy_tol = mca_column_energy_atol + mca_column_energy_rtol * abs(energy_target)
    boundary_energy_tol = mca_column_energy_atol + mca_boundary_energy_rtol * abs(energy_target)
 
    !- Find temperature at the lowest adjusted node
 
    tbase_mid = 0.0_rp
    found_base = .false.
 
    do ke_z = 1, nez
      do pz = 1, npz
        if ( mix_mask(pz,ke_z) ) then
          tbase_mid = temp(pz,ke_z)
          found_base = .true.
          exit
        end if
      end do
      if ( found_base ) exit
    end do
 
    if ( .not. found_base ) then
      profile_status = mca_profile_no_adjusted_node
      if ( debug_flag ) then
        write(*,*) "solve_moist_neutral_adjustment: No adjusted DG node was found."
      end if
      return
    end if
 
    !-  Calculate the total water mass in the column before adjustment
 
    water_mass_ini = 0.0_rp
 
    do ke_z = 1, nez
    do pz = 1, npz
      if ( mix_mask(pz,ke_z) ) then
        water_mass_ini = water_mass_ini + int_weight(pz,ke_z) * dens(pz,ke_z) * qvap(pz,ke_z)
      end if
    end do
    end do
 
    water_tol = mca_water_atol + mca_water_rtol * max(abs(water_mass_ini),1.0_rp)    
 
    != Step 1:
    ! Determine the upper limit of the water-feasible base temperature.
    ! A trial moist-neutral profile is water-feasible when its vapor mass does not exceed the initial vapor mass in the adjustment region.
 
    twater_lo = tbase_mid - mca_tbase_range
    twater_hi = tbase_mid + mca_tbase_range
 
    !- Evaluate the water mass at the lower bound of the temperature range
 
    pres_trial(:,:) = pres(:,:)
    call evaluate_trial_moist_neutral_profile( temp_trial, qvap_trial, water_tmin, energy_tmin, adiabat_profile_ok, & ! (out)
      pres_trial,                                                                                                 & ! (inout)
      twater_lo, temp, qvap, dens, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag  ) ! (in)
 
    if ( .not. adiabat_profile_ok ) then
      profile_status = mca_profile_adiabat_failure
      if ( debug_flag ) then
        write(*,*) "solve_moist_neutral_adjustment: Moist-adiabat integration failed. tbase_lo=", twater_lo
      end if
      return
    end if
 
    water_lo_feasible = ( water_tmin <= water_mass_ini + water_tol )
    if ( .not. water_lo_feasible ) then
      profile_status = mca_profile_no_water_feasible_state
 
      if ( debug_flag ) then
        write(*,*) "No water-feasible saturated profile exists."
        write(*,*) "  tbase_lo=", twater_lo, "water_mass_lo=", water_tmin, "water_mass_ini =", water_mass_ini
      end if
      return
    end if  
 
    ! Initialize the feasible side of the water-search bracket.
    water_wlo  = water_tmin
    energy_wlo = energy_tmin    
 
    
    !- Evaluate the water mass at the upper bound of the temperature range
 
    pres_trial(:,:) = pres(:,:)
    call evaluate_trial_moist_neutral_profile( temp_trial, qvap_trial, water_whi, energy_whi, adiabat_profile_ok, & ! (out)
      pres_trial,                                                                                     & ! (inout)
      twater_hi, temp, qvap, dens, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag  ) ! (in)
 
    if ( adiabat_profile_ok ) then
      water_hi_feasible = ( water_whi <= water_mass_ini + water_tol )
    else
      water_hi_feasible = .false.
      water_whi = huge(1.0_rp)
 
      if ( debug_flag ) then
        write(*,*) "High-temperature trial profile is not admissible."
        write(*,*) "  Treating it as the infeasible upper bracket: twater_hi=", twater_hi 
      end if
    end if
 
    !-
 
    if ( debug_flag ) then
      write(*,*) "Water feasibility at initial bounds:"
      write(*,*) "  t_lo, water_lo  =", twater_lo, water_tmin
      write(*,*) "  t_hi, water_hi  =", twater_hi, water_whi
      write(*,*) "  lower feasible  =", water_lo_feasible
      write(*,*) "  upper feasible  =", water_hi_feasible      
      write(*,*) "  water_mass_ini  =", water_mass_ini
      write(*,*) "  water_tol       =", water_tol
    end if
 
 
    if ( .not. water_hi_feasible ) then
      ! Search for the maximum base temperature for which the moist-adiabat construction succeeds 
      ! and the resulting vapor mass remains feasible.      
      ! The lower endpoint is always:
      !   - thermodynamically admissible
      !   - water feasible
      !
      ! The upper endpoint is either:
      !   - water infeasible, or
      !   - thermodynamically inadmissible
 
      pres_trial(:,:) = pres(:,:)
 
      do iter_water = 1, mca_max_water_iter
        twater_mid = 0.5_rp * ( twater_lo + twater_hi )
 
        call evaluate_trial_moist_neutral_profile( temp_trial, qvap_trial, water_mid, energy_mid, adiabat_profile_ok, & ! (out)
          pres_trial,                                                                                                   & ! (inout)
          twater_mid, temp, qvap, dens, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag         ) ! (in)
 
        if ( .not. adiabat_profile_ok ) then
          ! Thermodynamically inadmissible midpoint:
          ! move the upper boundary downward.
          twater_hi = twater_mid
          water_whi  = huge(1.0_rp)
          energy_whi = huge(1.0_rp)
 
          ! if ( debug_flag ) then
          !   write(*,*) "water iter=", iter_water
          !   write(*,*) "  inadmissible midpoint=", twater_mid, "  new bracket=", twater_lo, twater_hi
          ! end if
 
          if ( abs(twater_hi - twater_lo) <= mca_temp_tol ) exit
          cycle
        end if
 
        if ( water_mid <= water_mass_ini + water_tol ) then
          ! Admissible and water feasible.
          twater_lo = twater_mid
          water_wlo  = water_mid
          energy_wlo = energy_mid
        else
          ! Admissible but water infeasible.
          twater_hi = twater_mid
          water_whi  = water_mid
          energy_whi = energy_mid
        end if
 
        ! if ( debug_flag ) then
        !   write(*,*) "water iter=", iter_water, ", tbase=", twater_lo, twater_mid, twater_hi, ", water=", water_wlo, water_mid, water_whi, ", target=", water_mass_ini
        ! end if
 
        if ( abs(twater_hi - twater_lo) <= mca_temp_tol ) exit
        if ( abs(water_mid - water_mass_ini) <= water_tol ) exit
      end do
    end if
 
    ! Lower endpoint of the energy root search.
    tbase_lo = tbase_mid - mca_tbase_range
    water_lo = water_tmin
    energy_lo = energy_tmin    
 
    if ( water_hi_feasible ) then
      ! Initial upper bound itself is water feasible.
      tbase_hi = twater_hi
      water_hi = water_whi
      energy_hi = energy_whi
    else
      ! The final feasible lower side of the water bracket
      ! is the maximum admissible/water-feasible base temperature.
      tbase_hi = twater_lo
      water_hi = water_wlo
      energy_hi = energy_wlo
    end if
 
    !- Step 2:
    ! Evaluate the energy residual at both ends of the water-feasible interval.
 
    residual_lo = energy_lo - energy_target
    residual_hi = energy_hi - energy_target
 
    if ( debug_flag ) then
      write(*,*) "Energy-root feasibility check:"
      write(*,*) "  energy_target       =", energy_target
      write(*,*) "  energy_tol          =", energy_tol
      write(*,*) "  tbase_lo            =", tbase_lo
      write(*,*) "  energy_lo           =", energy_lo
      write(*,*) "  residual_lo         =", residual_lo
      write(*,*) "  tbase_water_max     =", tbase_hi
      write(*,*) "  energy_water_max    =", energy_hi
      write(*,*) "  residual_water_max  =", residual_hi
      write(*,*) "  water_at_upper      =", water_hi
      write(*,*) "  water_mass_ini      =", water_mass_ini
    end if    
 
    ! Check whether either endpoint is already an energy root.
    if ( abs(residual_lo) <= boundary_energy_tol ) then
      pres_trial(:,:) = pres(:,:)
      call evaluate_trial_moist_neutral_profile( temp_adj, qvap_adj, water_mass_adj, energy_trial, adiabat_profile_ok, & ! (out)
        pres_trial,                                                                                                    & ! (inout)
        tbase_lo, temp, qvap, dens, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag            ) ! (in)
 
      if ( .not. adiabat_profile_ok ) then
        profile_status = mca_profile_adiabat_failure
        return
      end if
      converged = .true.
      profile_status = mca_profile_success
 
    else if ( abs(residual_hi) <= boundary_energy_tol ) then
      pres_trial(:,:) = pres(:,:)
      call evaluate_trial_moist_neutral_profile( temp_adj, qvap_adj, water_mass_adj, energy_trial, adiabat_profile_ok, & ! (out)
        pres_trial,                                                                                                    & ! (inout)
        tbase_hi, temp, qvap, dens, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag                  ) ! (in)
 
      if ( .not. adiabat_profile_ok ) then
        profile_status = mca_profile_adiabat_failure
        return
      end if
 
      converged = .true.
      profile_status = mca_profile_success
 
    else if ( residual_lo * residual_hi > 0.0_rp ) then
      profile_status = mca_profile_no_energy_root
      if ( debug_flag ) then
        write(*,*) "No energy root exists in the water-feasible interval."
        write(*,*) "  tbase_lo           =", tbase_lo
        write(*,*) "  tbase_water_max    =", tbase_hi
        write(*,*) "  residual_lo        =", residual_lo
        write(*,*) "  residual_upper     =", residual_hi
        write(*,*) "  water_mass_ini     =", water_mass_ini
        write(*,*) "  water_mass_upper   =", water_hi
      end if
      return
    end if
 
    !- Step 3:
    ! If the endpoint check found a valid bracket, solve the column-energy constraint by bisection.
    
    if ( .not. converged ) then
      
      pres_trial(:,:) = pres(:,:)
 
      do iter=1, mca_max_energy_iter
 
        tbase_mid = 0.5_rp * ( tbase_lo + tbase_hi )
 
        call evaluate_trial_moist_neutral_profile( temp_adj, qvap_adj, water_trial, energy_trial, adiabat_profile_ok, & ! (out)
          pres_trial,                                                                                                 & ! (inout)
          tbase_mid, temp, qvap, dens, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag              ) ! (in)
 
        if ( .not. adiabat_profile_ok ) then
          profile_status = mca_profile_adiabat_failure
 
          if ( debug_flag ) then
            write(*,*) "Moist-adiabat integration failed during"
            write(*,*) "energy-root search: iteration=", iter, ", tbase=", tbase_mid
          end if
          return
        end if
 
        ! The energy search must remain inside the previously determined water-feasible interval.
 
        if ( water_trial > water_mass_ini + water_tol ) then
          profile_status = mca_profile_no_energy_root
          if ( debug_flag ) then
            write(*,*) "Energy bisection entered a water-infeasible state."
            write(*,*) "  tbase          =", tbase_mid
            write(*,*) "  water_trial    =", water_trial
            write(*,*) "  water_mass_ini =", water_mass_ini
          end if
          return
        end if
 
        residual = energy_trial - energy_target
 
        ! if ( debug_flag ) then
        !   write(*,*) "energy iter=", iter
        !   write(*,*) "  tbase=", tbase_lo, tbase_mid, tbase_hi
        !   write(*,*) "  energy_trial=", energy_trial, ", energy_target=", energy_target
        !   write(*,*) "  residual=", residual
        !   write(*,*) "  water_trial=", water_trial
        ! end if
 
        if ( abs(residual) <= energy_tol ) then
          water_mass_adj = water_trial
          converged = .true.
          profile_status = mca_profile_success
          exit
        end if
 
        ! General bisection update based on the residual signs.
 
        if ( residual_lo * residual <= 0.0_rp ) then
          tbase_hi   = tbase_mid
          residual_hi = residual
        else
          tbase_lo   = tbase_mid
          residual_lo = residual
        end if
      
      end do ! end loop for iteration
    end if
 
    if ( .not. converged ) then
      profile_status = mca_profile_energy_maxiter
 
      if ( debug_flag ) then
        write(*,*) "Energy bisection did not converge."
        write(*,*) "  tbase_lo=", tbase_lo, ", tbase_hi=", tbase_hi
        write(*,*) "  residual_lo=", residual_lo, ", residual_hi=", residual_hi
        write(*,*) "  energy_tol=", energy_tol
      end if
      return
    end if
 
    ! Step 4:
    ! Convert the converged vapor mixing ratio to density form and diagnose precipitation.
 
    rhoqvap_adj(:,:) = dens(:,:) * qvap_adj(:,:)
 
    precip_mass_total = max( 0.0_rp, water_mass_ini - water_mass_adj )
    positive_cond_mass = 0.0_rp
    do ke_z=1, nez
    do pz=1, npz
      if ( mix_mask(pz,ke_z) ) then
        precip_mass_local = max( 0.0_rp, dens(pz,ke_z) * qvap(pz,ke_z) - rhoqvap_adj(pz,ke_z) )
 
        positive_cond_mass = positive_cond_mass + precip_mass_local * int_weight(pz,ke_z)
      end if
    end do
    end do
 
    ! Distribute the column-integrated removed water mass in proportion to positive local vapor-density reductions.
 
    if ( positive_cond_mass > eps ) then
      precip_scale = precip_mass_total / positive_cond_mass
      do ke_z=1, nez
      do pz=1, npz
        if ( mix_mask(pz,ke_z) ) then
          precip_mass_local = max( 0.0_rp, dens(pz,ke_z) * qvap(pz,ke_z) - rhoqvap_adj(pz,ke_z) )
          rhoprecip_adj(pz,ke_z) = precip_scale * precip_mass_local
        end if
      end do
      end do
    end if
 
    return
  end subroutine solve_moist_neutral_adjustment
 
!OCL SERIAL
  subroutine evaluate_trial_moist_neutral_profile( &
    temp_work, qvap_work, water_mass, energy, profile_ok,  & ! (out)
    pres_work,                                             & ! (inout)
    tbase, temp, qvap, dens, zlev, int_weight,             & ! (in)
    adj_mask, forced_saturation_mask, nez, npz, debug_flag ) ! (in)
    implicit none
 
    integer, intent(in) :: nez
    integer, intent(in) :: npz
    real(rp), intent(out) :: temp_work(npz,nez)
    real(rp), intent(out) :: qvap_work(npz,nez)
    real(rp), intent(out) :: water_mass
    real(rp), intent(out) :: energy
    logical,  intent(out) :: profile_ok
    real(rp), intent(inout)  :: pres_work(npz,nez)
    real(rp), intent(in)  :: tbase
    real(rp), intent(in)  :: temp(npz,nez)
    real(rp), intent(in)  :: qvap(npz,nez)
    real(rp), intent(in)  :: zlev(npz,nez)
    real(rp), intent(in)  :: dens(npz,nez)
    real(rp), intent(in)  :: int_weight(npz,nez)
    logical,  intent(in)  :: adj_mask(npz,nez)
    logical,  intent(in)  :: forced_saturation_mask(npz,nez)
    logical,  intent(in)  :: debug_flag
 
    real(rp) :: qsat_work(npz,nez)
 
real(rp) :: pres_new(npz,nez)
real(rp) :: qdry_work
real(rp) :: rtot_work
real(rp) :: pres_err
integer, parameter :: mca_max_pres_iter = 5
real(rp), parameter :: mca_pres_rtol = 1.0e-6_rp
integer :: iter_pres
 
    integer :: ke_z, pz
    !------------------------------------------------
 
    temp_work(:,:) = temp(:,:)
    qvap_work(:,:) = qvap(:,:)
    qsat_work(:,:) = 0.0_rp  
 
do iter_pres = 1, mca_max_pres_iter
! Reset trial fields before rebuilding the profile.
temp_work(:,:) = temp(:,:)
qvap_work(:,:) = qvap(:,:)
qsat_work(:,:) = 0.0_rp
! Preserve pressure outside adj_mask.
pres_new(:,:) = pres_work(:,:)
 
    !- Evaluate the reference temperature profile over adj_mask.
 
    call build_moist_adiabat_temperature( &
      temp_work, qsat_work, profile_ok,                      & ! (inout,out)
      tbase, pres_work, zlev, adj_mask, npz, nez, debug_flag ) ! (in)
 
    if ( .not. profile_ok ) then
      water_mass = huge(1.0_rp)
      energy     = huge(1.0_rp)
      return
    end if
 
    !- Humidity construction:
    ! 1. Nearly saturated nodes in the diagnosed unstable core:
    !      qv = qsat.
    ! 2. Other nodes inside the adjustment region:
    !      retain the input qv unless the trial temperature produces supersaturation.
    ! 3. Outside the adjustment region:
    !      retain the input temperature and qv.
 
    do ke_z = 1, nez
    do pz = 1, npz
      if ( adj_mask(pz,ke_z) ) then
        if ( forced_saturation_mask(pz,ke_z) ) then
          qvap_work(pz,ke_z) = qsat_work(pz,ke_z)
        else
          qvap_work(pz,ke_z) = min(qvap(pz,ke_z), qsat_work(pz,ke_z))
        end if
      end if
      qvap_work(pz,ke_z) = max(0.0_rp, qvap_work(pz,ke_z))      
    end do
    end do    
 
 
! 4. Pressure convergence check.
!------------------------------------------------------------
  pres_err = 0.0_rp
  do ke_z = 1, nez
  do pz = 1, npz
    if ( adj_mask(pz,ke_z) ) then
      qdry_work = 1.0_rp - qvap_work(pz,ke_z)
      rtot_work = rdry * qdry_work &
                + rvap * qvap_work(pz,ke_z)
      pres_new(pz,ke_z) = dens(pz,ke_z) &
                        * rtot_work         &
                        * temp_work(pz,ke_z)
 
      pres_err = max( pres_err, &
        abs( pres_new(pz,ke_z) - pres_work(pz,ke_z) ) &
        / max( abs(pres_work(pz,ke_z)), 1.0_rp ) )
    end if
  end do
  end do
 
pres_work(:,:) = pres_new(:,:)
if ( pres_err <= mca_pres_rtol ) exit 
if ( iter_pres == mca_max_pres_iter ) then
  profile_ok = .false.
  water_mass = huge(1.0_rp)
  energy     = huge(1.0_rp)
  if ( debug_flag ) then
    write(*,*) "evaluate_trial_moist_neutral_profile: Pressure iteration failed to converge."
    write(*,*) "  pres_err=", pres_err, "MCA_PRES_RTOL=", mca_pres_rtol
  end if
  return
end if
 
end do
 
    !- Column water over the full adjustment region.
 
    water_mass = 0.0_rp
    do ke_z = 1, nez
    do pz = 1, npz
      if ( adj_mask(pz,ke_z) ) then
        water_mass = water_mass + int_weight(pz,ke_z) * dens(pz,ke_z) * qvap_work(pz,ke_z)
      end if
    end do
    end do
 
    !- Column MSE over the full adjustment region.
 
    call integ_masked_column_mse( energy, & ! (out)
      temp_work, qvap_work, zlev, dens,   & ! (in)
      int_weight, adj_mask, npz, nez      ) ! (in)
 
    return
  end subroutine evaluate_trial_moist_neutral_profile
 
  ! The column-MSE constraint is imposed while holding the density and pressure profiles fixed. 
  ! Thermodynamic composition includes dry air and water vapor only. 
  ! Condensate sensible heat and the energy carried away by precipitation are neglected.  
!OCL SERIAL
  subroutine integ_masked_column_mse( energy, & ! (out)
    temp, qv, zlev, dens,                     & ! (in)
    int_weight, integ_mask, npz, nez          ) ! (in)
    implicit none
    integer, intent(in) :: npz
    integer, intent(in) :: nez
    real(rp), intent(out) :: energy
    real(rp), intent(in) :: temp(npz,nez)
    real(rp), intent(in) :: qv(npz,nez)
    real(rp), intent(in) :: zlev(npz,nez)
    real(rp), intent(in) :: dens(npz,nez)
    real(rp), intent(in) :: int_weight(npz,nez)
    logical, intent(in) :: integ_mask(npz,nez)
 
    integer :: ke_z, pz
    real(rp) :: qdry
    real(rp) :: cptot
    !------------------------------------------------
 
    energy = 0.0_rp
 
    do ke_z = 1, nez
    do pz = 1, npz
      if ( integ_mask(pz,ke_z) ) then
        qdry = 1.0_rp - qv(pz,ke_z) ! - qcon_(pz_,ke_z_)
 
        cptot = cpdry    * qdry         &
              + cp_vapor * qv(pz,ke_z)
                                        ! Liquid-water sensible heat is neglected.
        
        energy = energy + int_weight(pz,ke_z) * dens(pz,ke_z) * &
          ( cptot * temp(pz,ke_z) + grav * zlev(pz,ke_z)        &
          + lhv * qv(pz,ke_z)                                   )
      end if
    end do
    end do
    return
  end subroutine integ_masked_column_mse
 
  !----
!OCL SERIAL
  subroutine build_moist_adiabat_temperature( &
    temp_profile, qsat_profile, profile_converged, & ! (out)
    tbase, pres, zlev, mask, npz, nez, debug_flag  ) ! (in)
    implicit none
    integer, intent(in) :: npz
    integer, intent(in) :: nez
    real(rp), intent(inout) :: temp_profile(npz,nez)
    real(rp), intent(inout) :: qsat_profile(npz,nez)
    logical, intent(out) :: profile_converged
    real(rp), intent(in) :: tbase
    real(rp), intent(in) :: pres(npz,nez)
    real(rp), intent(in) :: zlev(npz,nez)
    logical, intent(in) :: mask(npz,nez)
    logical, intent(in) :: debug_flag
 
    integer :: ke_z, pz
 
    real(rp) :: t_prev, p_prev, z_prev
    real(rp) :: t_now, p_now, z_now
 
    logical :: started
    logical :: is_converge_local
    !----------------------------------------------------
 
    started = .false.
    profile_converged = .true.
 
    do ke_z = 1, nez
    do pz = 1, npz
      if ( .not. mask(pz,ke_z) ) cycle
 
      z_now = zlev(pz,ke_z)
      p_now = pres(pz,ke_z)
 
      if ( .not. started ) then
        t_now = tbase
        started = .true.
      ! else if ( abs(z_now-z_prev) <= MCA_Z_TOL ) then
      !   t_now = t_prev
      else
        call solve_next_moist_adiabat_temp( t_now, is_converge_local, & ! (out)
          t_prev, p_prev, z_prev, p_now, z_now                        ) ! (in)
        
        if ( .not. is_converge_local ) then
          profile_converged = .false.
          return
        end if
      end if
 
      temp_profile(pz,ke_z) = t_now
      call atmos_saturation_pres2qsat_liq( t_now, p_now, & ! (in)
        qsat_profile(pz,ke_z)                            ) ! (out)
 
      t_prev = t_now 
      p_prev = p_now
      z_prev = z_now
    end do
    end do
 
    return
  end subroutine build_moist_adiabat_temperature
 
  !> Solve the temperature at (pres1,zlev1) such that the saturated MSE equals that at (temp0,pres0,zlev0).  
!OCL SERIAL
  subroutine solve_next_moist_adiabat_temp( temp1, converged, & ! (out)
    temp0, pres0, zlev0, pres1, zlev1                         ) ! (in)
    implicit none
    real(rp), intent(out) :: temp1
    logical,  intent(out) :: converged
    real(rp), intent(in)  :: temp0
    real(rp), intent(in)  :: pres0
    real(rp), intent(in)  :: zlev0
    real(rp), intent(in)  :: pres1
    real(rp), intent(in)  :: zlev1
 
    real(rp) :: temp_lo, temp_mid, temp_hi
    real(rp) :: temp_cur, temp_new
    real(rp) :: f_lo, f_hi, f_cur
    real(rp) :: hs_target, hs_cur
    real(rp) :: dhs_dtemp
    real(rp) :: temp_tol,  energy_tol
 
    integer :: iter
 
    integer, parameter :: max_iter = 60
    real(rp), parameter :: temp_range = 40.0_rp
 
    logical :: use_newton
    !----------------------------------------------------
 
    converged = .false.
 
    call saturated_mse_point( temp0, pres0, zlev0, &
      hs_target )
 
    temp_lo = temp0 - temp_range
    temp_hi = temp0 + temp_range
 
    call saturated_mse_point( temp_lo, pres1, zlev1, &
      hs_cur )
    f_lo = hs_cur - hs_target
 
    call saturated_mse_point( temp_hi, pres1, zlev1, &
      hs_cur )
    f_hi = hs_cur - hs_target
 
    if ( f_lo * f_hi > 0.0_rp ) then
      return
    end if
 
    energy_tol = mca_local_hmse_atol + mca_local_hmse_rtol * abs(hs_target)
    temp_tol   = mca_temp_tol
 
    temp_cur = min( max(temp0, temp_lo), temp_hi )
 
    do iter = 1, max_iter
      
      call saturated_mse_point_with_derivative( temp_cur, pres1, zlev1, &
        hs_cur, dhs_dtemp )
      
      f_cur = hs_cur - hs_target
      if ( abs(f_cur) <= energy_tol ) then
        temp1 = temp_cur
        converged = .true.
        return
      end if
 
      ! Update bracket using current function value.      
      if ( f_lo * f_cur <= 0.0_rp ) then
        temp_hi = temp_cur
        f_hi    = f_cur
      else
        temp_lo = temp_cur
        f_lo    = f_cur
      end if
 
      ! A sufficiently small bracket is also convergence.
      if ( abs(temp_hi - temp_lo) <= temp_tol ) then
        temp1 = 0.5_rp * (temp_lo + temp_hi)
        converged = .true.
        return
      end if      
 
      ! Safeguarded Newton step.
 
      use_newton = .false.
      if ( abs(dhs_dtemp) > eps ) then
        temp_new = temp_cur - f_cur / dhs_dtemp
        if ( temp_new > temp_lo .and. temp_new < temp_hi ) then
          use_newton = .true.
        end if
      end if
 
      if ( .not. use_newton ) then
        temp_new = 0.5_rp * (temp_lo + temp_hi)
      end if
 
      temp_cur = temp_new
    end do
 
    temp1 = temp_cur
    return
  contains
    subroutine saturated_mse_point( &
      temp, pres, zlev, hmse_sat )
      implicit none
      real(rp), intent(in)  :: temp
      real(rp), intent(in)  :: pres
      real(rp), intent(in)  :: zlev
      real(rp), intent(out) :: hmse_sat
 
      real(rp) :: qsat
      real(rp) :: cptot
      !------------------------------------------------
      call atmos_saturation_pres2qsat_liq( &
        temp, pres, qsat )
 
      cptot = cpdry * ( 1.0_rp - qsat ) + cp_vapor * qsat
 
      hmse_sat = cptot * temp + grav * zlev  + lhv * qsat
      return
    end subroutine saturated_mse_point
  end subroutine solve_next_moist_adiabat_temp
 
  subroutine saturated_mse_point_with_derivative( &
      temp, pres, zlev, hmse_sat, dhmse_dtemp )
      use scale_atmos_saturation, only: atmos_saturation_dqs_dtem_dpre_liq
      implicit none
      real(rp), intent(in)  :: temp
      real(rp), intent(in)  :: pres
      real(rp), intent(in)  :: zlev
      real(rp), intent(out) :: hmse_sat
      real(rp), intent(out) :: dhmse_dtemp
 
      real(rp) :: qsat
      real(rp) :: dqsat_dtemp
      real(rp) :: dqsat_dpres
      real(rp) :: cptot
      real(rp) :: qdry_dummy
      !--------------------------------------------------
 
      call atmos_saturation_pres2qsat_liq( temp, pres, &
        qsat )
      call atmos_saturation_dqs_dtem_dpre_liq( temp, pres, qdry_dummy, &
        dqsat_dtemp, dqsat_dpres )
      
      cptot = cpdry * (1.0_rp - qsat) &
            + cp_vapor * qsat
 
      hmse_sat = cptot * temp + grav * zlev  + lhv * qsat
 
      ! h_s = cp(qs) T + gz + Lv qs
      !   where cp(qs) = Cp_d + (Cp_v-Cp_d) qs
      !
      ! * dh_s/dT = cp(qs) + T (Cp_v-Cp_d) dqs/dT + Lv dqs/dT
 
      dhmse_dtemp = cptot &
                  + ( ( cp_vapor - cpdry ) * temp + lhv ) * dqsat_dtemp
      return
    end subroutine saturated_mse_point_with_derivative  
end module scale_atm_phy_cp_dgm_mconv_adjustment