scale_meshutil_vcoord.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
#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 & procedure
  !  
  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 & procedure
  !
  !-----------------------------------------------------------------------------
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !  
  !-----------------------------------------------------------------------------
 
contains
!OCL SERIAL
  function meshutil_get_vcoord_typeid( vcoord_type ) result( vcoord_id )
    implicit none
 
    character(len=*), intent(in) :: vcoord_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
 
!OCL SERIAL
  subroutine meshutil_vcoord_getmetric( G13, G23, zlev, GsqrtV, &
    topo, zTop, vcoord_id, lcmesh, elem, lcmesh2D, elem2D,      &
    Dx2D, Dy2D, Lift2D )
    implicit none
 
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem   
    class(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
    real(rp) :: del_flux(elem2d%nfptot,lcmesh2d%ne,2)
    real(rp) :: fx2d(elem2d%np), fy2d(elem2d%np), liftdelflux2d(elem2d%np)
    real(rp) :: gradzs(elem2d%np,lcmesh2d%ne,2)
    real(rp) :: coef3d(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)
 
      call cal_del_flux( del_flux, &
        topo, lcmesh2d%normal_fn(:,:,1), lcmesh2d%normal_fn(:,:,2), &
        lcmesh2d%VMapM, lcmesh2d%VMapP, lcmesh2d, elem2d            )
      
      !$omp parallel private(ke2D, ke,        &
      !$omp Fx2D, Fy2D, LiftDelFlux2D, coef3D )
 
      !$omp do
      do ke2d=1, lcmesh2d%Ne
        call sparsemat_matmul( dx2d, topo(:,ke2d), fx2d )
        call sparsemat_matmul( lift2d, lcmesh2d%Fscale(:,ke2d) * del_flux(:,ke2d,1), liftdelflux2d)
        gradzs(:,ke2d,1) = lcmesh2d%Escale(:,ke2d,1,1) * fx2d(:) + liftdelflux2d(:)
 
        call sparsemat_matmul( dy2d, topo(:,ke2d), fy2d )
        call sparsemat_matmul( lift2d, lcmesh2d%Fscale(:,ke2d) * del_flux(:,ke2d,2), liftdelflux2d)
        gradzs(:,ke2d,2) = lcmesh2d%Escale(:,ke2d,2,2) * fy2d(:) + liftdelflux2d(:)
      end do
      !$omp end do
 
      !$omp do
      do ke=1, lcmesh%Ne
        ke2d = lcmesh%EMap3Dto2D(ke)
        coef3d(:) = 1.0_rp - lcmesh%pos_en(:,ke,3) / ztop
        zlev(:,ke) = lcmesh%pos_en(:,ke,3)                  &
                   + coef3d(:) * topo(elem%IndexH2Dto3D,ke2d)
 
        gsqrtv(:,ke) = 1.0_rp - topo(elem%IndexH2Dto3D,ke2d) / ztop ! dz/dzeta
        coef3d(:) = - coef3d(:) / gsqrtv(:,ke)
        g13(:,ke) = coef3d(:) * gradzs(elem%IndexH2Dto3D(:),ke2d,1)
        g23(:,ke) = coef3d(:) * gradzs(elem%IndexH2Dto3D(:),ke2d,2)
      end do
      !$omp end do
      !$omp end parallel
    else
      log_error("Mesh_VCoord_GetMetric",*) "vcoord_id is inappropriate. Check!", vcoord_id
      call prc_abort        
    end if
 
    return
  end subroutine meshutil_vcoord_getmetric
 
!OCL SERIAL
  subroutine cal_del_flux( del_flux, &
    topo, nx, ny, vmapM, vmapP, lmesh, elem )
    implicit none
    type(localmesh2d), intent(in) :: lmesh
    type(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 )
    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