mod_atmos_phy_sfc.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module Atmosphere / Physics / surface process
!!
!! @par Description
!!          Module for surface process
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module mod_atmos_phy_sfc
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !  
  use scale_precision
  use scale_prc
  use scale_io
  use scale_prof
  use scale_const, only: &
    undef8 => const_undef8  
  use scale_atmos_hydrometeor, only: &
    atmos_hydrometeor_dry
      
  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, elementbase2d, elementbase3d
    
  use scale_meshfield_base, only: meshfieldbase
  use scale_localmeshfield_base, only: localmeshfieldbase
 
  use scale_model_mesh_manager, only: modelmeshbase
  use scale_model_var_manager, only: modelvarmanager
  use scale_model_component_proc, only:  modelcomponentproc
  
  use mod_atmos_mesh, only: atmosmesh
  use mod_atmos_phy_sfc_vars, only: atmosphysfcvars
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
  !> Derived type to manage a component of surface process
  !!
  type, extends(modelcomponentproc), public :: atmosphysfc
    class(atmosmesh), pointer :: mesh    !< Pointer to a object to manage the mesh with atmospheric component
 
    integer :: atm_var_container_typeid !< Type ID of variable container for surface process
 
    integer :: sfcflx_typeid            !< Type id of surface scheme
    type(atmosphysfcvars) :: vars       !< An object to manage variables with surface component
 
    logical :: cpl_sw  !< Switch to indicate if the coupler is activated
  contains
    procedure, public :: setup => atmosphysfc_setup
    procedure, public :: set_coupler_flag => atmosphysfc_set_coupler_flag
    procedure, public :: calc_tendency => atmosphysfc_calc_tendency
    procedure, public :: update => atmosphysfc_update
    procedure, public :: finalize => atmosphysfc_finalize
  end type atmosphysfc
 
  !-----------------------------------------------------------------------------
  !++ Public parameters & variables
  !-----------------------------------------------------------------------------
 
  integer, parameter :: sfcflx_typeid_const           = 1 !< Type ID of constant surface flux scheme
  integer, parameter :: sfcflx_typeid_simple          = 2 !< Type ID of simple bulk surface flux scheme
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
 
  private :: cal_tend_from_sfcflx
  private :: cal_del_flux
  private :: convert_uv2localorthvec
  private :: convert_localorth2uvvec
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
    
contains
 
!> Setup an object to manage a component of surface process
!!
!! @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
!!
!OCL SERIAL
  subroutine atmosphysfc_setup( this, model_mesh, tm_parent_comp )
    use scale_time_manager, only: time_manager_component
 
    use scale_atmos_phy_sf_const, only: &
       atmos_phy_sf_const_setup
    use scale_atm_phy_sf_bulk_simple, only: &
       atmos_phy_sf_simple_setup
 
    use mod_atmos_vars, only: atm_vars_container_primary_id
 
    implicit none
    class(atmosphysfc), 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 surface process
    character(len=H_SHORT) :: time_dt_unit    = 'SEC'    !< Unit of timestep
    character(len=H_MID) :: sfcflx_type       = "CONST"  !< Type of surface flux scheme
    integer :: atm_var_container_typeid                  !< Type ID of variable container for surface process
    real(rp) :: default_sfc_temp              = 300.0_rp !< Default value of surface temperature
    namelist /param_atmos_phy_sfc/ &
      time_dt,                  &
      time_dt_unit,             &
      sfcflx_type,              &
      atm_var_container_typeid, &
      default_sfc_temp
    
    integer :: ierr
 
    integer :: ldomid
    class(localmesh2d), pointer :: lcmesh2d
    integer :: ke
    !--------------------------------------------------
 
    if (.not. this%IsActivated()) return
 
    log_newline
    log_info("ATMOS_PHY_SFC_setup",*) 'Setup'
 
    atm_var_container_typeid = atm_vars_container_primary_id
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_sfc,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("ATMOS_PHY_SFC_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
      log_error("ATMOS_PHY_SFC_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_SFC. Check!'
      call prc_abort
    endif
    log_nml(param_atmos_phy_sfc)
 
    this%atm_var_container_typeid = atm_var_container_typeid
 
    !--- Regist this compoent in the time manager
    
    call tm_parent_comp%Regist_process( 'ATMOS_PHY_SFC', time_dt, time_dt_unit, & ! (in)
      this%tm_process_id )                                                        ! (out)
 
    !- initialize the variables 
    call this%vars%Init( model_mesh, default_sfc_temp )      
 
    !--- Set the type of surface flux scheme
 
    select case( sfcflx_type )
    case ('CONST')
      call atmos_phy_sf_const_setup()
      this%SFCFLX_TYPEID = sfcflx_typeid_const
    case ('SIMPLE')
      call atmos_phy_sf_simple_setup()
      this%SFCFLX_TYPEID = sfcflx_typeid_simple
    case default
      log_error("ATMOS_PHY_SFC_setup",*) 'Not appropriate names of SFCFLX_TYPE in namelist PARAM_PHY_SFC. Check!'
      call prc_abort
    end select
 
    !--
    select type(model_mesh)
    class is (atmosmesh)
      this%mesh => model_mesh
    end select
 
    return
  end subroutine atmosphysfc_setup
 
  !> Set switch to indicate if the coupler is activated
  subroutine atmosphysfc_set_coupler_flag( this, cpl_sw )
    implicit none
    class(atmosphysfc), intent(inout) :: this
    logical, intent(in) :: cpl_sw
    !--------------------------------------------------
    this%CPL_sw = cpl_sw
    return
  end subroutine atmosphysfc_set_coupler_flag
 
  !> Calculate tendencies associated with a surface model
  !!
  !!
  !! @param model_mesh Object to manage computational mesh of atmospheric model 
  !! @param prgvars_list Object to mange prognostic variables with atmospheric dynamical core
  !! @param trcvars_list Object to mange auxiliary variables 
  !! @param forcing_list Object to mange forcing terms
  !! @param is_update Flag to specify whether the tendencies are updated in this call
  !!
!OCL SERIAL
  subroutine atmosphysfc_calc_tendency( &
    this, model_mesh, prgvars_list, trcvars_list, &
    auxvars_list, forcing_list, is_update         )
 
    use scale_tracer, only: &
      tracer_inq_id
 
    use mod_atmos_vars, only: &
      atmosvars_getlocalmeshprgvars,    &
      atmosvars_getlocalmeshphytends,   &
      atmosvars_getlocalmeshphyauxvars, &
      atmosvars_getlocalmeshqtrc_qv,    &
      atmosvars_getlocalmeshqtrcphytend
    use mod_atmos_phy_sfc_vars, only: &
      atmosphysfcvars_getlocalmeshfields
    implicit none
    
    class(atmosphysfc), 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, qv
    class(localmeshfieldbase), pointer :: dens_hyd, pres_hyd
    class(localmeshfieldbase), pointer :: rtot, cvtot, cptot
    class(localmeshfieldbase), pointer :: pres, pt
    class(localmeshfieldbase), pointer :: dens_tp, momx_tp, momy_tp, momz_tp, rhot_tp, rhoh_p, rhoqv_tp
    class(localmeshfieldbase), pointer :: sflx_mu, sflx_mv, sflx_mw, sflx_sh, sflx_lh, sflx_qv
    class(localmeshfieldbase), pointer :: sfc_temp
 
    integer :: n
    integer :: iq
    !--------------------------------------------------
    if (.not. this%IsActivated()) return
 
    call model_mesh%GetModelMesh( mesh )
 
    if (is_update) then
 
      do n=1, mesh%LOCAL_MESH_NUM
        call prof_rapstart( 'ATM_PHY_SFC_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_getlocalmeshphytends( n,        &
          mesh, forcing_list,                          &
          dens_tp, momx_tp, momy_tp, momz_tp, rhot_tp, &
          rhoh_p  )
        call atmosvars_getlocalmeshqtrc_qv( n, mesh, trcvars_list, forcing_list, &
          qv, rhoqv_tp )
        call atmosvars_getlocalmeshphyauxvars( n,      &
          mesh, auxvars_list,                          &
          pres, pt )
 
        call atmosphysfcvars_getlocalmeshfields( n,    &
          mesh, this%vars%SFCVARS_manager, this%vars%SFCFLX_manager, &
          sfc_temp, sflx_mu, sflx_mv, sflx_mw, sflx_sh, sflx_lh, sflx_qv )
        call prof_rapend( 'ATM_PHY_SFC_get_localmesh_ptr', 2)   
 
        call prof_rapstart( 'ATM_PHY_SFC_cal_tend', 2)
        call cal_tend_from_sfcflx( this, is_update,                                     &
          dens_tp%val, momx_tp%val, momy_tp%val, momz_tp%val, rhoh_p%val, rhoqv_tp%val, &
          sflx_mu%val, sflx_mv%val, sflx_mw%val, sflx_sh%val, sflx_lh%val, sflx_qv%val, &
          ddens%val, momx%val, momy%val, momz%val, drhot%val, qv%val,                   &
          dens_hyd%val, pres_hyd%val,                                                   &
          pres%val, pt%val, rtot%val, sfc_temp%val,                                     &
          model_mesh%DOptrMat(3), model_mesh%SOptrMat(3), model_mesh%LiftOptrMat,       &
          lcmesh, lcmesh%refElem3D, lcmesh%lcmesh2D, lcmesh%lcmesh2D%refElem2D          )
        call prof_rapend( 'ATM_PHY_SFC_cal_tend', 2)
      end do
 
    end if
 
    return  
  end subroutine atmosphysfc_calc_tendency
 
!> Update variables in a surface model
!!
!! @param model_mesh Object to manage computational mesh of atmospheric model 
!! @param prgvars_list Object to mange prognostic variables with atmospheric dynamical core
!! @param trcvars_list Object to mange auxiliary variables 
!! @param forcing_list Object to mange forcing terms
!! @param is_update Flag to speicfy whether the tendencies are updated in this call
!!
!OCL SERIAL
  subroutine atmosphysfc_update( this, model_mesh,           &
    prgvars_list, trcvars_list, auxvars_list,  forcing_list, &
    is_update )
    implicit none
 
    class(atmosphysfc), 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 atmosphysfc_update
 
!> Finalize a component of surface process
!!
!OCL SERIAL
  subroutine atmosphysfc_finalize( this )
    implicit none
    class(atmosphysfc), intent(inout) :: this
 
    !--------------------------------------------------
    if (.not. this%IsActivated()) return
 
    call this%vars%Final()
 
    return
  end subroutine atmosphysfc_finalize
 
!- private ------------------------------------------------
 
  !> Calculate tendencies associated with a surface model
!OCL SERIAL  
  subroutine cal_tend_from_sfcflx( this, is_update_sflx, &
    DENS_tp, MOMX_tp, MOMY_tp, MOMZ_tp, RHOH_p, RHOQV_tp, &
    SFLX_MU, SFLX_MV, SFLX_MW, SFLX_SH, SFLX_LH, SFLX_QV, &
    DDENS, MOMX, MOMY, MOMZ, DRHOT,                       &
    QV,                                                   &
    DENS_hyd, PRES_hyd, PRES, PT, Rtot,                   &
    SFC_TEMP,                                             & 
    Dz, Sz, Lift,                                         &
    lcmesh, elem, lcmesh2D, elem2D )
 
    use scale_const, only: &
      rdry => const_rdry,  &
      cpdry => const_cpdry, &
      cvdry => const_cvdry, &
      pres00 => const_pre00, &
      rplanet => const_radius    
    use scale_atmos_hydrometeor, only: &
      lhv, &
      i_qv
    use scale_tracer, only: &
      tracer_cv, qa
    use scale_atmos_phy_sf_const, only: &
       atmos_phy_sf_const_flux
    use scale_atm_phy_sf_bulk_simple, only: &
       atmos_phy_sf_simple_flux
       
    use scale_sparsemat, only: &
      sparsemat, &
      sparsemat_matmul    
    implicit none
 
    class(atmosphysfc), intent(inout) :: this
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem
    class(localmesh2d), intent(in) :: lcmesh2d
    class(elementbase2d), intent(in) :: elem2d
    logical, intent(in) :: is_update_sflx
    real(rp), intent(inout) :: dens_tp(elem%np,lcmesh%nea)
    real(rp), intent(inout) :: momx_tp(elem%np,lcmesh%nea)
    real(rp), intent(inout) :: momy_tp(elem%np,lcmesh%nea)
    real(rp), intent(inout) :: momz_tp(elem%np,lcmesh%nea)
    real(rp), intent(inout) :: rhoh_p (elem%np,lcmesh%nea)
    real(rp), intent(inout) :: rhoqv_tp(elem%np,lcmesh%nea)
    real(rp), intent(inout) :: sflx_mu(elem2d%np,lcmesh2d%nea)
    real(rp), intent(inout) :: sflx_mv(elem2d%np,lcmesh2d%nea)
    real(rp), intent(inout) :: sflx_mw(elem2d%np,lcmesh2d%nea)
    real(rp), intent(inout) :: sflx_sh(elem2d%np,lcmesh2d%nea)
    real(rp), intent(inout) :: sflx_lh(elem2d%np,lcmesh2d%nea)
    real(rp), intent(inout) :: sflx_qv(elem2d%np,lcmesh2d%nea)
    real(rp), intent(in) :: ddens(elem%np,lcmesh%nea)
    real(rp), intent(in) :: momx(elem%np,lcmesh%nea)
    real(rp), intent(in) :: momy(elem%np,lcmesh%nea)
    real(rp), intent(in) :: momz(elem%np,lcmesh%nea)
    real(rp), intent(in) :: drhot(elem%np,lcmesh%nea)
    real(rp), intent(in) :: qv(elem%np,lcmesh%nea)
    real(rp), intent(in) :: pres_hyd(elem%np,lcmesh%nea)
    real(rp), intent(in) :: dens_hyd(elem%np,lcmesh%nea)
    real(rp), intent(in) :: pres(elem%np,lcmesh%nea)
    real(rp), intent(in) :: pt(elem%np,lcmesh%nea)
    real(rp), intent(in) :: rtot(elem%np,lcmesh%nea)
    real(rp), intent(in) :: sfc_temp(elem2d%np,lcmesh2d%nea)
    type(sparsemat), intent(in) :: dz, sz, lift
 
    real(rp) :: atm_w   (elem2d%np,lcmesh2d%ne)
    real(rp) :: atm_u   (elem2d%np,lcmesh2d%ne)
    real(rp) :: atm_v   (elem2d%np,lcmesh2d%ne)
    real(rp) :: atm_temp(elem2d%np,lcmesh2d%ne)
    real(rp) :: atm_pres(elem2d%np,lcmesh2d%ne)
    real(rp) :: atm_qv  (elem2d%np,lcmesh2d%ne)
    real(rp) :: dz1     (elem2d%np,lcmesh2d%ne)
    real(rp) :: z1      (elem2d%np,lcmesh2d%ne)
    real(rp) :: sfc_pres(elem2d%np,lcmesh2d%ne)
    real(rp) :: sfc_dens(elem2d%np,lcmesh2d%ne)
 
    real(rp) :: sflx_qtrc(elem2d%np,lcmesh2d%ne,max(1,qa))
    real(rp) :: sflx_engi(elem2d%np,lcmesh2d%ne)
 
    ! dummy
    real(rp) :: u10(elem2d%np,lcmesh2d%ne)
    real(rp) :: v10(elem2d%np,lcmesh2d%ne)
 
    integer :: ke
    integer :: ke2d
    integer :: hslicez0, hslicez1
    integer :: ij
    integer :: iq
    real(rp) :: dens
    real(rp) :: liftdelflx(elem%np)
    real(rp) :: del_flux(elem%nfptot,lcmesh%ne,5)
 
    real(rp) :: sflx_mu1, sflx_mv1
 
    !-------------------------------------------------
 
    if ( is_update_sflx .and. (.not. this%CPL_sw) ) then
      !$omp parallel do collapse(2) private( &
      !$omp ke, hSliceZ0, hsliceZ1,          &
      !$omp dens                             )
      do ke2d=lcmesh2d%NeS, lcmesh2d%NeE
      do ij=1, elem2d%Np
 
        ke = ke2d
        hslicez0 = elem%Hslice(ij,1)
        hslicez1 = elem%Hslice(ij,2)
 
        dens = dens_hyd(hslicez1,ke) + ddens(hslicez1,ke)
        atm_u(ij,ke2d) = momx(hslicez1,ke) / dens
        atm_v(ij,ke2d) = momy(hslicez1,ke) / dens
        atm_w(ij,ke2d) = 0.0_rp
        atm_pres(ij,ke2d) = pres(hslicez1,ke)
!        ATM_TEMP(ij,ke2D) = ATM_PRES(ij,ke2D) / ( Rtot(hsliceZ1,ke) * dens )
 
        atm_qv(ij,ke2d) = qv(hslicez1,ke)
 
        sfc_dens(ij,ke2d) = dens_hyd(hslicez0,ke) + ddens(hslicez0,ke)
        sfc_pres(ij,ke2d) = pres(hslicez0,ke)
        ! Tentative: 
        atm_temp(ij,ke2d) = sfc_pres(ij,ke2d) / ( rtot(hslicez0,ke) * sfc_dens(ij,ke2d) )
        atm_qv(ij,ke2d) = qv(hslicez0,ke)
 
        z1(ij,ke2d) = lcmesh%zlev(hslicez1,ke)
        dz1(ij,ke2d) = z1(ij,ke2d) - lcmesh%zlev(hslicez0,ke)
        z1(ij,ke2d) = z1(ij,ke2d) + rplanet
      end do
      end do
 
      call convert_uv2localorthvec( &
        this%mesh%ptr_mesh, lcmesh2d%pos_en(:,:,1), lcmesh2d%pos_en(:,:,2), z1(:,:), elem2d%Np*lcmesh2d%Ne, &
        atm_u, atm_v ) ! (inout)
      
      !---
      select case ( this%SFCFLX_TYPEID )
      case ( sfcflx_typeid_const )
        call atmos_phy_sf_const_flux( &
          elem2d%Np, 1, elem2d%Np, lcmesh2d%NeA, 1, lcmesh2d%Ne, & ! (in)
          atm_w(:,:), atm_u(:,:), atm_v(:,:), sfc_temp(:,:),     & ! (in)
          dz1(:,:), sfc_dens(:,:),                               & ! (in)
          sflx_mw(:,:), sflx_mu(:,:), sflx_mv(:,:),              & ! (out)
          sflx_sh(:,:), sflx_lh(:,:), sflx_qv(:,:),              & ! (out)
          u10(:,:), v10(:,:)                                     ) ! (out)
      case ( sfcflx_typeid_simple )
        call atmos_phy_sf_simple_flux( &
          elem2d%Np, 1, elem2d%Np, lcmesh2d%NeA, 1, lcmesh2d%Ne,            & ! (in)
          atm_w(:,:), atm_u(:,:), atm_v(:,:), atm_temp(:,:), atm_pres(:,:), & ! (in) Note: ATM_PRES is not used
          atm_qv(:,:),                                                      & ! (in)
          sfc_dens(:,:), sfc_temp(:,:), sfc_pres(:,:),                      & ! (in)
          dz1(:,:),                                                         & ! (in)
          sflx_mw(:,:), sflx_mu(:,:), sflx_mv(:,:),                         & ! (out)
          sflx_sh(:,:), sflx_lh(:,:), sflx_qv(:,:),                         & ! (out)
          u10(:,:), v10(:,:)                                                ) ! (out)
      end select
 
      call convert_localorth2uvvec( &
        this%mesh%ptr_mesh, lcmesh2d%pos_en(:,:,1), lcmesh2d%pos_en(:,:,2), z1(:,:), elem2d%Np*lcmesh2d%Ne, &
        sflx_mu(:,lcmesh2d%NeS:lcmesh2d%NeE), sflx_mv(:,lcmesh2d%NeS:lcmesh2d%NeE) ) ! (inout)
 
    end if
 
    !- Add the tendency due to the surface flux
 
    call cal_del_flux( del_flux,            &
      sflx_mu, sflx_mv, sflx_mw, sflx_sh,   &
      sflx_qv,                              &
      sfc_temp,                             &
      lcmesh%normal_fn(:,:,3),              &
      lcmesh, elem, lcmesh2d, elem2d        )
    
    !$omp parallel do private(ke, iq, LiftDelFlx)
    do ke2d=1, lcmesh2d%Ne
      ke = ke2d
 
      call sparsemat_matmul( lift, lcmesh%Fscale(:,ke)*del_flux(:,ke,1), liftdelflx )
      momx_tp(:,ke) = momx_tp(:,ke) - liftdelflx(:) 
 
      call sparsemat_matmul( lift, lcmesh%Fscale(:,ke)*del_flux(:,ke,2), liftdelflx )
      momy_tp(:,ke) = momy_tp(:,ke) - liftdelflx(:) 
 
      call sparsemat_matmul( lift, lcmesh%Fscale(:,ke)*del_flux(:,ke,3), liftdelflx )
      momz_tp(:,ke) = momz_tp(:,ke) - liftdelflx(:) 
      
      call sparsemat_matmul( lift, lcmesh%Fscale(:,ke)*del_flux(:,ke,4), liftdelflx )
      rhoh_p(:,ke) = rhoh_p(:,ke) - liftdelflx(:) 
 
      if ( .not. atmos_hydrometeor_dry ) then
        sflx_qtrc(:,ke2d,i_qv) = sflx_qv(:,ke2d)
        sflx_engi(:,ke2d) = sflx_qv(:,ke2d) * ( tracer_cv(i_qv) * sfc_temp(:,ke2d) + lhv )
 
        call sparsemat_matmul( lift, lcmesh%Fscale(:,ke)*del_flux(:,ke,5), liftdelflx )
        rhoqv_tp(:,ke) = rhoqv_tp(:,ke) - liftdelflx(:)
        dens_tp(:,ke) = dens_tp(:,ke) - liftdelflx(:)
      end if 
    end do
 
    return
  end subroutine cal_tend_from_sfcflx
 
!OCL SERIAL  
  subroutine cal_del_flux( del_flux,      &
    sflx_mu, sflx_mv, sflx_mw, sflx_sh,   &
    sflx_qv,                              &
    SFC_TEMP,                             &
    nz,                                   &
    lmesh, elem, lmesh2D, elem2D          )
 
    use scale_atmos_hydrometeor, only: &
      lhv, &
      i_qv
    use scale_tracer, only: &
      tracer_cv, qa    
    implicit none
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem  
    class(localmesh2d), intent(in) :: lmesh2d
    class(elementbase2d), intent(in) :: elem2d   
    real(rp), intent(out) ::  del_flux(elem%nfptot*lmesh%ne,5)
    real(rp), intent(in) :: sflx_mu(elem2d%np,lmesh2d%nea)
    real(rp), intent(in) :: sflx_mv(elem2d%np,lmesh2d%nea)
    real(rp), intent(in) :: sflx_mw(elem2d%np,lmesh2d%nea)
    real(rp), intent(in) :: sflx_sh(elem2d%np,lmesh2d%nea)
    real(rp), intent(in) :: sflx_qv(elem2d%np,lmesh2d%nea)
    real(rp), intent(in) :: sfc_temp(elem2d%np,lmesh2d%nea)
    real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
 
    integer :: ke2d, p
    integer :: i
    !------------------------------------------------------------------------
 
    !$omp parallel private(i)
    !$omp workshare
    del_flux(:,:) = 0.0_rp
    !$omp end workshare
 
    !$omp do collapse(2)
    do ke2d=1, lmesh2d%Ne
    do p=1, elem2d%Np
      i = elem%Nfaces_h*elem%Nfp_h + p + (ke2d-1)*elem%NfpTot
      del_flux(i,1) = sflx_mu(p,ke2d) * nz(i)
      del_flux(i,2) = sflx_mv(p,ke2d) * nz(i)
      del_flux(i,3) = sflx_mw(p,ke2d) * nz(i)
      del_flux(i,4) = sflx_sh(p,ke2d) * nz(i)
      del_flux(i,5) = sflx_qv(p,ke2d) * nz(i)
    end do
    end do
    !$omp end do
 
    if ( .not. atmos_hydrometeor_dry ) then
      !$omp do collapse(2)
      do ke2d=1, lmesh2d%Ne
      do p=1, elem2d%Np
        i = elem%Nfaces_h*elem%Nfp_h + p + (ke2d-1)*elem%NfpTot
        del_flux(i,4) = del_flux(i,4) &
          + sflx_qv(p,ke2d) * ( tracer_cv(i_qv) * sfc_temp(p,ke2d) + 0.0_rp * lhv ) * nz(i)
      end do
      end do
    end if     
    !$omp end parallel
 
    return
  end subroutine cal_del_flux
 
!OCL SERIAL
  subroutine convert_uv2localorthvec( mesh, x1, x2, r, N, & ! (in)
    u, v ) ! (inout)
    use scale_cubedsphere_coord_cnv, only: &
      cubedspherecoordcnv_cs2localorthvec_alpha
    use scale_mesh_cubedspheredom3d, only: meshcubedspheredom3d
 
    implicit none
    integer, intent(in) :: n
    class(meshbase3d), intent(in) :: mesh
    real(rp), intent(in) :: x1(n)
    real(rp), intent(in) :: x2(n)
    real(rp), intent(in) :: r(n)
    real(rp), intent(inout) :: u(n)
    real(rp), intent(inout) :: v(n)
 
    select type(mesh)
    type is (meshcubedspheredom3d)
      call cubedspherecoordcnv_cs2localorthvec_alpha( &
        x1, x2, r, n, & ! (in)
        u, v          ) ! (inout)
    end select
    return
  end subroutine convert_uv2localorthvec
 
!OCL SERIAL
  subroutine convert_localorth2uvvec( mesh, x1, x2, r, N, & ! (in)
    u, v ) ! (inout)
    use scale_cubedsphere_coord_cnv, only: &
      cubedspherecoordcnv_localorth2csvec_alpha
    use scale_mesh_cubedspheredom3d, only: meshcubedspheredom3d
 
    implicit none
    integer, intent(in) :: n
    class(meshbase3d), intent(in) :: mesh
    real(rp), intent(in) :: x1(n)
    real(rp), intent(in) :: x2(n)
    real(rp), intent(in) :: r(n)
    real(rp), intent(inout) :: u(n)
    real(rp), intent(inout) :: v(n)
 
    select type(mesh)
    type is (meshcubedspheredom3d)
      call cubedspherecoordcnv_localorth2csvec_alpha( &
        x1, x2, r, n, & ! (in)
        u, v          ) ! (inout)
    end select
    return
  end subroutine convert_localorth2uvvec 
end module mod_atmos_phy_sfc