scale_mesh_topography.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / Topography
!!
!! @par Description
!!          A module to manage vertical coordinate with topography 
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_mesh_topography
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
 
  use scale_localmesh_3d, only: localmesh3d
  use scale_localmesh_2d, only: localmesh2d
  use scale_mesh_base2d, only: meshbase2d
  use scale_mesh_base3d, only: meshbase3d
  use scale_element_base, only: elementbase3d
  use scale_meshfield_base, only: &
    meshfield2d, meshfield3d
  use scale_meshfieldcomm_base, only: &
    meshfieldcommbase, meshfieldcontainer
 
  use scale_sparsemat, only: sparsemat
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
  !> Derived type to manage topography datat and setup vertical coordinate metric
  type, public :: meshtopography
    type(meshfield2d) :: topo 
  contains
    procedure :: init => meshtopography_init
    procedure :: final => meshtopography_final
    procedure :: setvcoordinate => meshtopography_set_vcoordinate
  end type meshtopography
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
!OCL SERIAL
  subroutine meshtopography_init( this, varname, mesh )
    implicit none
 
    class(meshtopography), intent(inout) :: this
    character(len=*), intent(in) :: varname
    class(meshbase2d), intent(in), target :: mesh
 
    integer :: n
    !-----------------------------------------------------------------------------
 
    call this%topo%Init( varname, "m", mesh )
 
    do n=1, mesh%LOCAL_MESH_NUM
!$omp parallel workshare
      this%topo%local(n)%val(:,:) = 0.0_rp
!$omp end parallel workshare
      !$acc update device(this%topo%local(n)%val)
    end do
 
    return
  end subroutine meshtopography_init
 
  !> Finalize an object of MeshTopography
!OCL SERIAL
  subroutine meshtopography_final( this )
    implicit none
 
    class(meshtopography), intent(inout) :: this
    !-----------------------------------------------------------------------------
          
    call this%topo%Final()
    return
  end subroutine meshtopography_final
  
  !> Setup vertical coordinate metric with topography data
!OCL SERIAL
  subroutine meshtopography_set_vcoordinate( this, mesh3D, &
      vcoord_id, zTop, comm3D, comm2D                      )
    use scale_meshfieldcomm_cubedspheredom3d, only: meshfieldcommcubedspheredom3d
    implicit none
 
    class(meshtopography), target, intent(inout) :: this   !< Object of MeshTopography to manage topography data and vertical coordinate metric
    class(meshbase3d), intent(inout), target :: mesh3D     !< 3D mesh object
    integer, intent(in) :: vcoord_id                       !< Vertical coordinate ID
    real(RP), intent(in) :: zTop                           !< Height of the model domain
    class(meshfieldcommbase), intent(inout) :: comm3D      !< Object for 3D data communication
    class(meshfieldcommbase), intent(inout) :: comm2D      !< Object for 2D data communication
 
    type(sparsemat) :: Dx2D, Dy2D, Lift2D
    class(localmesh3d), pointer :: lcmesh
    class(localmesh2d), pointer :: lcmesh2D
    class(elementbase3d), pointer :: elem3D
 
    integer :: n
    integer :: ke, ke2D
 
    type(meshfieldcontainer) :: comm2d_varlist(1)
    type(meshfieldcontainer) :: comm3d_varlist(4)
 
    type(meshfield3d), target :: zlev, GsqrtV, G13, G23
    type(meshfield3d), target :: tmp_G13, tmp_G23
 
    logical :: flag_covariantvec
    !-------------------------------------------------------------
 
    lcmesh => mesh3d%lcmesh_list(1)
    lcmesh2d => lcmesh%lcmesh2D
    call dx2d  %Init( lcmesh2d%refElem2D%Dx1  )
    call dy2d  %Init( lcmesh2d%refElem2D%Dx2  )
    call lift2d%Init( lcmesh2d%refElem2D%Lift )
 
    ! Exchange topography data to fill halo
 
    comm2d_varlist(1)%field2d => this%topo
    call comm2d%Put(comm2d_varlist, 1)
    call comm2d%Exchange()
    call comm2d%Get(comm2d_varlist, 1)
 
    ! Calculate metric factors associated with general vertical coordinates
 
    call zlev%Init( "zlev", "", mesh3d )
    call gsqrtv%Init( "GsqrtV", "", mesh3d )
    call g13%Init( "G13", "", mesh3d )
    call g23%Init( "G23", "", mesh3d )
 
    call tmp_g13%Init( "tmp_G13", "", mesh3d )
    call tmp_g23%Init( "tmp_G23", "", mesh3d )
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      lcmesh => mesh3d%lcmesh_list(n)
      lcmesh2d => lcmesh%lcmesh2D
      
      call calc_vcoordinate_metrics_lc( &
        g13%local(n)%val, g23%local(n)%val, tmp_g13%local(n)%val, tmp_g23%local(n)%val,     & ! (out)
        zlev%local(n)%val, gsqrtv%local(n)%val,                                             & ! (out)
        this%topo%local(n)%val, ztop, vcoord_id,                                            & ! (in)
        lcmesh, lcmesh%refElem3D, lcmesh2d, lcmesh2d%refElem2D,                             & ! (in)
        lcmesh%GIJ(:,:,1,1), lcmesh%GIJ(:,:,1,2), lcmesh%GIJ(:,:,2,1), lcmesh%GIJ(:,:,2,2), & ! (in)
        dx2d, dy2d, lift2d )                                                                  ! (in)
    end do
 
    ! Exchange metric data to fill halo
 
    comm3d_varlist(1)%field3d => zlev
    comm3d_varlist(2)%field3d => gsqrtv
    comm3d_varlist(3)%field3d => tmp_g13
    comm3d_varlist(4)%field3d => tmp_g23
 
    flag_covariantvec = .false.
    select type(comm3d)
    type is (meshfieldcommcubedspheredom3d)
      call comm3d%SetCovariantVec( 1, g13, g23 )
      flag_covariantvec = .true.
    end select
 
    call comm3d%Put(comm3d_varlist, 1)
    call comm3d%Exchange()
    call comm3d%Get(comm3d_varlist, 1)
 
    ! Update metric data managed by local mesh
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      lcmesh => mesh3d%lcmesh_list(n)
      elem3d => lcmesh%refElem3D
 
      if ( flag_covariantvec ) then
        call mesh3d%Set_geometric_with_vcoord( n, &
          gsqrtv%local(n)%val, zlev%local(n)%val, g13%local(n)%val, g23%local(n)%val )
      else
        call mesh3d%Set_geometric_with_vcoord( n, &
          gsqrtv%local(n)%val, zlev%local(n)%val, tmp_g13%local(n)%val, tmp_g23%local(n)%val )
      end if
    end do
 
    !---
    call zlev%Final()
    call gsqrtv%Final()
    call g13%Final()
    call g23%Final()
    call tmp_g13%Final()
    call tmp_g23%Final()
 
    call dx2d%Final()
    call dy2d%Final()
    call lift2d%Final()
 
    return
  end subroutine meshtopography_set_vcoordinate
 
!-- prviate ----
 
!OCL SERIAL
  subroutine calc_vcoordinate_metrics_lc( G13, G23, tmp_G13, tmp_G23, zlev, GsqrtV, &
    topo, zTop, vcoord_id, lcmesh, elem3D, lcmesh2D, elem2D,      &
    G11, G12, G21, G22, Dx2D, Dy2D, Lift2D )
    use scale_element_base, only: elementbase2d
    use scale_meshutil_vcoord, only: meshutil_vcoord_getmetric
    implicit none
    type(localmesh3d), intent(in) :: lcmesh
    type(localmesh2d), intent(in) :: lcmesh2D
    class(elementbase3d), intent(in) :: elem3D
    class(elementbase2d), intent(in) :: elem2D
    real(RP), intent(out) :: G13(elem3D%Np,lcmesh%NeA)
    real(RP), intent(out) :: G23(elem3D%Np,lcmesh%NeA)
    real(RP), intent(out) :: tmp_G13(elem3D%Np,lcmesh%NeA)
    real(RP), intent(out) :: tmp_G23(elem3D%Np,lcmesh%NeA)
    real(RP), intent(out) :: zlev(elem3D%Np,lcmesh%NeA)
    real(RP), intent(out) :: GsqrtV(elem3D%Np,lcmesh%NeA)
    real(RP), intent(in) :: topo(elem2D%Np,lcmesh2D%NeA)
    real(RP), intent(in) :: zTop
    integer, intent(in) :: vcoord_id
    real(RP), intent(in) :: G11(elem2D%Np,lcmesh2D%Ne)
    real(RP), intent(in) :: G12(elem2D%Np,lcmesh2D%Ne)
    real(RP), intent(in) :: G21(elem2D%Np,lcmesh2D%Ne)
    real(RP), intent(in) :: G22(elem2D%Np,lcmesh2D%Ne)
    type(sparsemat), intent(in) :: Dx2D, Dy2D, Lift2D
 
    integer :: IndexH2Dto3D(elem3D%Np)
    integer :: ke, ke2D, p
    !------------------------------------------------------------------------------
 
    call meshutil_vcoord_getmetric( &
      g13, g23, zlev, gsqrtv,                                 & ! (out)
      topo, ztop, vcoord_id,                                  & ! (in)
      lcmesh, lcmesh%refElem3D, lcmesh2d, lcmesh2d%refElem2D, & ! (in)
      dx2d, dy2d, lift2d                                      ) ! (in)
 
    indexh2dto3d(:) = elem3d%IndexH2Dto3D(:)
    !$omp parallel do private(ke2D)
    !$acc parallel loop private(ke2D) collapse(2) &
    !$acc   present(G11, G12, G21, G22, G13, G23, tmp_G13, tmp_G23, lcmesh, elem3D) copyin(IndexH2Dto3D)
    do ke=lcmesh%NeS, lcmesh%NeE
    do p=1, elem3d%Np
      ke2d = lcmesh%EMap3Dto2D(ke)
      tmp_g13(p,ke) = g11(indexh2dto3d(p),ke2d) * g13(p,ke) + g12(indexh2dto3d(p),ke2d) * g23(p,ke)
      tmp_g23(p,ke) = g21(indexh2dto3d(p),ke2d) * g13(p,ke) + g22(indexh2dto3d(p),ke2d) * g23(p,ke)
    end do
    end do
    return
  end subroutine calc_vcoordinate_metrics_lc
 
end module scale_mesh_topography