mod_atmos_dyn.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module mod_atmos_dyn
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !  
  use scale_precision
  use scale_prc
  use scale_io
  use scale_prof
  use scale_const, only: &
    undef8 => const_undef8
 
  use scale_sparsemat, only: sparsemat
  use scale_timeint_rk, only: timeint_rk
  
  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, localmeshfieldbaselist
 
  use scale_model_mesh_manager, only: modelmeshbase
  use scale_model_var_manager, only: modelvarmanager
  use scale_model_component_proc, only:  modelcomponentproc
 
  use scale_atm_dyn_dgm_driver_nonhydro3d, only: atmdyndgmdriver_nonhydro3d
  use scale_atm_dyn_dgm_driver_trcadv3d, only: atmdyndgmdriver_trcadv3d
 
  use scale_atm_dyn_dgm_nonhydro3d_numdiff, only: atmdyn_nonhydro3d_numdiff
  
  use mod_atmos_mesh, only: atmosmesh
  use mod_atmos_vars, only: &
    atmosvars_getlocalmeshprgvar,        &
    atmosvars_getlocalmeshprgvars,       &
    atmosvars_getlocalmeshqtrcvar,       &
    atmosvars_getlocalmeshqtrcphytend
  use mod_atmos_dyn_vars, only: &
    atmosdynvars,                      &
    atmosdynauxvars_getlocalmeshfields
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
  type, extends(modelcomponentproc), public :: atmosdyn
    type(atmdyndgmdriver_nonhydro3d) :: dyncore_driver
    type(atmdyndgmdriver_trcadv3d) :: trcadv_driver
 
    type(atmosdynvars) :: dyn_vars
 
    ! explicit numerical diffusion
    logical :: calc_numdiff_flag
    type(atmdyn_nonhydro3d_numdiff) :: numdiff
 
  contains
    procedure, public :: setup => atmosdyn_setup 
    procedure, public :: calc_tendency => atmosdyn_calc_tendency
    procedure, public :: update => atmosdyn_update
    procedure, public :: finalize => atmosdyn_finalize
  end type atmosdyn
 
  !-----------------------------------------------------------------------------
  !++ Public parameters & variables
  !-----------------------------------------------------------------------------
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: setup_coriolis_parameter
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
 
contains
 
!OCL SERIAL
  subroutine atmosdyn_setup( this, model_mesh, tm_parent_comp )
    use mod_atmos_mesh, only: atmosmesh
    use scale_time_manager, only: time_manager_component
 
    implicit none
 
    class(atmosdyn), intent(inout) :: this
    class(modelmeshbase), target, intent(in) :: model_mesh
    class(time_manager_component), intent(inout) :: tm_parent_comp
 
    character(len=H_MID) :: EQS_TYPE             = "NONHYDRO3D_HEVE"
    character(len=H_SHORT) :: TINTEG_TYPE        = 'ERK_SSP_3s3o'
    character(len=H_SHORT) :: TINTEG_TYPE_TRACER = 'ERK_SSP_3s3o'
    real(DP) :: TIME_DT                          = undef8            
    character(len=H_SHORT) :: TIME_DT_UNIT       = 'SEC'
    
    logical :: MODALFILTER_FLAG           = .false. 
    logical :: NUMDIFF_FLAG               = .false. 
    logical :: SPONGELAYER_FLAG           = .false. 
    logical :: ONLY_TRACERADV_FLAG        = .false. 
    logical :: TRACERADV_DISABLE_LIMITER  = .false. 
    logical :: TRACERADV_MODALFILTER_FLAG = .false. 
    logical :: HIDE_MPI_COMM_FLAG         = .false.
    
    namelist / param_atmos_dyn /       &
      eqs_type,                        &
      tinteg_type,                     &
      tinteg_type_tracer,              &      
      time_dt,                         &
      time_dt_unit,                    &
      modalfilter_flag,                &
      numdiff_flag,                    &
      spongelayer_flag,                &
      only_traceradv_flag,             &
      traceradv_disable_limiter,       &
      traceradv_modalfilter_flag,      &
      hide_mpi_comm_flag
    
    class(atmosmesh), pointer     :: atm_mesh
    class(meshbase), pointer      :: ptr_mesh
    class(localmeshbase), pointer :: ptr_lcmesh
    class(elementbase3d), pointer :: elem3D
    integer :: n
    real(DP) :: dtsec
 
    class(meshbase3d), pointer :: mesh3D
 
    integer :: ierr
    !--------------------------------------------------
 
    if (.not. this%IsActivated()) return
    log_info('AtmosDyn_setup',*)
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_dyn,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("ATMOS_DYN_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
      log_error("ATMOS_DYN_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_DYN. Check!'
      call prc_abort
    endif
    log_nml(param_atmos_dyn)
    
    !- get mesh --------------------------------------------------
 
    call model_mesh%GetModelMesh( ptr_mesh )
    select type(model_mesh)
    class is (atmosmesh)
      atm_mesh => model_mesh
    end select
    mesh3d => atm_mesh%ptr_mesh
 
    !- Setup the temporal integrator
 
    call tm_parent_comp%Regist_process( 'ATMOS_DYN', time_dt, time_dt_unit, & ! (in)
      this%tm_process_id )                                                    ! (out)
 
    dtsec = tm_parent_comp%process_list(this%tm_process_id)%dtsec
 
    !- initialize the variables 
    call this%dyn_vars%Init( model_mesh )
 
    call setup_coriolis_parameter( this%dyn_vars, atm_mesh )
 
    !- Initialize a module for 3D dynamical core 
    call this%dyncore_driver%Init( eqs_type, &
      tinteg_type, dtsec,                    &
      spongelayer_flag, modalfilter_flag,    &
      hide_mpi_comm_flag, &
      atm_mesh )
 
    !- Initialize a module for tracer equations
    call this%trcadv_driver%Init( "TRCADV3D_HEVE", &
      tinteg_type_tracer, dtsec,                                        &
      traceradv_modalfilter_flag, traceradv_disable_limiter,            &
      atm_mesh, this%dyncore_driver%boundary_cond, only_traceradv_flag  )
 
    !- Setup the numerical diffusion
    this%CALC_NUMDIFF_FLAG = numdiff_flag
    if (this%CALC_NUMDIFF_FLAG) call this%numdiff%Init( atm_mesh, dtsec )
 
    return
  end subroutine atmosdyn_setup
 
 
!OCL SERIAL  
  subroutine atmosdyn_calc_tendency( this, model_mesh, prgvars_list, trcvars_list, auxvars_list, forcing_list, is_update )
    implicit none
    
    class(atmosdyn), 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
    !--------------------------------------------------
    if (.not. this%IsActivated()) return
    !LOG_INFO('AtmosDyn_tendency',*)
 
    return  
  end subroutine atmosdyn_calc_tendency
 
 
!OCL SERIAL
  subroutine atmosdyn_update( this, model_mesh, prgvars_list, trcvars_list, auxvars_list, forcing_list, is_update )
    use scale_tracer, only: &
      qa
    use scale_const, only: &
      grav => const_grav,  &
      rdry => const_rdry,  &
      cpdry => const_cpdry, &
      cvdry => const_cvdry, &
      pres00 => const_pre00
 
    use scale_meshfield_base, only: &
      meshfield2d, meshfield3d
    use scale_atm_dyn_dgm_nonhydro3d_common, only: &
      prgvar_num, &
      prgvar_ddens_id
    use scale_atm_dyn_dgm_driver_trcadv3d, only: &
      trcddens_id => trcvars3d_dens_id, trcddens0_id => trcvars3d_dens0_id, &
      massflx_z_tavg => massflx_z_id, massflx_x_tavg => massflx_x_id, massflx_y_tavg => massflx_y_id
 
    use mod_atmos_dyn_vars, only: &
      atmos_dyn_auxvars2d_coriolis_id
      
    implicit none
 
    class(atmosdyn), 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(meshbase3d), pointer :: mesh3D
    !--------------------------------------------------
    
    call prof_rapstart( 'ATM_DYN_update', 1)   
 
    call model_mesh%GetModelMesh( mesh )
    select type(mesh)
    class is (meshbase3d)
      mesh3d => mesh
    end select
 
    if ( .not. this%trcadv_driver%ONLY_TRACERADV_FLAG ) then
 
      call prof_rapstart( 'ATM_DYN_core', 2)
 
      call this%dyncore_driver%Update( &
        prgvars_list, auxvars_list, forcing_list,                                  & ! (inout)
        this%trcadv_driver%TRCVARS3D(trcddens_id),                                 & ! (inout)
        this%trcadv_driver%TRCVARS3D(trcddens0_id),                                & ! (inout)
        this%trcadv_driver%AUXTRC_FLUX_VARS3D(massflx_x_tavg),                     & ! (inout)
        this%trcadv_driver%AUXTRC_FLUX_VARS3D(massflx_y_tavg),                     & ! (inout)
        this%trcadv_driver%AUXTRC_FLUX_VARS3D(massflx_z_tavg),                     & ! (inout)
        this%trcadv_driver%alphaDensM, this%trcadv_driver%alphaDensP,              & ! (inout)
        this%dyn_vars%AUX_VARS2D(atmos_dyn_auxvars2d_coriolis_id),                 & ! (in)
        model_mesh%element3D_operation,                                            & ! (in)
        model_mesh%DOptrMat(1), model_mesh%DOptrMat(2), model_mesh%DOptrMat(3),    & ! (in)
        model_mesh%SOptrMat(1), model_mesh%SOptrMat(2), model_mesh%SOptrMat(3),    & ! (in)
        model_mesh%LiftOptrMat, mesh3d                                             ) ! (in)
 
      call prof_rapend( 'ATM_DYN_core', 2)
 
    end if
 
    !-- Tracer advection (prepair) ------------------------------------------------
 
    if ( qa > 0 ) then
      call prof_rapstart( 'ATM_DYN_qtracer', 2)
 
      call this%trcadv_driver%Update( &
        trcvars_list, prgvars_list, auxvars_list, forcing_list,                 & ! (inout)
        model_mesh%element3D_operation,                                         & ! (in)
        model_mesh%DOptrMat(1), model_mesh%DOptrMat(2), model_mesh%DOptrMat(3), & ! (in)
        model_mesh%SOptrMat(1), model_mesh%SOptrMat(2), model_mesh%SOptrMat(3), & ! (in)
        model_mesh%LiftOptrMat, mesh3d,                                         & ! (in)
        this%dyncore_driver                                                     ) ! (in)
      
      call prof_rapend( 'ATM_DYN_qtracer', 2)     
    end if
 
    !-- numerical diffusion for dynamical variables -----------------------------
 
    if ( this%CALC_NUMDIFF_FLAG ) then      
      call prof_rapstart( 'ATM_DYN_numfilter', 2)
      call this%numdiff%Apply( prgvars_list, &
        auxvars_list, this%dyncore_driver%boundary_cond, &
        model_mesh%DOptrMat(1), model_mesh%DOptrMat(2), model_mesh%DOptrMat(3),    & ! (in)
        model_mesh%LiftOptrMat, mesh3d                                             ) ! (in)
      call prof_rapend( 'ATM_DYN_numfilter', 2)
    end if
 
    !---------------------------
    call prof_rapend( 'ATM_DYN_update', 1)
 
    return  
  end subroutine atmosdyn_update
 
!OCL SERIAL
  subroutine atmosdyn_finalize( this )
    implicit none
    class(atmosdyn), intent(inout) :: this
 
    integer :: n
    !--------------------------------------------------
    if (.not. this%IsActivated()) return
    log_info('AtmosDyn_finalize',*)
 
    call this%dyncore_driver%Final()
    call this%trcadv_driver%Final()
 
    if (this%CALC_NUMDIFF_FLAG) call this%numdiff%Final()
 
    call this%dyn_vars%Final()
 
    return  
  end subroutine atmosdyn_finalize  
 
  !--- private ---------------
 
!OCL SERIAL
  subroutine setup_coriolis_parameter( this, atm_mesh )
 
    use scale_coriolis_param, only: get_coriolis_parameter
    use scale_mesh_cubedom3d, only: meshcubedom3d
    implicit none
 
    class(atmosdynvars), target, intent(inout) :: this
    class(atmosmesh), target, intent(in) :: atm_mesh
 
    class(localmeshfieldbase), pointer :: coriolis
    class(localmesh3d), pointer :: lcmesh3D
    class(localmesh2d), pointer :: lcmesh2D
    integer :: n
 
    character(len=H_SHORT) :: CORIOLIS_type
    real(RP) :: CORIOLIS_f0   = 0.0_rp
    real(RP) :: CORIOLIS_beta = 0.0_rp
    real(RP) :: CORIOLIS_y0
 
    namelist /param_atmos_dyn_coriolis/ &
      coriolis_type,                         &                
      coriolis_f0, coriolis_beta, coriolis_y0
        
    class(meshbase3d), pointer :: mesh3D
    class(meshcubedom3d), pointer :: meshCube
    integer :: ierr
    !---------------------------------------------------------------
 
    mesh3d => atm_mesh%ptr_mesh
 
    coriolis_type = 'NONE'
 
    select type(mesh3d)
    type is (meshcubedom3d)
      coriolis_y0 = 0.5_rp*(mesh3d%ymax_gl +  mesh3d%ymin_gl)
    end select
 
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_dyn_coriolis,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("ATMOS_DYN_setup_coriolis",*) 'Not found namelist. Default used.'
    else if( ierr > 0 ) then !--- fatal error
      log_error("ATMOS_DYN_setup_coriolis",*) 'Not appropriate names in namelist PARAM_ATMOS_DYN_CORIOLIS. Check!'
      call prc_abort
    end if
    log_nml(param_atmos_dyn_coriolis)
 
    do n = 1, mesh3d%LOCAL_MESH_NUM
      call atmosdynauxvars_getlocalmeshfields( n, mesh3d, this%AUXVARS2D_manager, &
        coriolis, lcmesh3d )
      lcmesh2d => lcmesh3d%lcmesh2D
 
      call get_coriolis_parameter( &
        coriolis%val(:,lcmesh2d%NeS:lcmesh2d%NeE),                       & ! (out)
        coriolis_type, lcmesh2d%refElem2D%Np * lcmesh2d%Ne,              & ! (in)
        lcmesh2d%pos_en(:,:,2), coriolis_f0, coriolis_beta, coriolis_y0, & ! (in)
        lcmesh3d%lat2D(:,:)                                              ) ! (in)
    end do
 
    return
  end subroutine setup_coriolis_parameter
 
!--------
 
!   subroutine cal_MOMZ_tend( &
!     MOMZ_t, MOMZ_t_advx, MOMZ_t_advY, MOMZ_t_advZ, MOMZ_t_lift, MOMZ_t_buoy,     & ! (out)
!      DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd,                    & ! (in)
!      Dx, Dy, Dz, Sx, Sy, Sz, Lift, lmesh, elem, lmesh2D, elem2D )
 
!      use scale_element_base
!      use scale_sparsemat
!      use scale_const, only: &
!       GRAV => CONST_GRAV,  &
!       Rdry => CONST_Rdry,  &
!       CPdry => CONST_CPdry, &
!       CVdry => CONST_CVdry, &
!       PRES00 => CONST_PRE00
!      use scale_atm_dyn_nonhydro3d, only: IntrpMat_VPOrdM1
!      implicit none
 
!      class(LocalMesh3D), intent(in) :: lmesh
!      class(elementbase3D), intent(in) :: elem
!      class(LocalMesh2D), intent(in) :: lmesh2D
!      class(elementbase2D), intent(in) :: elem2D
!      type(SparseMat), intent(in) :: Dx, Dy, Dz, Sx, Sy, Sz, Lift
!      real(RP), intent(out) :: MOMZ_t(elem%Np,lmesh%NeA)
!      real(RP), intent(out) :: MOMZ_t_advx(elem%Np,lmesh%NeA)
!      real(RP), intent(out) :: MOMZ_t_advy(elem%Np,lmesh%NeA)
!      real(RP), intent(out) :: MOMZ_t_advz(elem%Np,lmesh%NeA)
!      real(RP), intent(out) :: MOMZ_t_lift(elem%Np,lmesh%NeA)
!      real(RP), intent(out) :: MOMZ_t_buoy(elem%Np,lmesh%NeA)
 
!      real(RP), intent(in)  :: DDENS_(elem%Np,lmesh%NeA)
!      real(RP), intent(in)  :: MOMX_(elem%Np,lmesh%NeA)
!      real(RP), intent(in)  :: MOMY_(elem%Np,lmesh%NeA)
!      real(RP), intent(in)  :: MOMZ_(elem%Np,lmesh%NeA)
!      real(RP), intent(in)  :: DRHOT_(elem%Np,lmesh%NeA)
!      real(RP), intent(in)  :: DENS_hyd(elem%Np,lmesh%NeA)
!      real(RP), intent(in)  :: PRES_hyd(elem%Np,lmesh%NeA)
 
!      real(RP) :: Fx(elem%Np), Fy(elem%Np), Fz(elem%Np), LiftDelFlx(elem%Np)
!      real(RP) :: del_flux(elem%NfpTot,lmesh%Ne)
!      real(RP) :: dens_(elem%Np), RHOT_(elem%Np), dpres_(elem%Np)
!      real(RP) :: pres_(elem%Np), u_(elem%Np), v_(elem%Np), w_(elem%Np)
 
!      integer :: ke
!      !------------------------------------------------------------------------
 
!      call cal_del_flux_dyn( del_flux,                                          & ! (out)
!        DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd,                & ! (in)
!        lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
!        lmesh%vmapM, lmesh%vmapP,                                               & ! (in)
!        lmesh, elem )                                                             ! (in)
  
!      !-----
!      !$omp parallel do private(RHOT_,pres_,dpres_,dens_,u_,v_,w_,Fx,Fy,Fz,LiftDelFlx)
!      do ke = lmesh%NeS, lmesh%NeE
!        !--
 
!        RHOT_(:) = PRES00/Rdry * (PRES_hyd(:,ke)/PRES00)**(CVdry/CPdry) + DRHOT_(:,ke)
!        pres_(:) = PRES00 * (Rdry*RHOT_(:)/PRES00)**(CPdry/Cvdry)
!        dpres_(:) = pres_(:) - PRES_hyd(:,ke)
!        dens_(:) = DDENS_(:,ke) + DENS_hyd(:,ke)
 
!        u_(:) = MOMX_(:,ke)/dens_(:)
!        v_(:) = MOMY_(:,ke)/dens_(:)
!        w_(:) = MOMZ_(:,ke)/dens_(:)
 
!        !-- MOMZ
!        call sparsemat_matmul(Dx, u_(:)*MOMZ_(:,ke), Fx)
!        call sparsemat_matmul(Dy, v_(:)*MOMZ_(:,ke), Fy)
!        call sparsemat_matmul(Dz, w_(:)*MOMZ_(:,ke), Fz)
!        MOMZ_t_advx(:,ke) = - lmesh%Escale(:,ke,1,1) * Fx(:)
!        MOMZ_t_advy(:,ke) = - lmesh%Escale(:,ke,2,2) * Fy(:)
!        MOMZ_t_advz(:,ke) = - lmesh%Escale(:,ke,3,3) * Fz(:)
 
!        call sparsemat_matmul(Dz, dpres_(:), Fz)
!        MOMZ_t_buoy(:,ke) = - lmesh%Escale(:,ke,3,3) * Fz(:) &
!                            - matmul(IntrpMat_VPOrdM1, DDENS_(:,ke)) * Grav
 
!        call sparsemat_matmul(Lift, lmesh%Fscale(:,ke)*del_flux(:,ke), LiftDelFlx)
!        MOMZ_t_lift(:,ke) = - LiftDelFlx(:)
 
!        MOMZ_t(:,ke) = MOMZ_t_advx(:,ke) +  MOMZ_t_advy(:,ke) +  MOMZ_t_advz(:,ke) &
!                     + MOMZ_t_lift(:,ke) + MOMZ_t_buoy(:,ke)
!      end do
 
!      return
!  end subroutine cal_MOMZ_tend
 
 
!  subroutine cal_del_flux_dyn( del_flux, &
!    DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd,   &
!    nx, ny, nz, vmapM, vmapP, lmesh, elem )
 
!    use scale_const, only: &
!     GRAV => CONST_GRAV,  &
!     Rdry => CONST_Rdry,  &
!     CPdry => CONST_CPdry, &
!     CVdry => CONST_CVdry, &
!     PRES00 => CONST_PRE00
  
!    implicit none
 
!    class(LocalMesh3D), intent(in) :: lmesh
!    class(elementbase3D), intent(in) :: elem  
!    real(RP), intent(out) ::  del_flux(elem%NfpTot*lmesh%Ne)
!    real(RP), intent(in) ::  DDENS_(elem%Np*lmesh%NeA)
!    real(RP), intent(in) ::  MOMX_(elem%Np*lmesh%NeA)  
!    real(RP), intent(in) ::  MOMY_(elem%Np*lmesh%NeA)  
!    real(RP), intent(in) ::  MOMZ_(elem%Np*lmesh%NeA)  
!    real(RP), intent(in) ::  DRHOT_(elem%Np*lmesh%NeA)  
!    real(RP), intent(in) ::  DENS_hyd(elem%Np*lmesh%NeA)
!    real(RP), intent(in) ::  PRES_hyd(elem%Np*lmesh%NeA)
!    real(RP), intent(in) :: nx(elem%NfpTot*lmesh%Ne)
!    real(RP), intent(in) :: ny(elem%NfpTot*lmesh%Ne)
!    real(RP), intent(in) :: nz(elem%NfpTot*lmesh%Ne)
!    integer, intent(in) :: vmapM(elem%NfpTot*lmesh%Ne)
!    integer, intent(in) :: vmapP(elem%NfpTot*lmesh%Ne)
   
!    integer :: i, iP, iM
!    real(RP) :: VelP, VelM, alpha
!    real(RP) :: uM, uP, vM, vP, wM, wP, presM, presP, dpresM, dpresP, densM, densP, rhotM, rhotP, rhot_hyd_M, rhot_hyd_P
!    real(RP) :: gamm, rgamm
!    !------------------------------------------------------------------------
 
!    gamm = CpDry/CvDry
!    rgamm = CvDry/CpDry
 
!    !$omp parallel do private( &
!    !$omp iM, iP, uM, VelP, VelM, alpha, &
!    !$omp uP, vM, vP, wM, wP, presM, presP, dpresM, dpresP, densM, densP, rhotM, rhotP, rhot_hyd_M, rhot_hyd_P)
!    do i=1, elem%NfpTot*lmesh%Ne
!      iM = vmapM(i); iP = vmapP(i)
 
!      rhot_hyd_M = PRES00/Rdry * (PRES_hyd(iM)/PRES00)**rgamm
!      rhot_hyd_P = PRES00/Rdry * (PRES_hyd(iP)/PRES00)**rgamm
     
!      rhotM = rhot_hyd_M + DRHOT_(iM)
!      presM = PRES00 * (Rdry*rhotM/PRES00)**gamm
!      dpresM = presM - PRES_hyd(iM)*abs(nz(i))
 
!      rhotP = rhot_hyd_P + DRHOT_(iP) 
!      presP = PRES00 * (Rdry*rhotP/PRES00)**gamm
!      dpresP = presP - PRES_hyd(iP)*abs(nz(i))
 
!      densM = DDENS_(iM) + DENS_hyd(iM)
!      densP = DDENS_(iP) + DENS_hyd(iP)
 
!      VelM = (MOMX_(iM)*nx(i) + MOMY_(iM)*ny(i) + MOMZ_(iM)*nz(i))/densM
!      VelP = (MOMX_(iP)*nx(i) + MOMY_(iP)*ny(i) + MOMZ_(iP)*nz(i))/densP
 
!      alpha = max( sqrt(gamm*presM/densM) + abs(VelM), sqrt(gamm*presP/densP) + abs(VelP)  )
 
 
!      del_flux(i) = 0.5_RP*(                &
!                    ( MOMZ_(iP)*VelP - MOMZ_(iM)*VelM)   &
!                    + ( dpresP - dpresM )*nz(i)          &                    
!                    - alpha*(MOMZ_(iP) - MOMZ_(iM))      )
!    end do
 
!    return
!  end subroutine cal_del_flux_dyn
 
end module mod_atmos_dyn