mod_atmos_phy_mp.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module Atmosphere / Physics / cloud microphysics
!!
!! @par Description
!!          Module for cloud microphysics
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#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
  !
 
  !> Derived type to manage a component of cloud microphysics
  !!
  type, extends(modelcomponentproc), public :: atmosphymp
    integer :: mp_typeid         !< Type id of cloud microphysics scheme
    type(atmosphympvars) :: vars !< Object to manage variables with cloud microphysics
 
    logical, private :: do_precipitation    !< Apply sedimentation (precipitation)?
    logical, private :: do_negative_fixer   !< Apply negative fixer?
    real(rp), private :: limit_negative     !< Abort if abs(fixed negative value) > abs(MP_limit_negative)
    integer, private :: ntmax_sedimentation !< Number of time step for sedimentation
    real(rp), private :: max_term_vel       !< Terminal velocity for calculate dt of sedimentation
    real(rp), private :: cldfrac_thleshold  !< Threshold for cloud fraction
 
    real(dp), private :: dtsec                !< Timestep with cloud microphysics component
    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
  integer, parameter :: mp_typeid_tomita08 = 2
  integer, parameter :: mp_typeid_sn14     = 3
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
    
  logical :: flag_lt = .false. ! Tentative
 
contains
 
!> Setup a component of cloud microphysics
!!
!! @param model_mesh Object to manage computational mesh of atmospheric model 
!! @param tm_parent_comp Object to mange a temporal scheme in a parent component
!!
  subroutine atmosphymp_setup( this, model_mesh, tm_parent_comp )
    use scale_const, only: &
      eps => const_eps
    use scale_tracer, only: &
      tracer_regist
    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_atmos_phy_mp_tomita08, only: &
      atmos_phy_mp_tomita08_setup,               &
      atmos_phy_mp_tomita08_ntracers,            &
      atmos_phy_mp_tomita08_nwaters,             &
      atmos_phy_mp_tomita08_nices,               &
      atmos_phy_mp_tomita08_tracer_names,        &
      atmos_phy_mp_tomita08_tracer_descriptions, &
      atmos_phy_mp_tomita08_tracer_units
    use scale_atmos_phy_mp_sn14, only: &
      atmos_phy_mp_sn14_setup,               &
      atmos_phy_mp_sn14_ntracers,            &
      atmos_phy_mp_sn14_nwaters,             &
      atmos_phy_mp_sn14_nices,               &
      atmos_phy_mp_sn14_tracer_names,        &
      atmos_phy_mp_sn14_tracer_descriptions, &
      atmos_phy_mp_sn14_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 !< Timestep for cloud microcloud physics
    character(len=H_SHORT) :: time_dt_unit          = 'SEC'  !< Unit of timestep
 
    character(len=H_MID) :: mp_type = 'KESSLER'              !< Type of cloud microcloud physics scheme
 
    logical :: do_precipitation     !< Flag whether sedimentation (precipitation) is applied
    logical :: do_negative_fixer    !< Flag whether negative fixer is applied
    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
    class(localmesh3d), pointer :: lcmesh3d
    class(elementbase3d), pointer :: elem3d
 
    integer :: ierr
 
    integer :: qs_mp, qe_mp, qa_mp
    integer :: qs2
 
    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
 
    !--- Register this component 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 ( 'TOMITA08' )
      this%MP_TYPEID = mp_typeid_tomita08
      call atmos_hydrometeor_regist( &
            atmos_phy_mp_tomita08_nwaters,                & ! (in)
            atmos_phy_mp_tomita08_nices,                  & ! (in)
            atmos_phy_mp_tomita08_tracer_names(:),        & ! (in)
            atmos_phy_mp_tomita08_tracer_descriptions(:), & ! (in)
            atmos_phy_mp_tomita08_tracer_units(:),        & ! (in)
            qs_mp                                         ) ! (out)
      qa_mp = atmos_phy_mp_tomita08_ntracers
    ! case( 'SN14' )
    !   this%MP_TYPEID = MP_TYPEID_SN14
    !   call ATMOS_HYDROMETEOR_regist( &
    !         ATMOS_PHY_MP_SN14_nwaters,                   & ! (in)
    !         ATMOS_PHY_MP_SN14_nices,                     & ! (in)
    !         ATMOS_PHY_MP_SN14_tracer_names(1:6),         & ! (in)
    !         ATMOS_PHY_MP_SN14_tracer_descriptions(1:6),  & ! (in)
    !         ATMOS_PHY_MP_SN14_tracer_units(1:6),         & ! (in)
    !         QS_MP                                        ) ! (out)
      
    !   call TRACER_regist( QS2,                           & ! (out)
    !         5,                                           & ! (in)
    !         ATMOS_PHY_MP_SN14_tracer_names(7:11),        & ! (in)
    !         ATMOS_PHY_MP_SN14_tracer_descriptions(7:11), & ! (in)
    !         ATMOS_PHY_MP_SN14_tracer_units(7:11)         ) ! (in)
    !   QA_MP = ATMOS_PHY_MP_SN14_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()
    case( mp_typeid_tomita08 )
      lcmesh3d => atm_mesh%ptr_mesh%lcmesh_list(1)
      elem3d => lcmesh3d%refElem3D
      call atmos_phy_mp_tomita08_setup( &
        elem3d%Nnode_v*lcmesh3d%NeZ, 1, elem3d%Nnode_v*lcmesh3d%NeZ, elem3d%Nnode_h1D**2, 1, elem3d%Nnode_h1D**2, lcmesh3d%Ne2D, 1, lcmesh3d%Ne2D, &
        flag_lt )    
    case( mp_typeid_sn14 )
    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
 
!> Calculate tendencies associated with cloud microphysics
!!
!!
!! @param model_mesh a object to manage computational mesh of atmospheric model 
!! @param prgvars_list a object to mange prognostic variables with atmospheric dynamical core
!! @param trcvars_list a object to mange auxiliary variables 
!! @param forcing_list a object to mange forcing terms
!! @param is_update Flag to speicfy whether the tendencies are updated in this call
!!
!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
 
!> Update variables in a component of cloud microphysics
!!
!! @param model_mesh a object to manage computational mesh of atmospheric model 
!! @param prgvars_list a object to mange prognostic variables with atmospheric dynamical core
!! @param trcvars_list a object to mange auxiliary variables 
!! @param forcing_list a object to mange forcing terms
!! @param is_update Flag to speicfy whether the tendencies are updated in this call
!!
!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
 
!> Finalize a component of cloud microphysics
!!
!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_terminal_velocity
    use scale_atmos_phy_mp_tomita08, only: &
      atmos_phy_mp_tomita08_terminal_velocity
    use scale_atmos_phy_mp_sn14, only: &
      atmos_phy_mp_sn14_tendency, &
      atmos_phy_mp_sn14_terminal_velocity
    use scale_sparsemat, only: sparsemat
    use scale_atm_phy_mp_dgm_common, only:  &
      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) :: dens (elem3d%np,lcmesh%nea)
 
    real(rp) :: rhoe_t(elem3d%np,lcmesh%nea)
 
    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) :: temp2_tmp(elem3d%nnode_v,lcmesh%nez)
    real(rp) :: pres2_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 lcmesh%GetVmapZ1D( vmapm, vmapp ) ! (out)
    
    call atm_phy_mp_dgm_common_gen_intweight( intweight, & ! (out)
      lcmesh                                             ) ! (in)
    
    !$omp parallel do private(ke)
    do ke = lcmesh%NeS, lcmesh%NeE
      dens(:,ke) = dens_hyd(:,ke) + ddens(:,ke)
    end do
 
    select case( this%MP_TYPEID )
    case( mp_typeid_kessler )
      call calc_tendency_kessler( this, &
        rhoq_t_mp, cptot_t, cvtot_t, rhoe_t, evaporate,  & ! (out)
        dens, qtrc, pres, dens_hyd, rtot, cvtot, cptot,  & ! (in)
        rdt_mp, lcmesh, elem3d )                           ! (in)
 
    case( mp_typeid_tomita08 )
      call calc_tendency_tomita08( this, &
        rhoq_t_mp, cptot_t, cvtot_t, rhoe_t, evaporate,  & ! (out)
        dens, qtrc, pres, dens_hyd, rtot, cvtot, cptot,  & ! (in)
        rdt_mp, lcmesh, elem3d )                           ! (in)
 
    case( mp_typeid_sn14 )
      call calc_tendency_sn14( this, &
        rhoq_t_mp, cptot_t, cvtot_t, rhoe_t, evaporate,  & ! (out)
        dens, qtrc, pres, dens_hyd, rtot, cvtot, cptot,  & ! (in)
        rdt_mp, lcmesh, elem3d )                           ! (in)
 
    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 
 
          do ke2d = 1, lcmesh%Ne2D
            !$omp parallel do private( &
            !$omp p2D, ke_z, iq, ColMask,                                   &
            !$omp DENS2_tmp, REF_DENS_tmp, RHOQ2_tmp, TEMP2_tmp, PRES2_tmp, &
            !$omp 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)
                temp2_tmp(:,ke_z) = temp2(colmask(:),ke_z,ke2d)
                pres2_tmp(:,ke_z) = pres2(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
 
              select case( this%MP_TYPEID )
              case( mp_typeid_kessler )
                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)
              case( mp_typeid_tomita08 )
                call atmos_phy_mp_tomita08_terminal_velocity( &
                  elem3d%Nnode_v*lcmesh%NeZ, 1, elem3d%Nnode_v*lcmesh%NeZ, & ! (in)
                  dens2_tmp(:,:), temp2_tmp(:,:), rhoq2_tmp(:,:,:),        & ! (in)
                  vterm_tmp(:,:,:)                                         ) ! (out)
              case( mp_typeid_sn14 )
                call atmos_phy_mp_sn14_terminal_velocity( &
                  elem3d%Nnode_v*lcmesh%NeZ, 1, elem3d%Nnode_v*lcmesh%NeZ,          & ! (in)
                  dens2_tmp(:,:), temp2_tmp(:,:), rhoq2_tmp(:,:,:), pres2_tmp(:,:), & ! (in)
                  vterm_tmp(:,:,:)                                                  ) ! (out)
              end select   
 
              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
 
        
        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
 
!OCL SERIAL
  subroutine calc_tendency_kessler( this, &
    RHOQ_t_MP, CPtot_t, CVtot_t, RHOE_t, EVAPORATE,  & ! (out)
    dens, qtrc, pres, dens_hyd, rtot, cvtot, cptot,  & ! (in)
    rdt_mp, lcmesh, elem3d )                           ! (in)
 
    use scale_atmos_phy_mp_kessler, only:    &
      atmos_phy_mp_kessler_adjustment
    implicit none
 
    class(atmosphymp), intent(in) :: this
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem3d
    type(localmeshfieldbaselist), intent(inout) :: rhoq_t_mp(this%vars%qs:this%vars%qe)
    real(rp), intent(out) :: cptot_t(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: cvtot_t(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: rhoe_t(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: evaporate(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: dens(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) :: 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)
    real(rp), intent(in) :: rdt_mp
 
    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)
 
    integer :: ke
    integer :: iq
    !------------------------------------------------------------
 
    !$omp parallel private(ke, iq)
    !$omp do
    do ke = lcmesh%NeS, lcmesh%NeE
      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
    return
  end subroutine calc_tendency_kessler
 
!OCL SERIAL
  subroutine calc_tendency_tomita08( this, &
    RHOQ_t_MP, CPtot_t, CVtot_t, RHOE_t, EVAPORATE,  & ! (out)
    dens, qtrc, pres, dens_hyd, rtot, cvtot, cptot,  & ! (in)
    rdt_mp, lcmesh, elem3d )                           ! (in)
 
    use scale_const, only: &
      undef => const_undef
    use scale_atmos_phy_mp_tomita08, only: &
      atmos_phy_mp_tomita08_adjustment
    implicit none
 
    class(atmosphymp), intent(in) :: this
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem3d
    type(localmeshfieldbaselist), intent(inout) :: rhoq_t_mp(this%vars%qs:this%vars%qe)
    real(rp), intent(out) :: cptot_t(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: cvtot_t(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: rhoe_t(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: evaporate(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: dens(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) :: 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)
    real(rp), intent(in) :: rdt_mp
 
    real(rp) :: dens_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: pres_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: ccn_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: temp1_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: cptot1_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: cvtot1_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: qtrc1_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d,this%vars%qs:this%vars%qe)
    real(rp) :: rhoe_t_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: evaporate_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
 
    integer :: ke, ke_xy, ke_z
    integer :: p,p_xy, p_z
    integer :: iq
    !------------------------------------------------------------
 
    !$omp parallel do collapse(2) private(p_z,p_xy,iq, ke,p)
    do ke_z=1, lcmesh%NeZ
    do ke_xy=1, lcmesh%Ne2D
      ke = ke_xy + (ke_z-1)*lcmesh%Ne2D
      do p_z=1, elem3d%Nnode_v
      do p_xy=1, elem3d%Nnode_h1D**2
        p = p_xy + (p_z-1)*elem3d%Nnode_h1D**2
        dens_z(p_z,ke_z,p_xy,ke_xy) = dens(p,ke)
        pres_z(p_z,ke_z,p_xy,ke_xy) = pres(p,ke)
        temp1_z(p_z,ke_z,p_xy,ke_xy) = pres(p,ke) / ( dens(p,ke) * rtot(p,ke) )
        cptot1_z(p_z,ke_z,p_xy,ke_xy) = cptot(p,ke)
        cvtot1_z(p_z,ke_z,p_xy,ke_xy) = cvtot(p,ke)
 
        ccn_z(p_z,ke_z,p_xy,ke_xy) = undef
      end do
      end do
      do iq = this%vars%QS, this%vars%QE
        do p_z=1, elem3d%Nnode_v
        do p_xy=1, elem3d%Nnode_h1D**2
          p = p_xy + (p_z-1)*elem3d%Nnode_h1D**2
          qtrc1_z(p_z,ke_z,p_xy,ke_xy,iq) = qtrc(iq)%ptr%val(p,ke)
        end do
        end do
      end do
    end do
    end do
 
    call atmos_phy_mp_tomita08_adjustment( &
      elem3d%Nnode_v*lcmesh%NeZ, 1, elem3d%Nnode_v*lcmesh%NeZ, elem3d%Nnode_h1D**2, 1, elem3d%Nnode_h1D**2, lcmesh%Ne2D, 1, lcmesh%Ne2D, & ! (in)
      dens_z, pres_z, ccn_z, this%dtsec,       & ! (in)
      temp1_z, qtrc1_z, cptot1_z, cvtot1_z,    & ! (inout)
      rhoe_t_z, evaporate_z                    ) ! (out)
  
    !$omp parallel do collapse(2) private(p_z,p_xy,iq, ke,p)
    do ke_z=1, lcmesh%NeZ
    do ke_xy=1, lcmesh%Ne2D
      ke = ke_xy + (ke_z-1)*lcmesh%Ne2D
      do p_z=1, elem3d%Nnode_v
      do p_xy=1, elem3d%Nnode_h1D**2
        p = p_xy + (p_z-1)*elem3d%Nnode_h1D**2
        cptot_t(p,ke) = ( cptot1_z(p_z,ke_z,p_xy,ke_xy) - cptot(p,ke) ) * rdt_mp
        cvtot_t(p,ke) = ( cvtot1_z(p_z,ke_z,p_xy,ke_xy) - cvtot(p,ke) ) * rdt_mp
        rhoe_t(p,ke) = rhoe_t_z(p_z,ke_z,p_xy,ke_xy)
        evaporate(p,ke) = evaporate_z(p_z,ke_z,p_xy,ke_xy)
      end do
      end do
      do iq = this%vars%QS, this%vars%QE
        do p_z=1, elem3d%Nnode_v
        do p_xy=1, elem3d%Nnode_h1D**2
          p = p_xy + (p_z-1)*elem3d%Nnode_h1D**2
          rhoq_t_mp(iq)%ptr%val(p,ke) = ( qtrc1_z(p_z,ke_z,p_xy,ke_xy,iq) - qtrc(iq)%ptr%val(p,ke) ) * dens(p,ke) * rdt_mp
        end do
        end do
      end do
    end do
    end do
 
    return
  end subroutine calc_tendency_tomita08
 
!OCL SERIAL
  subroutine calc_tendency_sn14( this, &
    RHOQ_t_MP, CPtot_t, CVtot_t, RHOE_t, EVAPORATE,  & ! (out)
    dens, qtrc, pres, dens_hyd, rtot, cvtot, cptot,  & ! (in)
    rdt_mp, lcmesh, elem3d )                           ! (in)
 
    use scale_const, only: &
      undef => const_undef
    use scale_atmos_phy_mp_sn14, only: &
       atmos_phy_mp_sn14_tendency
    implicit none
 
    class(atmosphymp), intent(in) :: this
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem3d
    type(localmeshfieldbaselist), intent(inout) :: rhoq_t_mp(this%vars%qs:this%vars%qe)
    real(rp), intent(out) :: cptot_t(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: cvtot_t(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: rhoe_t(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: evaporate(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: dens(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) :: 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)
    real(rp), intent(in) :: rdt_mp
 
    real(rp) :: dens_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: pres_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: ccn_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: temp1_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: cptot1_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: cvtot1_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: qtrc1_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d,this%vars%qs:this%vars%qe)
    real(rp) :: rhoq_t_mp_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d,this%vars%qs:this%vars%qe)
    real(rp) :: rhoe_t_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: evaporate_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: cptot_t_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp) :: cvtot_t_z(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
 
    integer :: ke, ke_xy, ke_z
    integer :: p,p_xy, p_z
    integer :: iq
    !------------------------------------------------------------
 
    !$omp parallel do collapse(2) private(p_z,p_xy,iq, ke,p)
    do ke_z=1, lcmesh%NeZ
    do ke_xy=1, lcmesh%Ne2D
      ke = ke_xy + (ke_z-1)*lcmesh%Ne2D
      do p_z=1, elem3d%Nnode_v
      do p_xy=1, elem3d%Nnode_h1D**2
        p = p_xy + (p_z-1)*elem3d%Nnode_h1D**2
        dens_z(p_z,ke_z,p_xy,ke_xy) = dens(p,ke)
        pres_z(p_z,ke_z,p_xy,ke_xy) = pres(p,ke)
        temp1_z(p_z,ke_z,p_xy,ke_xy) = pres(p,ke) / ( dens(p,ke) * rtot(p,ke) )
        cptot1_z(p_z,ke_z,p_xy,ke_xy) = cptot(p,ke)
        cvtot1_z(p_z,ke_z,p_xy,ke_xy) = cvtot(p,ke)
 
        ccn_z(p_z,ke_z,p_xy,ke_xy) = undef
      end do
      end do
      do iq = this%vars%QS, this%vars%QE
        do p_z=1, elem3d%Nnode_v
        do p_xy=1, elem3d%Nnode_h1D**2
          p = p_xy + (p_z-1)*elem3d%Nnode_h1D**2
          qtrc1_z(p_z,ke_z,p_xy,ke_xy,iq) = qtrc(iq)%ptr%val(p,ke)
        end do
        end do
      end do
    end do
    end do
 
    ! call  ATMOS_PHY_MP_sn14_tendency( &
    !   elem3D%Nnode_v*lcmesh%NeZ, 1, elem3D%Nnode_v*lcmesh%NeZ, elem3D%Nnode_h1D**2, 1, elem3D%Nnode_h1D**2, lcmesh%Ne2D, 1, lcmesh%Ne2D, & ! (in)
    !   DENS_z, W_z, QTRC1_z, PRES_z, TEMP_z, Qdry_z, CPtot1_z, CVtot1_z, CCN_z, this%dtsec, REAL_CZ, REAL_FZ,  & ! (in)
    !   RHOQ_t_MP_z, RHOE_t_z, CPtot_t_z, CVtot_t_z, EVAPORATE_z )                                                ! (out)
  
 
    !$omp parallel do collapse(2) private(p_z,p_xy,iq, ke,p)
    do ke_z=1, lcmesh%NeZ
    do ke_xy=1, lcmesh%Ne2D
      ke = ke_xy + (ke_z-1)*lcmesh%Ne2D
      do p_z=1, elem3d%Nnode_v
      do p_xy=1, elem3d%Nnode_h1D**2
        p = p_xy + (p_z-1)*elem3d%Nnode_h1D**2
        cptot_t(p,ke) = cptot_t_z(p_z,ke_z,p_xy,ke_xy)
        cvtot_t(p,ke) = cvtot_t_z(p_z,ke_z,p_xy,ke_xy)
      end do
      end do
      do iq = this%vars%QS, this%vars%QE
        do p_z=1, elem3d%Nnode_v
        do p_xy=1, elem3d%Nnode_h1D**2
          p = p_xy + (p_z-1)*elem3d%Nnode_h1D**2
          rhoq_t_mp(iq)%ptr%val(p,ke) = rhoq_t_mp_z(p_z,ke_z,p_xy,ke_xy,iq)
        end do
        end do
      end do
    end do
    end do
 
    return
  end subroutine calc_tendency_sn14  
end module mod_atmos_phy_mp