mod_mkinit_util.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module Utility for mkinit
!!
!! @par Description
!!          subroutines useful to prepare initial data
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module mod_mkinit_util
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_prc
 
  use scale_const, only: &
    pi => const_pi
 
  use scale_element_base, only: elementbase2d, elementbase3d
  use scale_element_quadrilateral, only: quadrilateralelement
  use scale_element_hexahedral, only: hexahedralelement
  use scale_localmesh_2d, only: localmesh2d
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_cubedom3d, only: meshcubedom3d  
  use scale_localmeshfield_base, only: localmeshfieldbase
  
    
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: mkinitutil_gen_gpmat 
  public :: mkinitutil_gen_vm1mat 
  public :: mkinitutil_calc_cosinebell
  public :: mkinitutil_calc_cosinebell_global
  public :: mkinitutil_galerkinprojection
  public :: mkinitutil_galerkinprojection_global
 
contains
 
!OCL SERIAL
  subroutine mkinitutil_gen_gpmat( GPMat, &
    elem_intrp, elem )
    implicit none
 
    class(elementbase3d), intent(in) :: elem_intrp
    class(elementbase3d), intent(in) :: elem
    real(rp), intent(out) :: gpmat(elem%np,elem_intrp%np)
 
    integer :: p1, p2, p3, p_
    integer :: p_intrp
 
    real(rp) :: invv_intrp(elem%np,elem_intrp%np)
    !---------------------------------------------
 
    invv_intrp(:,:) = 0.0_rp
    do p3=1, elem%PolyOrder_v+1
    do p2=1, elem%PolyOrder_h+1
    do p1=1, elem%PolyOrder_h+1
      p_ = p1 + (p2-1)*(elem%PolyOrder_h + 1) + (p3-1)*(elem%PolyOrder_h + 1)**2
      p_intrp = p1 + (p2-1)*(elem_intrp%PolyOrder_h + 1) + (p3-1)*(elem_intrp%PolyOrder_h + 1)**2
      invv_intrp(p_,:) = elem_intrp%invV(p_intrp,:)
    end do
    end do
    end do
    gpmat(:,:) = matmul(elem%V, invv_intrp)
 
    return
  end subroutine mkinitutil_gen_gpmat
 
!OCL SERIAL  
  subroutine mkinitutil_gen_vm1mat( Vm1Mat, &
    elem_intrp, elem )
    implicit none
 
    class(elementbase3d), intent(in) :: elem_intrp
    class(elementbase3d), intent(in) :: elem
    real(rp), intent(out) :: vm1mat(elem%np,elem_intrp%np)
 
    integer :: p1, p2, p3, p_
    integer :: p_intrp
 
    real(rp) :: invv_intrpvm1(elem%np,elem_intrp%np)
    !---------------------------------------------
 
    invv_intrpvm1(:,:) = 0.0_rp
    do p3=1, elem%PolyOrder_v
    do p2=1, elem%PolyOrder_h+1
    do p1=1, elem%PolyOrder_h+1
      p_ = p1 + (p2-1)*(elem%PolyOrder_h + 1) + (p3-1)*(elem%PolyOrder_h + 1)**2
      p_intrp = p1 + (p2-1)*(elem_intrp%PolyOrder_h + 1) + (p3-1)*(elem_intrp%PolyOrder_h + 1)**2
      invv_intrpvm1(p_,:) = elem_intrp%invV(p_intrp,:)
    end do
    end do
    end do
    vm1mat(:,:) = matmul(elem%V, invv_intrpvm1)
 
    return
  end subroutine mkinitutil_gen_vm1mat
 
  !> Calculate the distribution function of a cosine bell in regional domain
  !! 
  !! If the vertical dependence is considered, specify z_func_type and z_func_params.
  !!  For z_func_type = 'sin', the values of z_func_params is
  !!   1:  the vertical model, 2: the half of wavelength
  !! 
!OCL SERIAL  
  subroutine mkinitutil_calc_cosinebell( &
    q,                                           &
    qmax, rx, ry, rz, xc, yc, zc,                &
    x, y, z, lcmesh3D, elem,                     &
    IntrpPolyOrder_h, IntrpPolyOrder_v,          &
    z_func_type, z_func_params, cosbell_exponent )
 
    implicit none
 
    class(localmesh3d), intent(in) :: lcmesh3d
    class(elementbase3d), intent(in) :: elem
    real(rp), intent(out) :: q(elem%np,lcmesh3d%nea)
    real(rp), intent(in) :: qmax
    real(rp), intent(in) :: rx, ry, rz
    real(rp), intent(in) :: xc, yc, zc
    real(rp), intent(in) :: x(elem%np,lcmesh3d%ne)
    real(rp), intent(in) :: y(elem%np,lcmesh3d%ne)
    real(rp), intent(in) :: z(elem%np,lcmesh3d%ne)
    integer, intent(in) :: intrppolyorder_h
    integer, intent(in) :: intrppolyorder_v
    character(len=*), optional, intent(in) :: z_func_type
    real(rp), optional, intent(in) :: z_func_params(:)
    integer, intent(in), optional :: cosbell_exponent  !< parameter to ensure 2*cosbell_exponent-1 continuous derivatives
 
    integer :: ke
 
    type(hexahedralelement) :: elem_intrp
    real(rp), allocatable :: x_intrp(:,:), y_intrp(:,:), z_intrp(:,:)
    real(rp), allocatable :: r_intrp(:)
    real(rp), allocatable :: z_func(:,:)
    real(rp) :: vx(elem%nv), vy(elem%nv), vz(elem%nv)
 
    real(rp), allocatable :: intrpmat(:,:)
    real(rp), allocatable :: q_intrp(:)
 
    integer :: exponent    
    !-----------------------------------------------
 
    if ( present(cosbell_exponent) ) then
      exponent = cosbell_exponent
    else
      exponent = 1
    end if
 
    call elem_intrp%Init( intrppolyorder_h, intrppolyorder_v, .false. )
 
    allocate( intrpmat(elem%Np,elem_intrp%Np) )
    call mkinitutil_gen_gpmat( intrpmat, elem_intrp, elem )
 
    allocate( x_intrp(elem_intrp%Np,lcmesh3d%Ne), y_intrp(elem_intrp%Np,lcmesh3d%Ne), z_intrp(elem_intrp%Np,lcmesh3d%Ne) )
    allocate( z_func(elem_intrp%Np,lcmesh3d%Ne))
    allocate( r_intrp(elem_intrp%Np) )
    allocate( q_intrp(elem_intrp%Np) )
 
    !$omp parallel private( &
    !$omp q_intrp, vx, vy, vz,               &
    !$omp x_intrp, y_intrp, z_intrp, r_intrp )
 
    !$omp do
    do ke=lcmesh3d%NeS, lcmesh3d%NeE
 
      vx(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),1)
      vy(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),2)
      vz(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),3)
      x_intrp(:,ke) = vx(1) + 0.5_rp * ( elem_intrp%x1(:) + 1.0_rp ) * ( vx(2) - vx(1) ) 
      y_intrp(:,ke) = vy(1) + 0.5_rp * ( elem_intrp%x2(:) + 1.0_rp ) * ( vy(4) - vy(1) )
      z_intrp(:,ke) = vz(1) + 0.5_rp * ( elem_intrp%x3(:) + 1.0_rp ) * ( vz(5) - vz(1) )
 
      z_func(:,ke) = 1.0_rp
    end do
    !$omp end do
 
    ! Calculate the vertical function
    if ( present(z_func_type) ) then
      select case(z_func_type)
      case ('sin')
        !$omp do
        do ke=lcmesh3d%NeS, lcmesh3d%NeE
          z_func(:,ke) = sin( z_func_params(1) * pi * z_intrp(:,ke) / z_func_params(2) )
        end do
        !$omp end do
      end select
    end if
 
    !$omp do
    do ke=lcmesh3d%NeS, lcmesh3d%NeE
      r_intrp(:) = sqrt( &
          ( (x_intrp(:,ke) - xc) / rx )**2 &
        + ( (y_intrp(:,ke) - yc) / ry )**2 &
        + ( (z_intrp(:,ke) - zc) / rz )**2 )
      
      where( r_intrp(:) <= 1.0_rp ) 
        q_intrp(:) = qmax * ( 0.5_rp * (1.0_rp + cos( pi * r_intrp(:) ) ) )**exponent
      elsewhere
        q_intrp(:) = 0.0_rp
      end where
      q(:,ke) = matmul(intrpmat, q_intrp(:) * z_func(:,ke))
    end do
    !$omp end do
    !$omp end parallel
 
    call elem_intrp%Final()
 
    return
  end subroutine mkinitutil_calc_cosinebell
 
  !> Calculate the distribution function of a cosine bell in global domain
  !! 
  !! If the vertical dependence is considered, specify z_func_type and z_func_params.
  !!  For z_func_type = 'sin', the values of z_func_params is
  !!   1:  the vertical model, 2: the half of wavelength
  !! 
!OCL SERIAL
  subroutine mkinitutil_calc_cosinebell_global( &
    q,                                           & 
    qmax, rh, lonc, latc, rplanet,               &
    x, y, z, lcmesh3D, elem,                     &
    IntrpPolyOrder_h, IntrpPolyOrder_v,          &
    z_func_type, z_func_params, cosbell_exponent )
 
    use scale_cubedsphere_coord_cnv, only: &
      cubedspherecoordcnv_cs2lonlatpos
    implicit none
 
    class(localmesh3d), intent(in) :: lcmesh3d
    class(elementbase3d), intent(in) :: elem
    real(rp), intent(out) :: q(elem%np,lcmesh3d%nea)
    real(rp), intent(in) :: qmax
    real(rp), intent(in) :: rh
    real(rp), intent(in) :: lonc, latc
    real(rp), intent(in) :: rplanet
    real(rp), intent(in) :: x(elem%np,lcmesh3d%ne)
    real(rp), intent(in) :: y(elem%np,lcmesh3d%ne)
    real(rp), intent(in) :: z(elem%np,lcmesh3d%ne)
    integer, intent(in) :: intrppolyorder_h
    integer, intent(in) :: intrppolyorder_v
    character(len=*), optional, intent(in) :: z_func_type
    real(rp), optional, intent(in) :: z_func_params(:)
    integer, intent(in), optional :: cosbell_exponent      !< parameter to ensure 2*cosbell_exponent-1 continuous derivatives
 
    integer :: ke
 
    type(hexahedralelement) :: elem_intrp
    real(rp), allocatable :: x_intrp(:,:), y_intrp(:,:), z_intrp(:,:)
    real(rp), allocatable :: gam_intrp(:,:)
    real(rp), allocatable :: lon_intrp(:,:), lat_intrp(:,:)
    real(rp), allocatable :: z_func(:,:)
    real(rp), allocatable :: r_intrp(:)
    real(rp) :: vx(elem%nv), vy(elem%nv), vz(elem%nv)
 
    real(rp), allocatable :: intrpmat(:,:)
    real(rp), allocatable :: q_intrp(:)
 
    integer :: exponent
 
    !-----------------------------------------------
 
    if ( present(cosbell_exponent) ) then
      exponent = cosbell_exponent
    else
      exponent = 1
    end if
 
    call elem_intrp%Init( intrppolyorder_h, intrppolyorder_v, .false. )
 
    allocate( intrpmat(elem%Np,elem_intrp%Np) )
    call mkinitutil_gen_gpmat( intrpmat, elem_intrp, elem )
 
    allocate( x_intrp(elem_intrp%Np,lcmesh3d%Ne), y_intrp(elem_intrp%Np,lcmesh3d%Ne), z_intrp(elem_intrp%Np,lcmesh3d%Ne) )
    allocate( gam_intrp(elem_intrp%Np,lcmesh3d%Ne) )
    allocate( lon_intrp(elem_intrp%Np,lcmesh3d%Ne), lat_intrp(elem_intrp%Np,lcmesh3d%Ne) )
    allocate( z_func(elem_intrp%Np,lcmesh3d%Ne) )
    allocate( r_intrp(elem_intrp%Np) )
    allocate( q_intrp(elem_intrp%Np) )
 
 
    !$omp parallel do private(vx, vy, vz)
    do ke=lcmesh3d%NeS, lcmesh3d%NeE
      vx(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),1)
      vy(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),2)
      vz(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),3)
      x_intrp(:,ke) = vx(1) + 0.5_rp * ( elem_intrp%x1(:) + 1.0_rp ) * ( vx(2) - vx(1) ) 
      y_intrp(:,ke) = vy(1) + 0.5_rp * ( elem_intrp%x2(:) + 1.0_rp ) * ( vy(4) - vy(1) )
      z_intrp(:,ke) = vz(1) + 0.5_rp * ( elem_intrp%x3(:) + 1.0_rp ) * ( vz(5) - vz(1) )
      gam_intrp(:,ke) = 1.0_rp
 
      z_func(:,ke) = 1.0_rp
    end do
 
    call cubedspherecoordcnv_cs2lonlatpos( lcmesh3d%panelID, x_intrp, y_intrp, gam_intrp, & ! (in)
      elem_intrp%Np * lcmesh3d%Ne,                                                        & ! (in)
      lon_intrp(:,:), lat_intrp(:,:) )                                                      ! (out)
 
    ! Calculate the vertical function
    if ( present(z_func_type) ) then
      select case(z_func_type)
      case ('sin')
        !$omp parallel do
        do ke=lcmesh3d%NeS, lcmesh3d%NeE
          z_func(:,ke) = sin( z_func_params(1) * pi * z_intrp(:,ke) / z_func_params(2) )
        end do
      end select
    end if
 
    !$omp parallel do private( r_intrp, q_intrp )
    do ke=lcmesh3d%NeS, lcmesh3d%NeE
 
      ! Calculate the horizontal function
      r_intrp(:) = rplanet / rh * acos( sin(latc) * sin(lat_intrp(:,ke)) + cos(latc) * cos(lat_intrp(:,ke)) * cos(lon_intrp(:,ke) - lonc) )
      where( r_intrp(:) <= 1.0_rp ) 
        q_intrp(:) = qmax * ( 0.5_rp * (1.0_rp + cos( pi * r_intrp(:) ) ) )**exponent
      elsewhere
        q_intrp(:) = 0.0_rp
      end where
 
      ! Perform Galerkin projection
      q(:,ke) = matmul(intrpmat, q_intrp(:) * z_func(:,ke))
    end do
 
    call elem_intrp%Final()
 
    return
  end subroutine mkinitutil_calc_cosinebell_global
 
 
!> Apply the Galerkin projection to the user-defined function 
!! 
!OCL SERIAL  
  subroutine mkinitutil_galerkinprojection( q, &
    func, IntrpPolyOrder_h, IntrpPolyOrder_v,  &
    lcmesh3D, elem                             )
  
  implicit none
  class(localmesh3d), intent(in) :: lcmesh3d
  class(elementbase3d), intent(in) :: elem
  real(rp), intent(out) :: q(elem%np,lcmesh3d%nea)
  integer, intent(in) :: intrppolyorder_h
  integer, intent(in) :: intrppolyorder_v
 
  interface
    subroutine func( q_intrp, &
        x, y, z, elem_intrp   )
      import elementbase3d
      import rp
      class(elementbase3d), intent(in) :: elem_intrp
      real(rp), intent(out) :: q_intrp(elem_intrp%np)
      real(rp), intent(in) :: x(elem_intrp%np)
      real(rp), intent(in) :: y(elem_intrp%np)
      real(rp), intent(in) :: z(elem_intrp%np)
    end subroutine func
  end interface
 
  type(hexahedralelement) :: elem_intrp
  real(rp), allocatable :: x_intrp(:,:), y_intrp(:,:), z_intrp(:,:)
  real(rp) :: vx(elem%nv), vy(elem%nv), vz(elem%nv)
 
  real(rp), allocatable :: intrpmat(:,:)
  real(rp), allocatable :: q_intrp(:)
 
  integer :: ke
  !-----------------------------------------------
 
  call elem_intrp%Init( intrppolyorder_h, intrppolyorder_v, .false. )
 
  allocate( intrpmat(elem%Np,elem_intrp%Np) )
  call mkinitutil_gen_gpmat( intrpmat, elem_intrp, elem )
 
  allocate( x_intrp(elem_intrp%Np,lcmesh3d%Ne), y_intrp(elem_intrp%Np,lcmesh3d%Ne), z_intrp(elem_intrp%Np,lcmesh3d%Ne) )
  allocate( q_intrp(elem_intrp%Np) )
 
  !$omp parallel do private(vx, vy, vz)
  do ke=lcmesh3d%NeS, lcmesh3d%NeE
    vx(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),1)
    vy(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),2)
    vz(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),3)
    x_intrp(:,ke) = vx(1) + 0.5_rp * ( elem_intrp%x1(:) + 1.0_rp ) * ( vx(2) - vx(1) ) 
    y_intrp(:,ke) = vy(1) + 0.5_rp * ( elem_intrp%x2(:) + 1.0_rp ) * ( vy(4) - vy(1) )
    z_intrp(:,ke) = vz(1) + 0.5_rp * ( elem_intrp%x3(:) + 1.0_rp ) * ( vz(5) - vz(1) )
  end do
 
  !$omp parallel do private( q_intrp )
  do ke=lcmesh3d%NeS, lcmesh3d%NeE
 
    call func( q_intrp,                                & ! (out)
      x_intrp(:,ke), y_intrp(:,ke), z_intrp(:,ke),     & ! (in)
      elem_intrp                                       ) ! (in)
    
    ! Perform Galerkin projection
    q(:,ke) = matmul( intrpmat, q_intrp )
  end do
 
  call elem_intrp%Final()
 
  return
end subroutine mkinitutil_galerkinprojection
 
!> Apply the Galerkin projection to the user-defined function (for global model)
!!
!OCL SERIAL  
  subroutine mkinitutil_galerkinprojection_global( q, &
      func, IntrpPolyOrder_h, IntrpPolyOrder_v,       &
      lcmesh3D, elem, rplanet )
 
    use scale_cubedsphere_coord_cnv, only: &
      cubedspherecoordcnv_cs2lonlatpos
    
    implicit none
    class(localmesh3d), intent(in) :: lcmesh3d
    class(elementbase3d), intent(in) :: elem
    real(rp), intent(out) :: q(elem%np,lcmesh3d%nea)
    integer, intent(in) :: intrppolyorder_h
    integer, intent(in) :: intrppolyorder_v
    real(rp), intent(in) :: rplanet
 
    interface
      subroutine func( q_intrp,             &
          lon, lat, z, elem_intrp, rplanet_ )
        import elementbase3d
        import rp
        class(elementbase3d), intent(in) :: elem_intrp
        real(rp), intent(out) :: q_intrp(elem_intrp%np)
        real(rp), intent(in) :: lon(elem_intrp%np)
        real(rp), intent(in) :: lat(elem_intrp%np)
        real(rp), intent(in) :: z(elem_intrp%np)
        real(rp), intent(in) :: rplanet_
      end subroutine func
    end interface
 
    type(hexahedralelement) :: elem_intrp
    real(rp), allocatable :: x_intrp(:,:), y_intrp(:,:), z_intrp(:,:)
    real(rp), allocatable :: gam_intrp(:,:)
    real(rp), allocatable :: lon_intrp(:,:), lat_intrp(:,:)
    real(rp) :: vx(elem%nv), vy(elem%nv), vz(elem%nv)
 
    real(rp), allocatable :: intrpmat(:,:)
    real(rp), allocatable :: q_intrp(:)
 
    integer :: ke
    !-----------------------------------------------
 
    call elem_intrp%Init( intrppolyorder_h, intrppolyorder_v, .false. )
 
    allocate( intrpmat(elem%Np,elem_intrp%Np) )
    call mkinitutil_gen_gpmat( intrpmat, elem_intrp, elem )
 
    allocate( x_intrp(elem_intrp%Np,lcmesh3d%Ne), y_intrp(elem_intrp%Np,lcmesh3d%Ne), z_intrp(elem_intrp%Np,lcmesh3d%Ne) )
    allocate( gam_intrp(elem_intrp%Np,lcmesh3d%Ne) )
    allocate( lon_intrp(elem_intrp%Np,lcmesh3d%Ne), lat_intrp(elem_intrp%Np,lcmesh3d%Ne) )
    allocate( q_intrp(elem_intrp%Np) )
 
    !$omp parallel do private(vx, vy, vz)
    do ke=lcmesh3d%NeS, lcmesh3d%NeE
      vx(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),1)
      vy(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),2)
      vz(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),3)
      x_intrp(:,ke) = vx(1) + 0.5_rp * ( elem_intrp%x1(:) + 1.0_rp ) * ( vx(2) - vx(1) ) 
      y_intrp(:,ke) = vy(1) + 0.5_rp * ( elem_intrp%x2(:) + 1.0_rp ) * ( vy(4) - vy(1) )
      z_intrp(:,ke) = vz(1) + 0.5_rp * ( elem_intrp%x3(:) + 1.0_rp ) * ( vz(5) - vz(1) )
 
      gam_intrp(:,ke) = 1.0_rp
    end do
 
    call cubedspherecoordcnv_cs2lonlatpos( lcmesh3d%panelID, x_intrp, y_intrp, gam_intrp, & ! (in)
      elem_intrp%Np * lcmesh3d%Ne,                                                        & ! (in)
      lon_intrp(:,:), lat_intrp(:,:) )                                                      ! (out)
 
    !$omp parallel do private( q_intrp )
    do ke=lcmesh3d%NeS, lcmesh3d%NeE
 
      call func( q_intrp,                                & ! (out)
        lon_intrp(:,ke), lat_intrp(:,ke), z_intrp(:,ke), & ! (in)
        elem_intrp, rplanet                              ) ! (in)
      
      ! Perform Galerkin projection
      q(:,ke) = matmul( intrpmat, q_intrp )
    end do
 
    call elem_intrp%Final()
  
    return
  end subroutine mkinitutil_galerkinprojection_global
 
end module mod_mkinit_util