module FElib / Mesh / utility for general vertical coordinate More...

Functions/Subroutines
integer function, public	meshutil_get_vcoord_typeid (vcoord_type)
	Get a type ID of vertical coordinate.
subroutine, public	meshutil_vcoord_getmetric (g13, g23, zlev, gsqrtv, topo, ztop, vcoord_id, lcmesh, elem, lcmesh2d, elem2d, dx2d, dy2d, lift2d)
	Get metric terms of vertical coordinate transformation.

Variables
character(*), parameter, public	mesh_vcoord_terrain_following_name = "TERRAIN_FOLLOWING"
integer, parameter, public	mesh_vcoord_terrain_following_id = 1

Detailed Description

module FElib / Mesh / utility for general vertical coordinate

Description: A module useful for general vertical coordinate

Author: Yuta Kawai, Team SCALE

Function/Subroutine Documentation

◆ meshutil_get_vcoord_typeid()

integer function, public scale_meshutil_vcoord::meshutil_get_vcoord_typeid ( character(len=*), intent(in) vcoord_type )

Get a type ID of vertical coordinate.

Parameters

[in] vcoord_type Name of vertical coordinate type

Definition at line 62 of file scale_meshutil_vcoord.F90.

    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

References mesh_vcoord_terrain_following_id, and mesh_vcoord_terrain_following_name.

Referenced by mod_atmos_mesh_gm::atmosmeshgm_init(), mod_atmos_mesh_rm::atmosmeshrm_init(), mod_ocean_mesh_gm::oceanmeshgm_init(), and mod_ocean_mesh_rm::oceanmeshrm_init().

◆ meshutil_vcoord_getmetric()

subroutine, public scale_meshutil_vcoord::meshutil_vcoord_getmetric	(	real(rp), dimension(elem%np,lcmesh%nea), intent(out)	g13,
		real(rp), dimension(elem%np,lcmesh%nea), intent(out)	g23,
		real(rp), dimension(elem%np,lcmesh%nea), intent(out)	zlev,
		real(rp), dimension(elem%np,lcmesh%nea), intent(inout)	gsqrtv,
		real(rp), dimension(elem2d%np,lcmesh2d%nea), intent(in)	topo,
		real(rp), intent(in)	ztop,
		integer, intent(in)	vcoord_id,
		type(localmesh3d), intent(in)	lcmesh,
		class(elementbase3d), intent(in)	elem,
		type(localmesh2d), intent(in)	lcmesh2d,
		class(elementbase2d), intent(in)	elem2d,
		type(sparsemat), intent(in)	dx2d,
		type(sparsemat), intent(in)	dy2d,
		type(sparsemat), intent(in)	lift2d )

Get metric terms of vertical coordinate transformation.

Definition at line 81 of file scale_meshutil_vcoord.F90.

    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