mod_atmos_phy_mp.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module mod_atmos_phy_mp
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !  
  use scale_precision
  use scale_prc
  use scale_io
  use scale_prof
  use scale_const, only: &
    undef8 => const_undef8
  use scale_tracer, only: qa
 
  use scale_sparsemat, only: sparsemat
  use scale_element_hexahedral, only: hexahedralelement
 
  use scale_mesh_base, only: meshbase
  use scale_mesh_base2d, only: meshbase2d
  use scale_mesh_base3d, only: meshbase3d
 
  use scale_localmesh_base, only: localmeshbase
  use scale_localmesh_2d, only: localmesh2d
  use scale_localmesh_3d, only: localmesh3d
  use scale_element_base, only: elementbase, &
    elementbase1d, elementbase2d, elementbase3d
    
  use scale_meshfield_base, only: &
    meshfieldbase, meshfield3d
  use scale_localmeshfield_base, only: &
    localmeshfieldbase, localmeshfieldbaselist
 
  use scale_model_mesh_manager, only: modelmeshbase
  use scale_model_var_manager, only: modelvarmanager
  use scale_model_component_proc, only:  modelcomponentproc
  
  use mod_atmos_phy_mp_vars, only: atmosphympvars
 
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
  type, extends(modelcomponentproc), public :: atmosphymp
    integer :: mp_typeid
    type(atmosphympvars) :: vars
 
    logical, private :: do_precipitation
    logical, private :: do_negative_fixer
    real(rp), private :: limit_negative
    integer, private :: ntmax_sedimentation
    real(rp), private :: max_term_vel
    real(rp), private :: cldfrac_thleshold
 
    real(dp), private :: dtsec
    integer, private :: nstep_sedmientation
    real(rp), private :: rnstep_sedmientation
    real(dp), private :: dtsec_sedmientation
 
    type(sparsemat) :: dz, lift
    type(hexahedralelement) :: elem
  contains
    procedure, public :: setup => atmosphymp_setup 
    procedure, public :: calc_tendency => atmosphymp_calc_tendency
    procedure, public :: update => atmosphymp_update
    procedure, public :: finalize => atmosphymp_finalize
    procedure, private :: calc_tendency_core => atmosphymp_calc_tendency_core
  end type atmosphymp
 
  !-----------------------------------------------------------------------------
  !++ Public parameters & variables
  !-----------------------------------------------------------------------------
 
  integer, parameter :: mp_typeid_kessler = 1
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
    
contains
 
  subroutine atmosphymp_setup( this, model_mesh, tm_parent_comp )
    use scale_const, only: &
      eps => const_eps
    use scale_atmos_hydrometeor, only: &
      atmos_hydrometeor_regist
    use scale_atmos_phy_mp_kessler, only: &
      atmos_phy_mp_kessler_setup,               &
      atmos_phy_mp_kessler_ntracers,            &
      atmos_phy_mp_kessler_nwaters,             &
      atmos_phy_mp_kessler_nices,               &
      atmos_phy_mp_kessler_tracer_names,        &
      atmos_phy_mp_kessler_tracer_descriptions, &
      atmos_phy_mp_kessler_tracer_units
    !------------------------------
    use scale_time_manager, only: time_manager_component    
    use mod_atmos_mesh, only: atmosmesh
 
    implicit none
    class(atmosphymp), intent(inout) :: this
    class(modelmeshbase), target, intent(in) :: model_mesh
    class(time_manager_component), intent(inout) :: tm_parent_comp
 
    real(dp) :: time_dt                             = undef8 
    character(len=H_SHORT) :: time_dt_unit          = 'SEC'
 
    character(len=H_MID) :: mp_type = 'KESSLER'
 
    logical :: do_precipitation
    logical :: do_negative_fixer
    real(rp) :: limit_negative
    integer :: ntmax_sedimentation  
    real(rp) :: max_term_vel        
    real(rp) :: cldfrac_thleshold
 
    namelist /param_atmos_phy_mp/ &
      time_dt,             &
      time_dt_unit,        &
      mp_type,             &
      do_precipitation,    &
      do_negative_fixer,   &
      limit_negative,      &
      ntmax_sedimentation, &
      max_term_vel,        &
      cldfrac_thleshold
    
    class(atmosmesh), pointer     :: atm_mesh
    class(meshbase), pointer      :: ptr_mesh
 
    integer :: ierr
 
    integer :: qs_mp, qe_mp, qa_mp
 
    integer :: nstep_max
    !--------------------------------------------------
 
    if (.not. this%IsActivated()) return
 
    log_newline
    log_info("ATMOS_PHY_MP_setup",*) 'Setup'
 
    cldfrac_thleshold = eps
 
    do_precipitation    = .true.
    do_negative_fixer   = .true.
    limit_negative      = 0.1_rp
    ntmax_sedimentation = 1
    max_term_vel        = 10.0_rp
    
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_mp,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("ATMOS_PHY_MP_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
      log_error("ATMOS_PHY_MP_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_MP. Check!'
      call prc_abort
    endif
    log_nml(param_atmos_phy_mp)
 
    this%cldfrac_thleshold = cldfrac_thleshold
    this%do_precipitation  = do_precipitation
    this%do_negative_fixer = do_negative_fixer
    this%ntmax_sedimentation = ntmax_sedimentation
    this%max_term_vel        = max_term_vel
 
    log_newline
    log_info("ATMOS_PHY_MP_setup",*) 'Enable negative fixer?                    : ', this%do_negative_fixer
    log_info("ATMOS_PHY_MP_setup",*) 'Value limit of negative fixer (abs)       : ', abs(this%limit_negative)
    log_info("ATMOS_PHY_MP_setup",*) 'Enable sedimentation (precipitation)?     : ', this%do_precipitation
 
    !- get mesh --------------------------------------------------
 
    call model_mesh%GetModelMesh( ptr_mesh )
    select type(model_mesh)
    class is (atmosmesh)
      atm_mesh => model_mesh
    end select
 
    !--- Regist this compoent in the time manager
    
    call tm_parent_comp%Regist_process( 'ATMOS_PHY_MP', time_dt, time_dt_unit, & ! (in)
      this%tm_process_id )                                                       ! (out) 
 
    this%dtsec = tm_parent_comp%process_list(this%tm_process_id)%dtsec
    nstep_max = 0 ! ceiling( this%dtsec * this%max_term_vel / maxval( CDZ) )
    this%ntmax_sedimentation = max( ntmax_sedimentation, nstep_max )
 
    this%nstep_sedmientation = ntmax_sedimentation
    this%rnstep_sedmientation = 1.0_rp / real(ntmax_sedimentation,kind=rp)
    this%dtsec_sedmientation  = this%dtsec * this%rnstep_sedmientation
  
    log_newline
    log_info("ATMOS_PHY_MP_setup",*) 'Enable negative fixer?                    : ', this%do_negative_fixer
    !LOG_INFO("ATMOS_PHY_MP_setup",*) 'Value limit of negative fixer (abs)       : ', abs(MP_limit_negative)
    log_info("ATMOS_PHY_MP_setup",*) 'Enable sedimentation (precipitation)?     : ', this%do_precipitation
    log_info("ATMOS_PHY_MP_setup",*) 'Timestep of sedimentation is divided into : ', this%ntmax_sedimentation, 'step'
    log_info("ATMOS_PHY_MP_setup",*) 'DT of sedimentation                       : ', this%dtsec_sedmientation, '[s]'
 
    !--- Set the type of microphysics scheme
 
    select case( mp_type )
    case ('KESSLER')
      this%MP_TYPEID = mp_typeid_kessler
      call atmos_hydrometeor_regist( &
           atmos_phy_mp_kessler_nwaters,                & ! (in)
           atmos_phy_mp_kessler_nices,                  & ! (in)
           atmos_phy_mp_kessler_tracer_names(:),        & ! (in)
           atmos_phy_mp_kessler_tracer_descriptions(:), & ! (in)
           atmos_phy_mp_kessler_tracer_units(:),        & ! (in)
           qs_mp                                        ) ! (out)
      qa_mp = atmos_phy_mp_kessler_ntracers
  
    case default
      log_error("ATMOS_PHY_MP_setup",*) 'Not appropriate names of MP_TYPE in namelist PARAM_ATMOS_PHY_MP. Check!'
      call prc_abort
    end select
 
    qe_mp = qs_mp + qa_mp - 1
 
    !- initialize the variables 
    call this%vars%Init( model_mesh, qs_mp, qe_mp, qa_mp )     
 
    !- Setup a module for microcloud physics
 
    select case( this%MP_TYPEID )
    case( mp_typeid_kessler )
      call atmos_phy_mp_kessler_setup()
    end select
 
    !
    call this%elem%Init( atm_mesh%ptr_mesh%refElem3D%PolyOrder_h, atm_mesh%ptr_mesh%refElem3D%PolyOrder_v, .true. )
    call this%Dz%Init( this%elem%Dx3, storage_format='ELL' )
    call this%Lift%Init( this%elem%Lift, storage_format='ELL' )
 
    return
  end subroutine atmosphymp_setup
 
!OCL SERIAL
  subroutine atmosphymp_calc_tendency( &
    this, model_mesh, prgvars_list, trcvars_list, &
    auxvars_list, forcing_list, is_update         )
 
    use mod_atmos_vars, only: &
      atmosvars_getlocalmeshprgvars,     &
      atmosvars_getlocalmeshphyauxvars,  &
      atmosvars_getlocalmeshqtrcvarlist, & 
      atmosvars_getlocalmeshphytends
    use mod_atmos_phy_mp_vars, only:           &
      atmosphympvars_getlocalmeshfields_tend,  &
      atmosphympvars_getlocalmeshfields_sfcflx
    
    implicit none
    
    class(atmosphymp), intent(inout) :: this
    class(modelmeshbase), intent(in) :: model_mesh
    class(modelvarmanager), intent(inout) :: prgvars_list
    class(modelvarmanager), intent(inout) :: trcvars_list    
    class(modelvarmanager), intent(inout) :: auxvars_list
    class(modelvarmanager), intent(inout) :: forcing_list
    logical, intent(in) :: is_update
 
    class(meshbase), pointer :: mesh
    class(localmesh3d), pointer :: lcmesh
 
    class(localmeshfieldbase), pointer :: ddens, momx, momy, momz, drhot
    class(localmeshfieldbase), pointer :: rtot, cvtot, cptot
    type(localmeshfieldbaselist) :: qtrc(this%vars%qs:this%vars%qe)
    class(localmeshfieldbase), pointer :: dens_hyd, pres_hyd
    class(localmeshfieldbase), pointer :: pres, pt
    class(localmeshfieldbase), pointer :: dens_tp, momx_tp, momy_tp, momz_tp, rhot_tp, rhoh_p
    type(localmeshfieldbaselist) :: rhoq_tp(qa)
    class(localmeshfieldbase), pointer :: mp_dens_t, mp_momx_t, mp_momy_t, mp_momz_t, mp_rhot_t, mp_rhoh, mp_evap
    type(localmeshfieldbaselist) :: mp_rhoq_t(this%vars%qs:this%vars%qe)
    class(localmeshfieldbase), pointer :: sflx_rain, sflx_snow, sflx_engi
    
    integer :: n
    integer :: ke
    integer :: iq
    !--------------------------------------------------
 
    if (.not. this%IsActivated()) return
 
    log_progress(*) 'atmosphere / physics / cloud microphysics'
 
    call model_mesh%GetModelMesh( mesh )
 
    do n=1, mesh%LOCAL_MESH_NUM
      call prof_rapstart('ATM_PHY_MP_get_localmesh_ptr', 2)         
      call atmosvars_getlocalmeshprgvars( n,  &
        mesh, prgvars_list, auxvars_list,       &
        ddens, momx, momy, momz, drhot,         &
        dens_hyd, pres_hyd, rtot, cvtot, cptot, &
        lcmesh                                  )
      
      call atmosvars_getlocalmeshqtrcvarlist( n,  &
        mesh, trcvars_list, this%vars%QS, qtrc(:) )
      
      call atmosvars_getlocalmeshphyauxvars( n,  &
        mesh, auxvars_list,                      &
        pres, pt )
      
      call atmosvars_getlocalmeshphytends( n,        &
        mesh, forcing_list,                          &
        dens_tp, momx_tp, momy_tp, momz_tp, rhot_tp, &
        rhoh_p, rhoq_tp  )
 
      call atmosphympvars_getlocalmeshfields_tend( n, &
        mesh, this%vars%tends_manager,                                           &
        mp_dens_t, mp_momx_t, mp_momy_t, mp_momz_t, mp_rhot_t, mp_rhoh, mp_evap, &
        mp_rhoq_t                                                                )
 
      call atmosphympvars_getlocalmeshfields_sfcflx( n, &
        mesh, this%vars%auxvars2D_manager,              &
        sflx_rain, sflx_snow, sflx_engi                 )     
 
      call prof_rapend('ATM_PHY_MP_get_localmesh_ptr', 2)   
 
      call prof_rapstart('ATM_PHY_MP_cal_tend', 2)
      if (is_update) then
        call this%calc_tendency_core( &
          mp_dens_t%val, mp_momx_t%val, mp_momy_t%val, mp_momz_t%val, mp_rhoq_t,  & ! (out)
          mp_rhoh%val, mp_evap%val, sflx_rain%val, sflx_snow%val, sflx_engi%val,  & ! (out)
          ddens%val, momx%val, momy%val, momz%val, pt%val, qtrc,                  & ! (in)
          pres%val, pres_hyd%val, dens_hyd%val, rtot%val, cvtot%val, cptot%val,   & ! (in)
          model_mesh%DOptrMat(3), model_mesh%LiftOptrMat,                         & ! (in)
          lcmesh, lcmesh%refElem3D, lcmesh%lcmesh2D, lcmesh%lcmesh2D%refElem2D    ) ! (in)
      end if
      
      !$omp parallel private(ke, iq)
      !$omp do
      do ke = lcmesh%NeS, lcmesh%NeE
        dens_tp%val(:,ke) = dens_tp%val(:,ke) + mp_dens_t%val(:,ke)
        momx_tp%val(:,ke) = momx_tp%val(:,ke) + mp_momx_t%val(:,ke)
        momy_tp%val(:,ke) = momy_tp%val(:,ke) + mp_momy_t%val(:,ke)
        momz_tp%val(:,ke) = momz_tp%val(:,ke) + mp_momz_t%val(:,ke)
        rhoh_p %val(:,ke) = rhoh_p %val(:,ke) + mp_rhoh  %val(:,ke)
      end do
      !$omp end do
      !$omp do collapse(2)
      do iq = this%vars%QS, this%vars%QE
      do ke = lcmesh%NeS, lcmesh%NeE
        rhoq_tp(iq)%ptr%val(:,ke) = rhoq_tp(iq)%ptr%val(:,ke)  &
                                  + mp_rhoq_t(iq)%ptr%val(:,ke)
      end do
      end do 
      !$omp end do
      !$omp end parallel
 
      call prof_rapend('ATM_PHY_MP_cal_tend', 2)
    end do
 
    return  
  end subroutine atmosphymp_calc_tendency
 
!OCL SERIAL  
  subroutine atmosphymp_update( this, model_mesh, &
    prgvars_list, trcvars_list,                   &
    auxvars_list, forcing_list, is_update         )  
    
    implicit none
    class(atmosphymp), intent(inout) :: this
    class(modelmeshbase), intent(in) :: model_mesh
    class(modelvarmanager), intent(inout) :: prgvars_list
    class(modelvarmanager), intent(inout) :: trcvars_list    
    class(modelvarmanager), intent(inout) :: auxvars_list
    class(modelvarmanager), intent(inout) :: forcing_list
    logical, intent(in) :: is_update
    !--------------------------------------------------
 
    return
  end subroutine atmosphymp_update
 
!OCL SERIAL  
  subroutine atmosphymp_finalize( this )
    implicit none
    class(atmosphymp), intent(inout) :: this
 
    !--------------------------------------------------
    if (.not. this%IsActivated()) return
 
    select case ( this%MP_TYPEID )
    case( mp_typeid_kessler )
    end select
 
    call this%vars%Final()
 
    return
  end subroutine atmosphymp_finalize
  
!- private ------------------------------------------------
 
!OCL SERIAL
  subroutine atmosphymp_calc_tendency_core( this, &
    DENS_t_MP, RHOU_t_MP, RHOV_t_MP, MOMZ_t_MP, RHOQ_t_MP, & ! (out)
    rhoh_mp, evaporate, sflx_rain, sflx_snow, sflx_engi,   & ! (out)
    ddens, rhou, rhov, momz, pt, qtrc,                     & ! (in)
    pres, pres_hyd, dens_hyd, rtot, cvtot, cptot,          & ! (in)
    dz, lift,                                              & ! (in)
    lcmesh, elem3d, lcmesh2d, elem2d                       ) ! (in)
 
    use scale_const, only: &
      pre00 => const_pre00
    use scale_const, only: &
      cvdry => const_cvdry, &
      cpdry => const_cpdry
    use scale_atmos_hydrometeor, only: & 
      lhf,                             &
      qha, qla, qia
    use scale_atmos_phy_mp_kessler, only:    &
      atmos_phy_mp_kessler_adjustment,       &
      atmos_phy_mp_kessler_terminal_velocity
 
    use scale_sparsemat, only: sparsemat
    use scale_atm_phy_mp_dgm_common, only:  &
      atm_phy_mp_dgm_common_gen_vmap,              &
      atm_phy_mp_dgm_common_gen_intweight,         &
      atm_phy_mp_dgm_common_precipitation,         &
      atm_phy_mp_dgm_common_precipitation_momentum
 
    implicit none
 
    class(atmosphymp), intent(inout) :: this
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem3d
    class(localmesh2d), intent(in) :: lcmesh2d
    class(elementbase2d), intent(in) :: elem2d
    real(rp), intent(out) :: dens_t_mp(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: rhou_t_mp(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: rhov_t_mp(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: momz_t_mp(elem3d%np,lcmesh%nea)
    type(localmeshfieldbaselist), intent(inout) :: rhoq_t_mp(this%vars%qs:this%vars%qe)
    real(rp), intent(out) :: rhoh_mp(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: evaporate(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: sflx_rain(elem2d%np,lcmesh2d%nea) 
    real(rp), intent(out) :: sflx_snow(elem2d%np,lcmesh2d%nea) 
    real(rp), intent(out) :: sflx_engi(elem2d%np,lcmesh2d%nea) 
    real(rp), intent(in) :: ddens(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: rhou(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: rhov(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: momz(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: pt  (elem3d%np,lcmesh%nea)
    type(localmeshfieldbaselist), intent(in) :: qtrc(this%vars%qs:this%vars%qe)
    real(rp), intent(in) :: pres(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: pres_hyd(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: dens_hyd(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: rtot (elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: cvtot(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: cptot(elem3d%np,lcmesh%nea)
    class(sparsemat), intent(in) :: dz
    class(sparsemat), intent(in) :: lift
 
    real(rp) :: rhoe_t(elem3d%np,lcmesh%nea)
 
    real(rp) :: dens (elem3d%np,lcmesh%nea)
    real(rp) :: temp1(elem3d%np,lcmesh%nea)
    real(rp) :: cptot1(elem3d%np,lcmesh%nea)
    real(rp) :: cvtot1(elem3d%np,lcmesh%nea)
    real(rp) :: qtrc1(elem3d%np,lcmesh%nea,this%vars%qs:this%vars%qe)
 
    real(rp) :: cptot_t(elem3d%np,lcmesh%nea)
    real(rp) :: cvtot_t(elem3d%np,lcmesh%nea)
 
    real(rp) :: vterm(elem3d%np,lcmesh%nez,lcmesh%ne2d,this%vars%qs+1:this%vars%qe)
    real(rp) :: dens0(elem3d%np,lcmesh%nez,lcmesh%ne2d)    
    real(rp) :: dens2(elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: ref_dens(elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: rhou2(elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: rhov2(elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: momz2(elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: temp2(elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: pres2(elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: cptot2(elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: cvtot2(elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: rhoe (elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: rhoe2(elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: rhoq (elem3d%np,lcmesh%nez,lcmesh%ne2d,this%vars%qs+1:this%vars%qe)
    real(rp) :: rhoq2(elem3d%np,lcmesh%nez,lcmesh%ne2d,this%vars%qs+1:this%vars%qe)
    real(rp) :: rhoq2_tmp(elem3d%nnode_v,lcmesh%nez,this%vars%qs+1:this%vars%qe)
    real(rp) :: dens2_tmp(elem3d%nnode_v,lcmesh%nez)
    real(rp) :: ref_dens_tmp(elem3d%nnode_v,lcmesh%nez)
    real(rp) :: vterm_tmp(elem3d%nnode_v,lcmesh%nez,this%vars%qs+1:this%vars%qe)
    real(rp) :: flx_hydro(elem3d%np,lcmesh%nez,lcmesh%ne2d)
    real(rp) :: cp_t(elem3d%np), cv_t(elem3d%np)
 
    integer :: iq
    integer :: domid
    integer :: ke
    integer :: ke2d, ke_z
    integer :: p2d, pv1d, p
    integer :: colmask(elem3d%nnode_v)
 
    integer :: step
    real(rp) :: rdt_mp
 
    integer :: vmapm(elem3d%nfptot,lcmesh%nez)
    integer :: vmapp(elem3d%nfptot,lcmesh%nez)
    real(rp) :: intweight(elem3d%nfaces,elem3d%nfptot)
    real(rp) :: nz(elem3d%nfptot,lcmesh%nez,lcmesh%ne2d)
    !----------------------------------------------
 
    rdt_mp = 1.0_rp / this%dtsec
    domid  = lcmesh%lcdomID
 
    call atm_phy_mp_dgm_common_gen_vmap( vmapm, vmapp, & ! (out)
      lcmesh, elem3d                                   ) ! (in)
 
    call atm_phy_mp_dgm_common_gen_intweight( intweight, & ! (out)
      lcmesh                                             ) ! (in)
    
    select case( this%MP_TYPEID )
    case( mp_typeid_kessler )
 
      !$omp parallel private(ke, iq)
      !$omp do
      do ke = lcmesh%NeS, lcmesh%NeE
        dens(:,ke) = dens_hyd(:,ke) + ddens(:,ke)
        temp1(:,ke) = pres(:,ke) / ( dens(:,ke) * rtot(:,ke) )
        cptot1(:,ke) = cptot(:,ke)
        cvtot1(:,ke) = cvtot(:,ke)
      end do
      !$omp end do
      !$omp do collapse(2)
      do iq = this%vars%QS, this%vars%QE
      do ke = lcmesh%NeS, lcmesh%NeE
        qtrc1(:,ke,iq) = qtrc(iq)%ptr%val(:,ke)
      end do
      end do
      !$omp end do
      !$omp end parallel
 
      call atmos_phy_mp_kessler_adjustment( &
        elem3d%Np, 1, elem3d%Np, lcmesh%NeA, lcmesh%NeS, lcmesh%NeE, 1, 1, 1,  & ! (in)
        dens, pres, this%dtsec,                                                & ! (in)
        temp1, qtrc1, cptot1, cvtot1,                                          & ! (inout)
        rhoe_t, evaporate                                                      ) ! (out)
    
     !$omp parallel private(ke, iq)
     !$omp do collapse(2)
      do iq = this%vars%QS, this%vars%QE
      do ke = lcmesh%NeS, lcmesh%NeE
        rhoq_t_mp(iq)%ptr%val(:,ke) = ( qtrc1(:,ke,iq) - qtrc(iq)%ptr%val(:,ke) ) * dens(:,ke) * rdt_mp
      end do  
      end do
     !$omp end do
     !$omp do
      do ke = lcmesh%NeS, lcmesh%NeE
        cptot_t(:,ke) = ( cptot1(:,ke) - cptot(:,ke) ) * rdt_mp
        cvtot_t(:,ke) = ( cvtot1(:,ke) - cvtot(:,ke) ) * rdt_mp
      end do
     !$omp end do
     !$omp end parallel
 
    end select
 
    !$omp parallel do
    do ke = lcmesh%NeS, lcmesh%NeE
      rhoh_mp(:,ke) = rhoe_t(:,ke) &
        - ( cptot_t(:,ke) + log( pres(:,ke) / pre00 ) * ( cvtot(:,ke) / cptot(:,ke) * cptot_t(:,ke) - cvtot_t(:,ke) ) ) &
        * pres(:,ke) / rtot(:,ke)
    end do
 
    if ( this%do_precipitation ) then
 
      !$omp parallel private(ke,ke2D,ke_z,iq)
      !$omp do collapse(2)
      do ke2d = 1, lcmesh%Ne2D
      do ke_z = 1, lcmesh%NeZ
        ke = ke2d + (ke_z-1)*lcmesh%Ne2D
 
        dens0(:,ke_z,ke2d) = dens(:,ke)        
        dens2(:,ke_z,ke2d) = dens(:,ke)
        ref_dens(:,ke_z,ke2d) = dens_hyd(:,ke)
        rhou2(:,ke_z,ke2d) = rhou(:,ke)
        rhov2(:,ke_z,ke2d) = rhov(:,ke)
        momz2(:,ke_z,ke2d) = momz(:,ke)
        temp2(:,ke_z,ke2d) = pres(:,ke) / ( dens(:,ke) * rtot(:,ke) )
        pres2(:,ke_z,ke2d) = pres(:,ke)
        cptot2(:,ke_z,ke2d) = cptot(:,ke)
        cvtot2(:,ke_z,ke2d) = cvtot(:,ke)
        rhoe(:,ke_z,ke2d) = temp2(:,ke_z,ke2d) * cvtot(:,ke) * dens(:,ke)
        rhoe2(:,ke_z,ke2d) = rhoe(:,ke_z,ke2d)
 
        nz(:,ke_z,ke2d) = lcmesh%normal_fn(:,ke,3)
 
        flx_hydro(:,ke_z,ke2d) = 0.0_rp
      end do
      end do
      !$omp end do
      !$omp do collapse(3)
      do iq = this%vars%QS+1, this%vars%QE
      do ke2d = 1, lcmesh%Ne2D
      do ke_z = 1, lcmesh%NeZ
        ke = ke2d + (ke_z-1)*lcmesh%Ne2D
        rhoq(:,ke_z,ke2d,iq) = dens2(:,ke_z,ke2d) * qtrc(iq)%ptr%val(:,ke) &
                              + rhoq_t_mp(iq)%ptr%val(:,ke) * this%dtsec
        rhoq2(:,ke_z,ke2d,iq) = rhoq(:,ke_z,ke2d,iq)
      end do
      end do
      end do
      !$omp end do
      !$omp workshare
      sflx_rain(:,:) = 0.0_rp
      sflx_snow(:,:) = 0.0_rp
      sflx_engi(:,:) = 0.0_rp
      !$omp end workshare
      !$omp end parallel
 
      do iq = this%vars%QS+1, this%vars%QE
        if ( this%vars%vterm_hist_id(iq) > 0 ) then
          !$omp parallel do
          do ke = lcmesh%NeS, lcmesh%NeE
            this%vars%vterm_hist(iq)%local(domid)%val(:,ke) = 0.0_rp
          end do
        end if
      end do
 
      do step = 1, this%nstep_sedmientation
 
        !- Calculate terminal velocity 
 
        select case( this%MP_TYPEID )
        case( mp_typeid_kessler )
          do ke2d = 1, lcmesh%Ne2D
            !$omp parallel do private( &
            !$omp p2D, ke_z, iq, ColMask,                       &
            !$omp DENS2_tmp, REF_DENS_tmp, RHOQ2_tmp, vterm_tmp ) 
            do p2d = 1, elem3d%Nnode_h1D**2
              colmask(:) = elem3d%Colmask(:,p2d)
              do ke_z = 1, lcmesh%NeZ
                dens2_tmp(:,ke_z) = dens2(colmask(:),ke_z,ke2d)
                ref_dens_tmp(:,ke_z) = ref_dens(colmask(:),ke_z,ke2d)
              end do
              do iq = this%vars%QS+1, this%vars%QE
              do ke_z = 1, lcmesh%NeZ            
                rhoq2_tmp(:,ke_z,iq) = rhoq2(colmask(:),ke_z,ke2d,iq)
              end do
              end do
 
              call atmos_phy_mp_kessler_terminal_velocity( &
                elem3d%Nnode_v*lcmesh%NeZ, 1, elem3d%Nnode_v*lcmesh%NeZ,   & ! (in)
                dens2_tmp(:,:), rhoq2_tmp(:,:,:), ref_dens_tmp(:,:),       & ! (in)
                vterm_tmp(:,:,:)                                           ) ! (out)
 
              do iq = this%vars%QS+1, this%vars%QE
              do ke_z = 1, lcmesh%NeZ            
                vterm(colmask(:),ke_z,ke2d,iq) = vterm_tmp(:,ke_z,iq)
              end do
              end do
            end do
          end do
 
        end select   
        
        do iq = this%vars%QS+1, this%vars%QE
          if ( this%vars%vterm_hist_id(iq) > 0 ) then
            !$omp parallel do collapse(2) private(ke2D, ke_z, ke)
            do ke2d = 1, lcmesh%Ne2D
            do ke_z = 1, lcmesh%NeZ
              ke = ke2d + (ke_z-1)*lcmesh%Ne2D
              this%vars%vterm_hist(iq)%local(domid)%val(:,ke) = &
                this%vars%vterm_hist(iq)%local(domid)%val(:,ke) + vterm(:,ke_z,ke2d,iq) * this%rnstep_sedmientation
            end do
            end do
          end if
        end do 
 
        !- precipiation of hydrometers
 
        call atm_phy_mp_dgm_common_precipitation( &
          dens2, rhoq2, cptot2, cvtot2, rhoe2,                 & ! (inout)
          flx_hydro, sflx_rain, sflx_snow, sflx_engi,          & ! (inout)
          temp2, vterm,                                        & ! (in)
          this%dtsec_sedmientation, this%rnstep_sedmientation, & ! (in)
!          Dz, Lift, nz, vmapM, vmapP, IntWeight,               & ! (in)
          this%Dz, this%Lift, nz, vmapm, vmapp, intweight,               & ! (in)
          this%vars%QE - this%vars%QS, qla, qia,               & ! (in)
          lcmesh, elem3d )
 
        !$omp parallel do private(ke2D, ke_z, ke) collapse(2)
        do ke2d = 1, lcmesh%Ne2D
        do ke_z = 1, lcmesh%NeZ
          ke = ke2d + (ke_z-1)*lcmesh%Ne2D
          temp2(:,ke_z,ke2d) = rhoe2(:,ke_z,ke2d) / ( dens2(:,ke_z,ke2d) * cvtot2(:,ke_z,ke2d) )
        end do
        end do 
        
      end do ! end loop for step
 
      !$omp parallel private(ke2D, ke_z, iq, ke, CP_t, CV_t)
      !$omp workshare
      sflx_engi(:,:) = sflx_engi(:,:) - sflx_snow(:,:) * lhf ! moist internal energy
      !$omp end workshare
      !$omp do collapse(2)
      do ke2d = 1, lcmesh%Ne2D
      do ke_z = 1, lcmesh%NeZ
        ke = ke2d + (ke_z-1)*lcmesh%Ne2D
        dens_t_mp(:,ke) = ( dens2(:,ke_z,ke2d) - dens(:,ke) ) * rdt_mp
 
        cp_t(:) = ( cptot2(:,ke_z,ke2d) - cptot(:,ke) ) * rdt_mp
        cv_t(:) = ( cvtot2(:,ke_z,ke2d) - cvtot(:,ke) ) * rdt_mp
        rhoh_mp(:,ke) = rhoh_mp(:,ke) &
          + ( rhoe2(:,ke_z,ke2d) - rhoe(:,ke_z,ke2d) ) * rdt_mp &
          - ( cp_t(:) &
              + log( pres(:,ke) / pre00 ) * ( cvtot(:,ke) / cptot(:,ke) * cp_t(:) - cv_t(:) ) &
            ) * pres(:,ke) / rtot(:,ke) 
      end do
      end do
      !$omp end do
      !$omp do collapse(2)
      do iq = this%vars%QS+1, this%vars%QE
      do ke2d = 1, lcmesh%Ne2D
      do ke_z = 1, lcmesh%NeZ
        ke = ke2d + (ke_z-1)*lcmesh%Ne2D
        rhoq_t_mp(iq)%ptr%val(:,ke) = rhoq_t_mp(iq)%ptr%val(:,ke)       &
          + ( rhoq2(:,ke_z,ke2d,iq) - rhoq(:,ke_z,ke2d,iq) ) * rdt_mp
      end do
      end do
      end do
      !$omp end do
      !$omp end parallel
 
      !- precipiation of momentum
 
      call atm_phy_mp_dgm_common_precipitation_momentum( &
        rhou_t_mp, rhov_t_mp, momz_t_mp,          & ! (out)
        dens0, rhou2, rhov2, momz2, flx_hydro,    & ! (in)
!        Dz, Lift, nz, vmapM, vmapP,              & ! (in)
        this%Dz, this%Lift, nz, vmapm, vmapp,     & ! (in)
        lcmesh, elem3d                            ) ! (in)
 
    end if ! end if do_precipitation
 
    return
  end subroutine atmosphymp_calc_tendency_core
 
end module mod_atmos_phy_mp