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_calc_cosinebell
  public :: mkinitutil_calc_cosinebell_global
  public :: mkinitutil_galerkinprojection
  public :: mkinitutil_galerkinprojection_global
 
contains
  !> 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
    integer :: i, j
 
    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 :: l2projmat(:,:)
    real(rp), allocatable :: q_intrp(:)
    real(rp) :: s
 
    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( l2projmat(elem%Np,elem_intrp%Np) )
    call elem%Generate_L2ProjMat( elem_intrp, & ! (in)
      l2projmat ) ! (out)
 
    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( q_intrp, vx, vy, vz, r_intrp, s )
    !$acc data create( x_intrp, y_intrp, z_intrp, z_func) copyin(L2ProjMat)
 
    !$omp do
    !$acc parallel loop private(vx, vy, vz) present(lcmesh3D, elem_intrp, x_intrp, y_intrp, z_intrp, z_func)
    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
        !$acc parallel loop present(z_func, z_intrp) copyin(z_func_params)
        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
    !$acc parallel loop private( r_intrp, q_intrp ) present(x_intrp, y_intrp, z_intrp, z_func, q)
    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
 
      ! Perform Galerkin projection
      !$acc loop worker
      do i=1, elem%Np
        s = 0.0_rp
        !$acc loop vector reduction(+:s)
        do j=1, elem_intrp%Np
          s = s + l2projmat(i,j) * q_intrp(j) * z_func(j,ke)
        end do
        q(i,ke) = s
      end do
    end do
    !$omp end do
 
    !$acc end data
    !$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
    integer :: i, j
 
    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 :: l2projmat(:,:)
    real(rp), allocatable :: q_intrp(:)
    real(rp) :: s
 
    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( l2projmat(elem%Np,elem_intrp%Np) )
    call elem%Generate_L2ProjMat( elem_intrp, & ! (in)
      l2projmat ) ! (out)
 
    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 private( q_intrp, vx, vy, vz, r_intrp, s )
    !$acc data create( x_intrp, y_intrp, z_intrp, lon_intrp, lat_intrp, gam_intrp, z_func) copyin(L2ProjMat)
 
    !$omp do
    !$acc parallel loop private(vx, vy, vz) present(lcmesh3D, elem_intrp, x_intrp, y_intrp, z_intrp, z_func, gam_intrp)
    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
    !$omp 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 do
        !$acc parallel loop present(z_func, z_intrp) copyin(z_func_params)
        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 do
    !$acc parallel loop private(r_intrp, q_intrp) present(x_intrp, y_intrp, z_intrp, z_func, lon_intrp, lat_intrp, q)
    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
      !$acc loop worker
      do i=1, elem%Np
        s = 0.0_rp
        !$acc loop vector reduction(+:s)
        do j=1, elem_intrp%Np
          s = s + l2projmat(i,j) * q_intrp(j) * z_func(j,ke)
        end do
        q(i,ke) = s
      end do
    end do
    !$omp end do
 
    !$acc end data
    !$omp end parallel
 
    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, lcmesh3D, elem_intrp   )
      import localmesh3d
      import elementbase3d
      import rp
      class(localmesh3d), intent(in) :: lcmesh3d
      class(elementbase3d), intent(in) :: elem_intrp
      real(rp), intent(out) :: q_intrp(elem_intrp%np,lcmesh3d%ne)
      real(rp), intent(in) :: x(elem_intrp%np,lcmesh3d%ne)
      real(rp), intent(in) :: y(elem_intrp%np,lcmesh3d%ne)
      real(rp), intent(in) :: z(elem_intrp%np,lcmesh3d%ne)
    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 :: l2projmat(:,:)
  real(rp), allocatable :: q_intrp(:,:)
 
  integer :: ke, p
  integer :: i, j
  real(rp) :: s
  !-----------------------------------------------
 
  call elem_intrp%Init( intrppolyorder_h, intrppolyorder_v, .false. )
 
  allocate( l2projmat(elem%Np,elem_intrp%Np) )
  call elem%Generate_L2ProjMat( elem_intrp, & ! (in)
    l2projmat ) ! (out)
 
  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,lcmesh3d%Ne) )
 
  !$acc data create( x_intrp, y_intrp, z_intrp, q_intrp ) copyin(L2ProjMat)
 
  !$omp parallel do private(vx, vy, vz)
  !$acc parallel loop gang private(vx, vy, vz) present(x_intrp, y_intrp, z_intrp, lcmesh3D,elem_intrp)
  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)
    !$acc loop vector
    do p=1, elem_intrp%Np
      x_intrp(p,ke) = vx(1) + 0.5_rp * ( elem_intrp%x1(p) + 1.0_rp ) * ( vx(2) - vx(1) ) 
      y_intrp(p,ke) = vy(1) + 0.5_rp * ( elem_intrp%x2(p) + 1.0_rp ) * ( vy(4) - vy(1) )
      z_intrp(p,ke) = vz(1) + 0.5_rp * ( elem_intrp%x3(p) + 1.0_rp ) * ( vz(5) - vz(1) )
    end do
  end do
 
  call func( q_intrp,              & ! (out)
    x_intrp, y_intrp, z_intrp,     & ! (in)
    lcmesh3d, elem_intrp           ) ! (in)
 
  !$omp parallel do private( s )
  !$acc parallel loop gang private(s) present(q_intrp, q, L2ProjMat, lcmesh3D,elem_intrp,elem)
  do ke=lcmesh3d%NeS, lcmesh3d%NeE
    ! Perform Galerkin projection
    !$acc loop worker
    do i=1, elem%Np
      s = 0.0_rp
      !$acc loop vector reduction(+:s)
      do j=1, elem_intrp%Np
        s = s + l2projmat(i,j) * q_intrp(j,ke)
      end do
      q(i,ke) = s
    end do
  end do
 
  !$acc end data
 
  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, lcmesh3D, elem_intrp, rplanet   )
        import localmesh3d
        import elementbase3d
        import rp
        class(localmesh3d), intent(in) :: lcmesh3d
        class(elementbase3d), intent(in) :: elem_intrp
        real(rp), intent(out) :: q_intrp(elem_intrp%np,lcmesh3d%ne)
        real(rp), intent(in) :: lon(elem_intrp%np,lcmesh3d%ne)
        real(rp), intent(in) :: lat(elem_intrp%np,lcmesh3d%ne)
        real(rp), intent(in) :: z(elem_intrp%np,lcmesh3d%ne)
        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 :: l2projmat(:,:)
    real(rp), allocatable :: q_intrp(:,:)
 
    integer :: ke
    integer :: i, j
 
    real(rp) :: s
    !-----------------------------------------------
 
    call elem_intrp%Init( intrppolyorder_h, intrppolyorder_v, .false. )
 
    allocate( l2projmat(elem%Np,elem_intrp%Np) )
    call elem%Generate_L2ProjMat( elem_intrp, & ! (in)
      l2projmat ) ! (out)
 
    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,lcmesh3d%Ne) )
 
    !$acc data create( x_intrp, y_intrp, z_intrp, gam_intrp, lon_intrp, lat_intrp, q_intrp ) copyin(L2ProjMat)
 
    !$omp parallel do private(vx, vy, vz)
    !$acc parallel loop private(vx, vy, vz) present(x_intrp, y_intrp, z_intrp, gam_intrp, L2ProjMat)
    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)
 
    call func( q_intrp,              & ! (out)
      lon_intrp, lat_intrp, z_intrp, & ! (in)
      lcmesh3d, elem_intrp, rplanet  ) ! (in)
 
    !$omp parallel do private( s )
    !$acc parallel loop gang private(s) present(q_intrp, q, L2ProjMat, lcmesh3D,elem_intrp,elem)
    do ke=lcmesh3d%NeS, lcmesh3d%NeE
      ! Perform Galerkin projection
      !$acc loop worker
      do i=1, elem%Np
        s = 0.0_rp
        !$acc loop vector reduction(+:s)
        do j=1, elem_intrp%Np
          s = s + l2projmat(i,j) * q_intrp(j,ke)
        end do
        q(i,ke) = s
      end do
    end do
 
    !$acc end data
 
    call elem_intrp%Final()
    return
  end subroutine mkinitutil_galerkinprojection_global
 
end module mod_mkinit_util