scale_atm_phy_cp_dgm_mconv_adjustment Module Reference

module FElib / Atmosphere / Physics cumulus parameterization More...

Functions/Subroutines
subroutine, public	atm_phy_cp_dgm_mconv_adjustment_setup ()
	Setup a module for moist convective adjustment scheme.
subroutine, public	atm_phy_cp_dgm_mconv_adjustment_calc_tendency (dens_t, rhot_t, rhoqv_t, sflx_rain, sflx_engi, ddens, drhot, qv, pt, pres, dens_hyd, rtot, cptot, dtsec, lmesh, elem, elem1d)
	Calculate tendencies with moist convective adjustment scheme.
subroutine, public	atm_phy_cp_dgm_mconv_adjustment_finalize ()
	Finalize a module for moist convective adjustment scheme.

Detailed Description

module FElib / Atmosphere / Physics cumulus parameterization

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

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.

Function/Subroutine Documentation

atm_phy_cp_dgm_mconv_adjustment_setup()

subroutine, public scale_atm_phy_cp_dgm_mconv_adjustment::atm_phy_cp_dgm_mconv_adjustment_setup

Setup a module for moist convective adjustment scheme.

Definition at line 109 of file scale_atm_phy_cp_dgm_mconv_adjustment.F90.

    implicit none
    !----------------------------------------------------
    return

atm_phy_cp_dgm_mconv_adjustment_calc_tendency()

subroutine, public scale_atm_phy_cp_dgm_mconv_adjustment::atm_phy_cp_dgm_mconv_adjustment_calc_tendency	(	real(rp), dimension(elem%np,lmesh%nea), intent(out)	dens_t,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	rhot_t,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	rhoqv_t,
		real(rp), dimension(elem%nnode_h1d**2,lmesh%ne2da), intent(out)	sflx_rain,
		real(rp), dimension(elem%nnode_h1d**2,lmesh%ne2da), intent(out)	sflx_engi,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	ddens,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	drhot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	qv,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pt,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pres,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dens_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	rtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	cptot,
		real(rp), intent(in)	dtsec,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		class(elementbase1d), intent(in)	elem1d )

Calculate tendencies with moist convective adjustment scheme.

Definition at line 117 of file scale_atm_phy_cp_dgm_mconv_adjustment.F90.

 
    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

atm_phy_cp_dgm_mconv_adjustment_finalize()

subroutine, public scale_atm_phy_cp_dgm_mconv_adjustment::atm_phy_cp_dgm_mconv_adjustment_finalize

Finalize a module for moist convective adjustment scheme.

Definition at line 539 of file scale_atm_phy_cp_dgm_mconv_adjustment.F90.

    implicit none
    !----------------------------------------------------
    return