scale_atm_dyn_dgm_bnd.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_bnd
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
  use scale_const, only: &
    grav => const_grav,  &
    rdry => const_rdry,  &
    cpdry => const_cpdry, &
    cvdry => const_cvdry, &
    pres00 => const_pre00
 
  use scale_element_base, only: &
    elementbase, elementbase2d, elementbase3d
  use scale_mesh_base, only: meshbase
  use scale_localmesh_base, only: localmeshbase
 
  use scale_mesh_base3d, only: meshbase3d
  use scale_mesh_base2d, only: meshbase2d  
  use scale_localmesh_3d, only: localmesh3d
  use scale_localmesh_2d, only: localmesh2d  
 
  use scale_localmeshfield_base, only: localmeshfield3d
  use scale_meshfield_base, only: meshfield3d
 
  use scale_mesh_bndinfo, only: meshbndinfo
 
  use scale_atm_dyn_dgm_nonhydro3d_common, only: &
    dens_vid => prgvar_ddens_id, rhot_vid => prgvar_drhot_id, &
    momx_vid => prgvar_momx_id, momy_vid => prgvar_momy_id,   &
    momz_vid => prgvar_momz_id,                               &
    prgvar_num
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public parameter & type & procedures
  !
  type, public :: atmdynbnd
    type(meshbndinfo), allocatable :: velbc_list(:)
    type(meshbndinfo), allocatable :: thermalbc_list(:)
 
    integer, allocatable:: velbc_ids(:)
    integer, allocatable :: thermalbc_ids(:)
    real(rp), allocatable :: thermal_fixval(:)
  contains
    procedure :: init => atmos_dyn_bnd_setup
    procedure :: final => atmos_dyn_bnd_finalize
    procedure :: setbcinfo => atmos_dyn_bnd_setbcinfo
    procedure :: applybc_progvars_lc => atmos_dyn_bnd_applybc_prgvars_lc
    procedure :: applybc_numdiff_odd_lc => atmos_dyn_bnd_applybc_numdiff_odd_lc
    procedure :: applybc_numdiff_even_lc => atmos_dyn_bnd_applybc_numdiff_even_lc
    procedure :: applybc_grad_tbvars_lc => atmos_dyn_bnd_applybc_tbvars_lc
    procedure :: applybc_grad_tbstress_lc => atmos_dyn_bnd_applybc_tbstress_lc
    procedure :: inquire_bound_flag => atmos_dyn_bnd_inquire_bound_flag
  end type
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures
  !
  !-------------------
 
 
  integer, parameter :: dombnd_south_id = 1
  integer, parameter :: dombnd_east_id  = 2
  integer, parameter :: dombnd_north_id = 3
  integer, parameter :: dombnd_west_id  = 4
  integer, parameter :: dombnd_btm_id   = 5
  integer, parameter :: dombnd_top_id   = 6
  integer, parameter :: dom_bnd_num     = 6
 
contains
 
!OCL SERIAL
  subroutine atmos_dyn_bnd_setup( this )
    use scale_const, only: &
      undef8 => const_undef8    
    use scale_mesh_bndinfo, only: &
      bnd_type_nospec_name, &
      bnd_type_periodic_id, &
      bndtype_nametoid
    implicit none
    
    class(atmdynbnd), intent(inout) :: this
 
    character(len=H_SHORT) :: btm_vel_bc, top_vel_bc
    character(len=H_SHORT) :: north_vel_bc, south_vel_bc
    character(len=H_SHORT) :: east_vel_bc, west_vel_bc
    character(len=H_SHORT) :: btm_thermal_bc, top_thermal_bc
    character(len=H_SHORT) :: north_thermal_bc, south_thermal_bc
    character(len=H_SHORT) :: east_thermal_bc, west_thermal_bc
 
    real(rp) :: btm_thermal_fixval, top_thermal_fixval
    real(rp) :: north_thermal_fixval, south_thermal_fixval
    real(rp) :: east_thermal_fixval, west_thermal_fixval
 
    namelist /param_atmos_dyn_bnd/ &
      btm_vel_bc, top_vel_bc, north_vel_bc, south_vel_bc, east_vel_bc, west_vel_bc,                         &
      btm_thermal_bc, top_thermal_bc, north_thermal_bc, south_thermal_bc, east_thermal_bc, west_thermal_bc, &
      btm_thermal_fixval, top_thermal_fixval, north_thermal_fixval, south_thermal_fixval, east_thermal_fixval, west_thermal_fixval
    
    integer :: ierr
    !-----------------------------------------------
 
    btm_vel_bc   = bnd_type_nospec_name
    top_vel_bc   = bnd_type_nospec_name
    north_vel_bc  = bnd_type_nospec_name
    south_vel_bc = bnd_type_nospec_name
    east_vel_bc  = bnd_type_nospec_name
    west_vel_bc = bnd_type_nospec_name
 
    btm_thermal_bc   = bnd_type_nospec_name
    top_thermal_bc   = bnd_type_nospec_name
    north_thermal_bc  = bnd_type_nospec_name
    south_thermal_bc = bnd_type_nospec_name
    east_thermal_bc  = bnd_type_nospec_name
    west_thermal_bc = bnd_type_nospec_name
 
    btm_thermal_fixval = undef8
    top_thermal_fixval = undef8    
    north_thermal_fixval = undef8
    south_thermal_fixval = undef8    
    east_thermal_fixval = undef8
    west_thermal_fixval = undef8    
 
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_dyn_bnd,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_dyn_bnd_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_dyn_bnd_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_DYN_BND. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_dyn_bnd)
    
    !--
    allocate( this%velBC_ids(dom_bnd_num) )
    this%velBC_ids(dombnd_btm_id) = bndtype_nametoid(btm_vel_bc)
    this%velBC_ids(dombnd_top_id) = bndtype_nametoid(top_vel_bc)
    this%velBC_ids(dombnd_north_id) = bndtype_nametoid(north_vel_bc)
    this%velBC_ids(dombnd_south_id) = bndtype_nametoid(south_vel_bc)
    this%velBC_ids(dombnd_east_id) = bndtype_nametoid(east_vel_bc)
    this%velBC_ids(dombnd_west_id) = bndtype_nametoid(west_vel_bc)
 
    allocate( this%thermalBC_ids(dom_bnd_num) )
    this%thermalBC_ids(dombnd_btm_id) = bndtype_nametoid(btm_thermal_bc)
    this%thermalBC_ids(dombnd_top_id) = bndtype_nametoid(top_thermal_bc)
    this%thermalBC_ids(dombnd_north_id) = bndtype_nametoid(north_thermal_bc)
    this%thermalBC_ids(dombnd_south_id) = bndtype_nametoid(south_thermal_bc)    
    this%thermalBC_ids(dombnd_east_id) = bndtype_nametoid(east_thermal_bc)
    this%thermalBC_ids(dombnd_west_id) = bndtype_nametoid(west_thermal_bc)   
    
    allocate( this%thermal_fixval(dom_bnd_num) )
    this%thermal_fixval(dombnd_btm_id) = btm_thermal_fixval
    this%thermal_fixval(dombnd_top_id) = top_thermal_fixval
    this%thermal_fixval(dombnd_north_id) = north_thermal_fixval
    this%thermal_fixval(dombnd_south_id) = south_thermal_fixval
    this%thermal_fixval(dombnd_east_id) = east_thermal_fixval
    this%thermal_fixval(dombnd_west_id) = west_thermal_fixval
 
    !------
 
    return
  end subroutine atmos_dyn_bnd_setup
 
!OCL SERIAL
  subroutine atmos_dyn_bnd_finalize( this )
    implicit none
    class(atmdynbnd), intent(inout) :: this
 
    integer :: n
    !--------------------------------------
 
    if ( allocated(this%VelBC_list) ) then
      do n=1, size(this%VelBC_list)
        call this%VelBC_list(n)%Final()
        call this%ThermalBC_list(n)%Final()
      end do
      deallocate( this%VelBC_list, this%ThermalBC_list )
 
      deallocate( this%velBC_ids, this%thermalBC_ids )
    end if
 
    return
  end subroutine atmos_dyn_bnd_finalize  
 
!OCL SERIAL
  subroutine atmos_dyn_bnd_setbcinfo( this, mesh )
 
    implicit none
 
    class(atmdynbnd), intent(inout) :: this
    class(meshbase), target, intent(in) :: mesh
    
    integer :: n
    integer :: b
    class(localmeshbase), pointer :: ptr_lcmesh
    class(localmesh3d), pointer :: lcmesh3d
    !--------------------------------------------------
 
 
    allocate( this%VelBC_list(mesh%LOCAL_MESH_NUM) )
    allocate( this%ThermalBC_list(mesh%LOCAL_MESH_NUM) )
 
    nullify( lcmesh3d )
 
    do n=1, mesh%LOCAL_MESH_NUM
      call mesh%GetLocalMesh( n, ptr_lcmesh )
      select type (ptr_lcmesh)
      type is (localmesh3d)
        lcmesh3d => ptr_lcmesh
      end select
 
      call bnd_init_lc( &
        this%VelBC_list(n), this%ThermalBC_list(n),        & ! (inout)
        this%velBC_ids(:), this%thermalBC_ids(:),          & ! (in)
        this%thermal_fixval(:),                            & ! (in)
        lcmesh3d%VMapB, mesh, lcmesh3d, lcmesh3d%refElem3D ) ! (in)
    end do
 
    return
  end subroutine atmos_dyn_bnd_setbcinfo
 
!OCL SERIAL
  subroutine atmos_dyn_bnd_applybc_prgvars_lc( this,  &
    domID,                                              & ! (in)
    ddens, momx, momy, momz, therm,                     & ! (inout)
    dens_hyd, pres_hyd,                                 & ! (in)
    gsqrt, gsqrth, g11, g12, g22, g13, g23, nx, ny, nz, & ! (in)
    vmapm, vmapp, vmapb, lmesh, elem, lmesh2d, elem2d   ) ! (in)
 
    use scale_mesh_bndinfo, only: &
      bnd_type_slip_id, bnd_type_noslip_id, &
      bnd_type_adiabat_id
        
    implicit none
 
    class(atmdynbnd), intent(in) :: this    
    integer, intent(in) :: domid
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    class(localmesh2d), intent(in) :: lmesh2d
    class(elementbase2d), intent(in) :: elem2d
    real(rp), intent(inout) :: ddens(elem%np*lmesh%nea)
    real(rp), intent(inout) :: momx(elem%np*lmesh%nea)
    real(rp), intent(inout) :: momy(elem%np*lmesh%nea)
    real(rp), intent(inout) :: momz(elem%np*lmesh%nea)
    real(rp), intent(inout) :: therm(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) :: gsqrt(elem%np*lmesh%nea)
    real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g11(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g12(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g22(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g13(elem%np*lmesh%nea)
    real(rp), intent(in) ::  g23(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, intent(in) :: vmapb(:)
 
    integer :: p, ke, ke2d
    integer :: i, i_, im, ip
    real(rp) :: mom_normal
 
    real(rp) :: momw
    real(rp) :: gsqrtv, g11_, g12_, g22_
    real(rp) :: fac
    !-----------------------------------------------
 
    !$omp parallel do collapse(2) private( &
    !$omp ke, p, ke2D, i, i_, iM, iP,      &
    !$omp mom_normal, MOMW, GsqrtV, G11_, G12_, G22_, fac    )
    do ke=lmesh%NeS, lmesh%NeE
    do p=1, elem%NfpTot
      i = p + (ke-1)*elem%NfpTot
      ip = vmapp(i)
      i_ = ip - elem%Np*lmesh%NeE
      
      if (i_ > 0) then
        im = vmapm(i)
 
        select case( this%VelBC_list(domid)%list(i_) )
        case ( bnd_type_slip_id)
          ke2d = lmesh%EMap3Dto2D(ke)
 
          gsqrtv = gsqrt(im) / gsqrth(elem%IndexH2Dto3D_bnd(p),ke2d) 
          g11_ = g11(elem%IndexH2Dto3D_bnd(p),ke2d)
          g12_ = g12(elem%IndexH2Dto3D_bnd(p),ke2d)
          g22_ = g22(elem%IndexH2Dto3D_bnd(p),ke2d)
 
          momw = momz(im) / gsqrtv &
             + g13(im) * momx(im) + g23(im) * momy(im)
          fac = nz(i) * gsqrtv**2 / ( 1.0_rp + g11_ * ( gsqrtv * g13(im) )**2 + 2.0_rp * g12_ * ( gsqrtv**2 * g13(im) * g23(im) ) + g22_ * ( gsqrtv * g23(im) )**2 )
 
          mom_normal = momx(im) * nx(i) + momy(im) * ny(i) + momw * nz(i)
          momx(ip) = momx(im) - 2.0_rp * mom_normal * ( nx(i) + fac * ( g11_ * g13(im) + g12_ * g23(im) ) )
          momy(ip) = momy(im) - 2.0_rp * mom_normal * ( ny(i) + fac * ( g12_ * g13(im) + g22_ * g23(im) ) )
          momz(ip) = momz(im) - 2.0_rp * mom_normal * fac / gsqrtv
 
        case ( bnd_type_noslip_id )
          momx(ip) = - momx(im)
          momy(ip) = - momy(im)
          momz(ip) = - momz(im)          
        end select
      end if
 
    end do
    end do
    
    return
  end  subroutine atmos_dyn_bnd_applybc_prgvars_lc
  
!OCL SERIAL
  subroutine atmos_dyn_bnd_applybc_numdiff_odd_lc(  this, & ! (in)
    gxvar, gyvar, gzvar,                                  & ! (inout)
    is_bound,                                             & ! (out)
    varid, domid,                                         & ! (in)
    nx, ny, nz, vmapm, vmapp, vmapb, lmesh, elem )          ! (in)
 
    use scale_mesh_bndinfo, only: &
      bnd_type_slip_id, bnd_type_noslip_id, &
      bnd_type_adiabat_id
    
    implicit none
 
    class(atmdynbnd), intent(in) :: this
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem    
    real(rp), intent(inout) :: gxvar(elem%np*lmesh%nea)
    real(rp), intent(inout) :: gyvar(elem%np*lmesh%nea)
    real(rp), intent(inout) :: gzvar(elem%np*lmesh%nea)
    logical, intent(out) :: is_bound(elem%nfptot*lmesh%ne)
    integer, intent(in) :: varid
    integer, intent(in) :: domid
    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, intent(in) :: vmapb(:)
 
    integer :: i, i_, im, ip
    real(rp) :: grad_normal
    !-----------------------------------------------
 
    do i=1, elem%NfpTot*lmesh%Ne
      ip = vmapp(i)
      i_ = ip - elem%Np*lmesh%NeE
      is_bound(i) = .false.
      
      if (i_ > 0) then  
        im = vmapm(i)
        grad_normal = gxvar(im) * nx(i) + gyvar(im) * ny(i) + gzvar(im) * nz(i)
 
        if ( this%VelBC_list(domid)%list(i_) == bnd_type_slip_id ) then
          select case(varid)
          case(momx_vid)
            gyvar(ip) = gyvar(im) - 2.0_rp * grad_normal * ny(i)
            gzvar(ip) = gzvar(im) - 2.0_rp * grad_normal * nz(i)
          case(momy_vid)          
            gxvar(ip) = gxvar(im) - 2.0_rp * grad_normal * nx(i)
            gzvar(ip) = gzvar(im) - 2.0_rp * grad_normal * nz(i)
          case(momz_vid)          
            gxvar(ip) = gxvar(im) - 2.0_rp * grad_normal * nx(i)
            gyvar(ip) = gyvar(im) - 2.0_rp * grad_normal * ny(i)
          end select
          is_bound(i) = .true.
 
        end if
        if ( this%ThermalBC_list(domid)%list(i_) == bnd_type_adiabat_id ) then
          select case(varid)
          case(dens_vid, rhot_vid)
            gxvar(ip) = gxvar(im) - 2.0_rp * grad_normal * nx(i)
            gyvar(ip) = gyvar(im) - 2.0_rp * grad_normal * ny(i)
            gzvar(ip) = gzvar(im) - 2.0_rp * grad_normal * nz(i)
          end select
          is_bound(i) = .true.
 
        end if
      
      end if
 
    end do
 
    return
  end subroutine atmos_dyn_bnd_applybc_numdiff_odd_lc
 
!OCL SERIAL
  subroutine atmos_dyn_bnd_applybc_numdiff_even_lc(  this, & ! (in)
    var,                                                   & ! (inout)
    is_bound,                                              & ! (out)
    varid, domid,                                          & ! (in)
    nx, ny, nz, vmapm, vmapp, vmapb, lmesh, elem )           ! (in)
 
    use scale_mesh_bndinfo, only: &
      bnd_type_slip_id, bnd_type_noslip_id, &
      bnd_type_adiabat_id
    implicit none
 
    class(atmdynbnd), intent(in) :: this
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem    
    real(rp), intent(inout) :: var(elem%np*lmesh%nea)
    logical, intent(out) :: is_bound(elem%nfptot*lmesh%ne)
    integer, intent(in) :: varid
    integer, intent(in) :: domid
    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, intent(in) :: vmapb(:)
 
    integer :: i, i_, im, ip
    real(rp) :: grad_normal
    !-----------------------------------------------
 
    do i=1, elem%NfpTot*lmesh%Ne
      ip = vmapp(i)
      i_ = ip - elem%Np*lmesh%NeE
      is_bound(i) = .false.
 
      if (i_ > 0) then  
        im = vmapm(i)
 
        if ( this%VelBC_list(domid)%list(i_) == bnd_type_slip_id ) then
          select case(varid)
          case(momx_vid)
            grad_normal = var(im) * nx(i)         
            var(ip) = var(im) - 2.0_rp * grad_normal * nx(i)
          case(momy_vid) 
            grad_normal = var(im) * ny(i)         
            var(ip) = var(im) - 2.0_rp * grad_normal * ny(i)
          case(momz_vid)          
            grad_normal = var(im) * nz(i)         
            var(ip) = var(im) - 2.0_rp * grad_normal * nz(i)
          end select
          is_bound(i) = .true.
 
        else if ( this%VelBC_list(domid)%list(i_) == bnd_type_noslip_id ) then          
          select case(varid)
          case( momx_vid, momy_vid, momz_vid )
            var(ip) = - var(im)
          end select
          is_bound(i) = .true.
        end if
 
      end if
    end do
 
    return
  end subroutine atmos_dyn_bnd_applybc_numdiff_even_lc
 
 
!OCL SERIAL
  subroutine atmos_dyn_bnd_applybc_tbvars_lc( this,  &
    domID,                                              & ! (in)
    ddens, momx, momy, momz, pt, pres,                  & ! (inout)
    dens_hyd, pres_hyd, rtot, cptot,                    & ! (in)
    gsqrt, gsqrth, g11, g12, g22, g13, g23, nx, ny, nz, & ! (in)
    vmapm, vmapp, vmapb, lmesh, elem, lmesh2d, elem2d   ) ! (in)
 
    use scale_const, only: &
      pres00 => const_pre00   
    use scale_mesh_bndinfo, only: &
      bnd_type_slip_id, bnd_type_noslip_id, &
      bnd_type_adiabat_id, bnd_type_fixval_id
        
    implicit none
 
    class(atmdynbnd), intent(in) :: this    
    integer, intent(in) :: domid
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    class(localmesh2d), intent(in) :: lmesh2d
    class(elementbase2d), intent(in) :: elem2d
    real(rp), intent(inout) :: ddens(elem%np*lmesh%nea)
    real(rp), intent(inout) :: momx(elem%np*lmesh%nea)
    real(rp), intent(inout) :: momy(elem%np*lmesh%nea)
    real(rp), intent(inout) :: momz(elem%np*lmesh%nea)
    real(rp), intent(inout) :: pt(elem%np*lmesh%nea)
    real(rp), intent(inout) :: pres(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) :: rtot(elem%np*lmesh%nea)    
    real(rp), intent(in) :: cptot(elem%np*lmesh%nea)    
    real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
    real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g11(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g12(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g22(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g13(elem%np*lmesh%nea)
    real(rp), intent(in) ::  g23(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, intent(in) :: vmapb(:)
 
    integer :: p, ke, ke2d
    integer :: i, i_, im, ip
    real(rp) :: mom_normal
 
    real(rp) :: momw
    real(rp) :: temp_p, temp_b
    real(rp) :: gsqrtv, g11_, g12_, g22_
    real(rp) :: fac
    !-----------------------------------------------
 
    !$omp parallel do collapse(2) private( &
    !$omp ke, p, ke2D, i, i_, iM, iP,      &
    !$omp mom_normal, MOMW, GsqrtV, G11_, G12_, G22_, fac, &
    !$omp TEMP_P, TEMP_B )
    do ke=lmesh%NeS, lmesh%NeE
    do p=1, elem%NfpTot
      i = p + (ke-1)*elem%NfpTot
      ip = vmapp(i)
      i_ = ip - elem%Np*lmesh%NeE
      
      if (i_ > 0) then
        im = vmapm(i)
 
        select case( this%VelBC_list(domid)%list(i_) )
        case ( bnd_type_slip_id)
          ke2d = lmesh%EMap3Dto2D(ke)
 
          gsqrtv = gsqrt(im) / gsqrth(elem%IndexH2Dto3D_bnd(p),ke2d) 
          g11_ = g11(elem%IndexH2Dto3D_bnd(p),ke2d)
          g12_ = g12(elem%IndexH2Dto3D_bnd(p),ke2d)
          g22_ = g22(elem%IndexH2Dto3D_bnd(p),ke2d)
 
          momw = momz(im) / gsqrtv &
             + g13(im) * momx(im) + g23(im) * momy(im)
          fac = nz(i) * gsqrtv**2 / ( 1.0_rp + g11_ * ( gsqrtv * g13(im) )**2 + 2.0_rp * g12_ * ( gsqrtv**2 * g13(im) * g23(im) ) + g22_ * ( gsqrtv * g23(im) )**2 )
 
          mom_normal = momx(im) * nx(i) + momy(im) * ny(i) + momw * nz(i)
          momx(ip) = momx(im) - 2.0_rp * mom_normal * ( nx(i) + fac * ( g11_ * g13(im) + g12_ * g23(im) ) )
          momy(ip) = momy(im) - 2.0_rp * mom_normal * ( ny(i) + fac * ( g12_ * g13(im) + g22_ * g23(im) ) )
          momz(ip) = momz(im) - 2.0_rp * mom_normal * fac / gsqrtv
        
        case ( bnd_type_noslip_id )
          momx(ip) = - momx(im)
          momy(ip) = - momy(im)
          momz(ip) = - momz(im)          
        end select
 
        select case( this%ThermalBC_list(domid)%list(i_) )
        case ( bnd_type_fixval_id )
          temp_b =  this%ThermalBC_list(domid)%val(i_)
          temp_p = 2.0_rp * temp_b - pres(im) / ( ( dens_hyd(im) + ddens(im) ) * rtot(im) )
          pt(ip) = temp_p * ( pres00 / pres(im) )**(rtot(im)/cptot(im))
          ddens(ip) = pres(im) / ( rtot(im) * temp_p ) - dens_hyd(im)
        end select
      end if
 
    end do
    end do
    
    return
  end  subroutine atmos_dyn_bnd_applybc_tbvars_lc
 
!OCL SERIAL
  subroutine atmos_dyn_bnd_applybc_tbstress_lc( this,  &
    domID,                                              & ! (in)
    t11, t12, t13, t21, t22, t23, t31, t32, t33,        & ! (inout)
    df1, df2, df3,                                      & ! (in)
    gsqrt, gsqrth, g11, g12, g22, g13, g23, nx, ny, nz, & ! (in)
    vmapm, vmapp, vmapb, lmesh, elem, lmesh2d, elem2d   ) ! (in)
 
    use scale_mesh_bndinfo, only: &
      bnd_type_slip_id, bnd_type_noslip_id, &
      bnd_type_adiabat_id, bnd_type_fixval_id
        
    implicit none
 
    class(atmdynbnd), intent(in) :: this    
    integer, intent(in) :: domid
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    class(localmesh2d), intent(in) :: lmesh2d
    class(elementbase2d), intent(in) :: elem2d
    real(rp), intent(inout) :: t11(elem%np*lmesh%nea), t12(elem%np*lmesh%nea), t13(elem%np*lmesh%nea)
    real(rp), intent(inout) :: t21(elem%np*lmesh%nea), t22(elem%np*lmesh%nea), t23(elem%np*lmesh%nea)
    real(rp), intent(inout) :: t31(elem%np*lmesh%nea), t32(elem%np*lmesh%nea), t33(elem%np*lmesh%nea)
    real(rp), intent(inout) :: df1(elem%np*lmesh%nea), df2(elem%np*lmesh%nea), df3(elem%np*lmesh%nea)
    real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
    real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g11(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g12(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g22(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  g13(elem%np*lmesh%nea)
    real(rp), intent(in) ::  g23(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, intent(in) :: vmapb(:)
 
    integer :: p, ke, ke2d
    integer :: i, i_, im, ip
 
    real(rp) :: gsqrtv, g11_, g12_, g22_
    real(rp) :: fac, normal_vec(3)
 
    real(rp) :: stress_t_tmp, stress_t_tmpvec(3)
    real(rp) :: vecw
    real(rp) :: normal_flux
    !-----------------------------------------------
 
    !$omp parallel do collapse(2) private( &
    !$omp ke, p, ke2D, i, i_, iM, iP,                      &
    !$omp VECW, GsqrtV, G11_, G12_, G22_, fac, normal_vec, &
    !$omp stress_t_tmp, stress_t_tmpvec, normal_flux )
    do ke=lmesh%NeS, lmesh%NeE
    do p=1, elem%NfpTot
      i = p + (ke-1)*elem%NfpTot
      ip = vmapp(i)
      i_ = ip - elem%Np*lmesh%NeE
      
      if (i_ > 0) then
        ke2d = lmesh%EMap3Dto2D(ke)
        im = vmapm(i)
 
        gsqrtv = gsqrt(im) / gsqrth(elem%IndexH2Dto3D_bnd(p),ke2d) 
        g11_ = g11(elem%IndexH2Dto3D_bnd(p),ke2d)
        g12_ = g12(elem%IndexH2Dto3D_bnd(p),ke2d)
        g22_ = g22(elem%IndexH2Dto3D_bnd(p),ke2d)
        fac = nz(i) * gsqrtv**2 / ( 1.0_rp + g11_ * ( gsqrtv * g13(im) )**2 + 2.0_rp * g12_ * ( gsqrtv**2 * g13(im) * g23(im) ) + g22_ * ( gsqrtv * g23(im) )**2 )
        normal_vec(1) = nx(i) + fac * ( g11_ * g13(im) + g12_ * g23(im) )
        normal_vec(2) = ny(i) + fac * ( g12_ * g13(im) + g22_ * g23(im) )
        normal_vec(3) = fac / gsqrtv
 
        select case( this%VelBC_list(domid)%list(i_) )
        case ( bnd_type_slip_id)
          vecw = t13(im) / gsqrtv + g13(im) * t11(im) + g23(im) * t12(im)
          stress_t_tmpvec(1) = t11(im) * nx(i) + t12(im) * ny(i) + vecw * nz(i)
 
          vecw = t23(im) / gsqrtv + g13(im) * t21(im) + g23(im) * t22(im)
          stress_t_tmpvec(2) = t21(im) * nx(i) + t22(im) * ny(i) + vecw * nz(i)
 
          vecw = t33(im) / gsqrtv
          stress_t_tmpvec(3) = t31(im) * nx(i) + t32(im) * ny(i) + vecw * nz(i)
 
          stress_t_tmp = sum( stress_t_tmpvec(:) * normal_vec(:) )
          stress_t_tmpvec(:) = stress_t_tmpvec(:) - stress_t_tmp * normal_vec(:) 
 
          ! T1j
          t11(ip) = t11(im) - 2.0_rp * stress_t_tmpvec(1) * normal_vec(1)
          t12(ip) = t12(im) - 2.0_rp * stress_t_tmpvec(1) * normal_vec(2)
          t13(ip) = t13(im) - 2.0_rp * stress_t_tmpvec(1) * normal_vec(3)
 
          ! T2j
          t21(ip) = t21(im) - 2.0_rp * stress_t_tmpvec(2) * normal_vec(1)
          t22(ip) = t22(im) - 2.0_rp * stress_t_tmpvec(2) * normal_vec(2)
          t23(ip) = t23(im) - 2.0_rp * stress_t_tmpvec(2) * normal_vec(3)
 
          ! T3j
          t31(ip) = t31(im) - 2.0_rp * stress_t_tmpvec(3) * normal_vec(1)
          t32(ip) = t32(im) - 2.0_rp * stress_t_tmpvec(3) * normal_vec(2)
          t33(ip) = t33(im) - 2.0_rp * stress_t_tmpvec(3) * normal_vec(3)
 
        case ( bnd_type_noslip_id )
        end select
 
        select case( this%ThermalBC_list(domid)%list(i_) )
        case ( bnd_type_adiabat_id )
          normal_flux = df1(im) * nx(i) + df2(im) * ny(i) &
            + ( df3(im) / gsqrtv + g13(im) * df1(im) + g23(im) * df2(im) ) * nz(i) 
          df1(ip) = df1(im) - 2.0_rp * normal_flux * normal_vec(1)
          df2(ip) = df2(im) - 2.0_rp * normal_flux * normal_vec(2)
          df3(ip) = df3(im) - 2.0_rp * normal_flux * normal_vec(3)
        end select
      end if
    end do
    end do
    
    return
  end  subroutine atmos_dyn_bnd_applybc_tbstress_lc
 
!OCL SERIAL
  subroutine atmos_dyn_bnd_inquire_bound_flag(  this,      & ! (in)
    is_bound,                                              & ! (out)
    domid, vmapm, vmapp, vmapb, lmesh, elem                ) ! (in)
 
    use scale_mesh_bndinfo, only: &
      bnd_type_slip_id, bnd_type_noslip_id, &
      bnd_type_adiabat_id
    implicit none
 
    class(atmdynbnd), intent(in) :: this
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem    
    logical, intent(out) :: is_bound(elem%nfptot*lmesh%ne)
    integer, intent(in) :: domid
    integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
    integer, intent(in) :: vmapb(:)
 
    integer :: i, i_, im, ip
    !-----------------------------------------------
 
    do i=1, elem%NfpTot*lmesh%Ne
      ip = vmapp(i)
      i_ = ip - elem%Np*lmesh%NeE
      is_bound(i) = .false.
 
      if (i_ > 0) then  
        im = vmapm(i)
        if ( this%VelBC_list(domid)%list(i_) == bnd_type_slip_id ) then
          is_bound(i) = .true.
        else if ( this%VelBC_list(domid)%list(i_) == bnd_type_noslip_id ) then          
          is_bound(i) = .true.
        end if
 
      end if
    end do
 
    return
  end subroutine atmos_dyn_bnd_inquire_bound_flag
 
  !-- private -------------------------------------------------------
 
!OCL SERIAL
  subroutine bnd_init_lc(   &
    velBCInfo, thermalBCInfo,           &  ! (inout)
    velbc_ids, thermalbc_ids,           &  ! (in)
    thermal_fixval,                     &  ! (in)
    vmapb, mesh, lmesh, elem )             ! (in)
 
    use scale_mesh_bndinfo, only: bnd_type_nospec_id
    implicit none
    
    type(meshbndinfo), intent(inout) :: velbcinfo
    type(meshbndinfo), intent(inout) :: thermalbcinfo
    integer, intent(in) :: velbc_ids(dom_bnd_num)
    integer, intent(in) :: thermalbc_ids(dom_bnd_num)
    real(rp), intent(in) :: thermal_fixval(dom_bnd_num)
    class(meshbase), intent(in) :: mesh
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: vmapb(:)
 
    integer :: tileid
    integer :: dom_bnd_sizes(dom_bnd_num)
    integer :: bnd_buf_size
    integer :: b, is_, ie_
 
    !-----------------------------------------------
 
    dom_bnd_sizes(:) = &
        elem%Nfp_h*lmesh%NeZ*(/ lmesh%NeX, lmesh%NeY, lmesh%NeX, lmesh%NeY, 0, 0 /) &
      + elem%Nfp_v*lmesh%NeX*lmesh%NeY*(/ 0, 0, 0, 0, 1, 1 /)
    bnd_buf_size = sum(dom_bnd_sizes)
    
    call velbcinfo%Init( bnd_buf_size )
    call velbcinfo%Set(1, bnd_buf_size, bnd_type_nospec_id)
 
    call thermalbcinfo%Init( bnd_buf_size )
    call thermalbcinfo%Set(1, bnd_buf_size, bnd_type_nospec_id)
 
    tileid = lmesh%tileID
    is_ = 1
    do b=1, dom_bnd_num
      ie_ = is_ + dom_bnd_sizes(b) - 1
      if ( mesh%tileID_globalMap(b,tileid) == tileid      &
           .and. mesh%tileFaceID_globalMap(b,tileid) == b ) then
        call velbcinfo%Set( is_, ie_, velbc_ids(b) )
        call thermalbcinfo%Set( is_, ie_, thermalbc_ids(b), thermal_fixval(b) )
      end if
      is_ = ie_ + 1
    end do
 
    return
  end  subroutine bnd_init_lc
 
end module scale_atm_dyn_dgm_bnd