scale_meshutil_vcoord.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Mesh / utility for general vertical coordinate
!!
!! @par Description
!!          A module useful for general vertical coordinate
!!
!! @author Yuta Kawai, Team SCALE
!<
#include "scaleFElib.h"
module scale_meshutil_vcoord
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_const, only: &
    pi => const_pi,      &
    eps => const_eps
  use scale_precision
  use scale_prc
  use scale_io
  use scale_prc
 
  use scale_localmesh_2d, only: localmesh2d
  use scale_localmesh_3d, only: localmesh3d
  use scale_element_base, only: &
    elementbase2d, elementbase3d
  use scale_sparsemat, only:&
    sparsemat, sparsemat_matmul
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedures
  !  
  public :: meshutil_get_vcoord_typeid
  public :: meshutil_vcoord_getmetric
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  character(*), public, parameter :: mesh_vcoord_terrain_following_name = "TERRAIN_FOLLOWING"
  integer, public, parameter :: mesh_vcoord_terrain_following_id        = 1
 
  !-----------------------------------------------------------------------------
  !
  !++ Private type & procedures
  !
  !-----------------------------------------------------------------------------
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !  
  !-----------------------------------------------------------------------------
 
contains
  !> Get a type ID of vertical coordinate
!OCL SERIAL
  function meshutil_get_vcoord_typeid( vcoord_type ) result( vcoord_id )
    implicit none
 
    character(len=*), intent(in) :: vcoord_type !< Name of vertical coordinate type
    integer :: vcoord_id
 
    select case( vcoord_type )
    case( mesh_vcoord_terrain_following_name )
      vcoord_id = mesh_vcoord_terrain_following_id
    case default
      log_error("MeshUtil_VCoord_TypeID",*) "vcoord_type is inappropriate. Check!", vcoord_type
      call prc_abort
    end select
 
    return
  end function meshutil_get_vcoord_typeid
 
  !> Get metric terms of vertical coordinate transformation
!OCL SERIAL
  subroutine meshutil_vcoord_getmetric( G13, G23, zlev, GsqrtV, &
    topo, zTop, vcoord_id, lcmesh, elem, lcmesh2D, elem2D,      &
    Dx2D, Dy2D, Lift2D )
    implicit none
 
    type(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem   
    type(localmesh2d), intent(in) :: lcmesh2d
    class(elementbase2d), intent(in) :: elem2d
    real(rp), intent(out) :: g13(elem%np,lcmesh%nea)
    real(rp), intent(out) :: g23(elem%np,lcmesh%nea)
    real(rp), intent(out) :: zlev(elem%np,lcmesh%nea)
    real(rp), intent(inout) :: gsqrtv(elem%np,lcmesh%nea)
    real(rp), intent(in) :: topo(elem2d%np,lcmesh2d%nea)
    integer, intent(in) :: vcoord_id
    real(rp), intent(in) :: ztop
    type(sparsemat), intent(in) :: dx2d
    type(sparsemat), intent(in) :: dy2d
    type(sparsemat), intent(in) :: lift2d
 
    integer :: ke, ke2d, p
    real(rp) :: del_flux(elem2d%nfptot,lcmesh2d%ne,2)
    real(rp) :: fx2d(elem2d%np), fy2d(elem2d%np), liftdelflux2d(elem2d%np,2)
    real(rp) :: gradzs(elem2d%np,lcmesh2d%ne,2)
    real(rp) :: coef3d
 
    integer :: indexh2dto3d(elem%np)
    !------------------------------------------------
 
    if ( vcoord_id == mesh_vcoord_terrain_following_id ) then
      
      ! * z    = topo + (1 - topo / zTop ) * zeta
      ! * zeta = zTop * (z - topo)/(zTop - topo)
      ! * Gi3  = (dzeta(x1,x2,z)/dxi)_z = d (zeta,z) / d (xi,z) = - d(z,zeta)/ d(xi,zeta) *  d(xi,zeta)/d(xi,z)
      !                                 = - (dz/dxi)_zeta * dzeta/dz
      !                                 = (GsqrtV)^-1 * [ - 1 + zeta / zTop ] * d topo /dxi     (i=1, 2)
 
 
      indexh2dto3d(:) = elem%IndexH2Dto3D(:)
      
      !$acc data create( del_flux, GradZs ) copyin( IndexH2Dto3D )
      
      call cal_del_flux( del_flux, &
        topo, lcmesh2d%normal_fn(:,:,1), lcmesh2d%normal_fn(:,:,2), &
        lcmesh2d%VMapM, lcmesh2d%VMapP, lcmesh2d, elem2d            )
      
      !$omp parallel private(Fx2D, Fy2D, LiftDelFlux2D, coef3D )
      !$acc parallel loop gang private(ke2D, Fx2D, Fy2D, LiftDelFlux2D, coef3D ) &
      !$acc present(topo, del_flux, zlev, GsqrtV, G13, G23, lcmesh,elem,lcmesh2D,elem2D)
      !$omp do
      do ke2d=1, lcmesh2d%Ne
        call sparsemat_matmul( dx2d, topo(:,ke2d), fx2d )
        call sparsemat_matmul( dy2d, topo(:,ke2d), fy2d )
#ifdef _OPENACC
        call sparsemat_matmul( lift2d, lcmesh2d%Fscale(:,ke2d), del_flux(:,ke2d,1), liftdelflux2d(:,1))
        call sparsemat_matmul( lift2d, lcmesh2d%Fscale(:,ke2d), del_flux(:,ke2d,2), liftdelflux2d(:,2))
#else
        call sparsemat_matmul( lift2d, lcmesh2d%Fscale(:,ke2d) * del_flux(:,ke2d,1), liftdelflux2d(:,1))
        call sparsemat_matmul( lift2d, lcmesh2d%Fscale(:,ke2d) * del_flux(:,ke2d,2), liftdelflux2d(:,2))
#endif
        !$acc loop vector
        do p=1, elem2d%Np
          gradzs(p,ke2d,1) = lcmesh2d%Escale(p,ke2d,1,1) * fx2d(p) + liftdelflux2d(p,1)
          gradzs(p,ke2d,2) = lcmesh2d%Escale(p,ke2d,2,2) * fy2d(p) + liftdelflux2d(p,2)
        end do
      end do
      !$omp end do
 
      !$omp do
      !$acc parallel loop gang private(ke2D,coef3D) present(lcmesh,elem)
      do ke=1, lcmesh%Ne
        ke2d = lcmesh%EMap3Dto2D(ke)
        !$acc loop vector
        do p=1, elem%Np
          coef3d = 1.0_rp - lcmesh%pos_en(p,ke,3) / ztop
          zlev(p,ke) = lcmesh%pos_en(p,ke,3)                 &
                     + coef3d * topo(indexh2dto3d(p),ke2d)
 
          gsqrtv(p,ke) = 1.0_rp - topo(indexh2dto3d(p),ke2d) / ztop ! dz/dzeta
 
          coef3d = - coef3d / gsqrtv(p,ke)
          g13(p,ke) = coef3d * gradzs(indexh2dto3d(p),ke2d,1)
          g23(p,ke) = coef3d * gradzs(indexh2dto3d(p),ke2d,2)
        end do
      end do
      !$omp end do
      !$acc end parallel
      !$omp end parallel
 
      !$acc end data      
    else
      log_error("Mesh_VCoord_GetMetric",*) "vcoord_id is inappropriate. Check!", vcoord_id
      call prc_abort        
    end if
    return
  end subroutine meshutil_vcoord_getmetric
 
!- private subroutines ----------------------
 
!OCL SERIAL
  subroutine cal_del_flux( del_flux, &
    topo, nx, ny, vmapM, vmapP, lmesh, elem )
    implicit none
    type(localmesh2d), intent(in) :: lmesh
    class(elementbase2d), intent(in) :: elem
    real(rp), intent(out) :: del_flux(elem%nfptot*lmesh%ne,2)
    real(rp), intent(in) :: topo(elem%np*lmesh%nea)
    real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
    real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
    integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
 
    integer :: i
    integer :: ip, im
    real(rp) :: dtopo
    !-------------------------------------
 
    !$omp parallel do private( i, iM, iP, dtopo )
    !$acc parallel loop present(topo, nx, ny, vmapM, vmapP, del_flux, lmesh, elem)
    do i=1, elem%NfpTot * lmesh%Ne
      im = vmapm(i); ip = vmapp(i)
 
      dtopo = 0.5_rp * ( topo(ip) - topo(im) )
 
      del_flux(i,1) = dtopo * nx(i)
      del_flux(i,2) = dtopo * ny(i)
    end do
    return
  end subroutine cal_del_flux
end module scale_meshutil_vcoord