scale_atm_dyn_dgm_spongelayer.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_spongelayer
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
 
  use scale_element_base, only: &
    elementbase, elementbase3d
  use scale_localmesh_base, only: localmeshbase
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_base3d, only: meshbase3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
  type, public :: atmdynspongelayer
    real(rp) :: wdamp_tau
    real(rp) :: wdamp_height
    logical  :: hveldamp_flag
  contains
    procedure :: init => atm_dyn_dgm_spongelayer_init
    procedure :: final => atm_dyn_dgm_spongelayer_final
    procedure :: addtend => atm_dyn_dgm_spongelayer_add_tend
  end type atmdynspongelayer
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  private :: calc_wdampcoef
 
contains
 
!OCL SERIAL
  subroutine atm_dyn_dgm_spongelayer_init( this, mesh3D, dtsec )
    use scale_prc, only: prc_abort
    use scale_mesh_cubedom3d, only: meshcubedom3d
    use scale_mesh_cubedspheredom3d, only: meshcubedspheredom3d
    implicit none
 
    class(atmdynspongelayer), target, intent(inout) :: this
    class(meshbase3d), target, intent(in) :: mesh3D
    real(RP), intent(in) :: dtsec
 
    real(RP) :: SL_WDAMP_TAU        = -1.0_rp ! the maximum tau for Rayleigh damping of w [s]
    real(RP) :: SL_WDAMP_HEIGHT     = -1.0_rp ! the height to start apply Rayleigh damping [m]
    integer  :: SL_WDAMP_LAYER      = -1      ! the vertical number of finite element to start apply Rayleigh damping [num]
    logical  :: SL_HORIVELDAMP_FLAG = .false. ! Is the horizontal velocity damped? 
    
    namelist /param_atmos_dyn_spongelayer/ &
      sl_wdamp_tau,                        &                
      sl_wdamp_height,                     &
      sl_wdamp_layer,                      &
      sl_horiveldamp_flag
    
    class(localmesh3d), pointer :: lcmesh3D
    class(elementbase3d), pointer :: elem3D
  
    integer :: NeGZ
    integer :: ierr
    !---------------------------------------------------------------
 
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_dyn_spongelayer,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("ATMOS_DYN_setup_spongelayer",*) 'Not found namelist. Default used.'
    else if( ierr > 0 ) then !--- fatal error
      log_error("ATMOS_DYN_setup_spongelayer",*) 'Not appropriate names in namelist PARAM_ATMOS_DYN_SPONGELAYER. Check!'
      call prc_abort
    end if
    log_nml(param_atmos_dyn_spongelayer)
 
    this%wdamp_tau    = sl_wdamp_tau 
    this%wdamp_height = sl_wdamp_height
 
    lcmesh3d => mesh3d%lcmesh_list(1)
    elem3d => lcmesh3d%refElem3D
    select type(mesh3d)
    type is (meshcubedom3d)
      negz = mesh3d%NeGZ
    type is (meshcubedspheredom3d)
      negz = mesh3d%NeGZ
    end select
 
    if ( sl_wdamp_layer > negz ) then
      log_error("ATMOS_DYN_setup_spongelayer",*) 'SL_wdamp_layer should be less than total of vertical elements (NeGZ). Check!'
      call prc_abort
    else if( sl_wdamp_layer > 0 ) then
      this%wdamp_height = lcmesh3d%pos_en(1,1+(sl_wdamp_layer-1)*lcmesh3d%NeX*lcmesh3d%NeY,3)
    end if
    if ( this%wdamp_tau < 0.0_rp ) then
      this%wdamp_tau = dtsec * 10.0_rp
    else if ( this%wdamp_tau < dtsec ) then
      log_error("ATMOS_DYN_setup_spongelayer",*) 'SL_wdamp_tau should be larger than TIME_DT (ATMOS_DYN). Check!'
      call prc_abort
    end if
    
    this%hveldamp_flag = sl_horiveldamp_flag
 
    return
  end subroutine atm_dyn_dgm_spongelayer_init
 
!OCL SERIAL
  subroutine atm_dyn_dgm_spongelayer_final( this )
    implicit none
    class(atmdynspongelayer), target, intent(inout) :: this
    !---------------------------------------------------------------
 
    return
  end subroutine atm_dyn_dgm_spongelayer_final
 
  subroutine atm_dyn_dgm_spongelayer_add_tend( this, MOMX_dt, MOMY_dt, MOMZ_dt, &
    MOMX_, MOMY_, MOMZ_,                                                        &
    lmesh, elem   )
 
    implicit none
 
    class(atmdynspongelayer), intent(in) :: this
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    real(RP), intent(inout) :: MOMX_dt(elem%Np,lmesh%NeA)
    real(RP), intent(inout) :: MOMY_dt(elem%Np,lmesh%NeA)
    real(RP), intent(inout) :: MOMZ_dt(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)
 
    integer :: ke
    integer :: ke_x, ke_y, ke_z
    integer :: keZtop
    real(RP) :: wdamp_coef(elem%Np)
    real(RP) :: zTop(elem%Nnode_h1D**2)
    real(RP) :: s
    !-----------------------------------------------------------------
 
    if ( this%hveldamp_flag ) then
      s = 1.0_rp
    else
      s = 0.0_rp
    end if
 
    !$omp parallel do collapse(3) private(ke,keZtop,zTop,wdamp_coef)
    do ke_z = 1, lmesh%NeZ
    do ke_y = 1, lmesh%NeY
    do ke_x = 1, lmesh%NeX
      ke = ke_x + (ke_y-1)*lmesh%NeX + (ke_z-1)*lmesh%NeX*lmesh%NeY
      keztop =  ke_x + (ke_y-1)*lmesh%NeX + (lmesh%NeZ-1)*lmesh%NeX*lmesh%NeY
      ztop(:) = lmesh%pos_en(elem%Hslice(:,elem%Nnode_v),keztop,3)
 
      call calc_wdampcoef( &
        this%wdamp_tau, this%wdamp_height, lmesh%pos_en(:,ke,3), ztop(:), &
        elem%Nnode_h1D, elem%Nnode_v,                                     &
        wdamp_coef(:) )
 
      momx_dt(:,ke) = momx_dt(:,ke) - s * wdamp_coef(:) * momx_(:,ke)
      momy_dt(:,ke) = momy_dt(:,ke) - s * wdamp_coef(:) * momy_(:,ke)
      momz_dt(:,ke) = momz_dt(:,ke) - wdamp_coef(:) * momz_(:,ke)
    end do
    end do
    end do
 
    return
  end subroutine atm_dyn_dgm_spongelayer_add_tend
 
!-- private ------------------------------
 
!OCL SERIAL
  subroutine calc_wdampcoef( &
    wdamp_tau, wdamp_height, z, zTop, Nnode_h1D, Nnode_v, & ! (in)
    wdamp_coef                                            ) ! (out)
 
    use scale_const, only: &
      pi => const_pi
    implicit none
 
    integer, intent(in) :: Nnode_h1D
    integer, intent(in) :: Nnode_v
    real(RP), intent(out) :: wdamp_coef(Nnode_h1D**2,Nnode_v)
    real(RP), intent(in) :: wdamp_tau
    real(RP), intent(in) :: wdamp_height
    real(RP), intent(in) :: z(Nnode_h1D**2,Nnode_v)
    real(RP), intent(in) :: zTop(Nnode_h1D**2)
 
    integer :: p_z
    real(RP) :: sw(Nnode_h1D**2)
    real(RP) :: r_wdamp_tau
    !-----------------------------------------------------------------
 
    r_wdamp_tau = 1.0_rp / wdamp_tau
    do p_z=1, nnode_v
      wdamp_coef(:,p_z) = 0.25_rp * r_wdamp_tau                                        &
        * ( 1.0_rp + sign( 1.0_rp, z(:,p_z) - wdamp_height )                         ) &
        * ( 1.0_rp - cos( pi * (z(:,p_z) - wdamp_height)/(ztop(:) - wdamp_height) )  )
    end do
 
    return
  end subroutine calc_wdampcoef
 
end module scale_atm_dyn_dgm_spongelayer