!> \~japanese
!! \brief 剛体回転場のような典型的な物理場を生成する手続きを提供するモジュール.
!!
!! \par Field_Pattern_Builder モジュールについて
!!   - 本モジュールの手続きによって, 剛体回転場・惑星渦度場・コサインベル型の物理場などといった,
!!     地球流体力学における数値シミュレーションでよく用いられる物理場を生成することができる.
!!   - 生成された物理場データは, 正二十面格子上の物理場を管理する構造体 Field_IcGrid2D に格納される.
!!
!! \~english
!! \brief This module provides some subroutines to create a typical physical field, such as solid body rotation.
!!
!! \par About Field_Pattern_Builder module
!! - By the subroutines in this module, one can create the physical fields commonly used for numerical simulation in geophysical fluid dynamics,
!!   such as solid body rotation, planetary vorticity or cosin bell field.
!! - The created physical field data will be stored in derived type Field_IcGrid2D which manages a physical field on icosahedral grid.
!!
!! \~
!! \since IGMBaseLib 1.0
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module Field_Pattern_Builder

  ! モジュール引用 ; Use statements
  !

  !* gtool ********************
  !*

  ! 種類型パラメータ
  ! Kind type parameter 
  !
  use dc_types, only: DP  ! 倍精度実数型. Double precision. 
  
  !* IGMBaseLib ********************
  !*

  ! 数学
  ! Mathematic
  !
  use igmcore_math, only: PI

  ! 線形代数
  ! Linear algebra
  !
  use igmcore_linear_algebra, only: &
    & cross, rotateY

  ! 座標変換
  ! Coordinate conversion
  !
  use igmcore_coordinate_conversion, only: &
    & orth_to_geo_pos, geo_to_orth_pos

  ! 球面三角法
  ! Spherical trigometry
  !
  use igmcore_spherical_trigonometry, only: &
    & geodesic_arc_length
  
  ! 正二十面体格子の管理
  ! Manager of icosahedral grid
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & check_pole, &
    & get_EffSize_Min, get_EffSize_Max, get_IcRadius, &
    & get_LocalProcRectInfo, &
    & RC_REGIONS_NUM, NORTH_POLE_FLAG, SOUTH_POLE_FLAG, NOT_POLE_FLAG

  ! 正二十面格子における広域矩形領域のサブドメインである局所矩形領域(各 MPI プロセスと対応)の情報管理
  ! Manager of information about a local rectangle region (corresponding to each MPI process)
  !
  use LocalProcRect_Manager, only: &
    & LocalProcRectInfo, &
    & get_InteriorRegion_IndexInfo

  ! 正二十面体格子上の物理場データの管理
  ! Manager of a physical field data on icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & get_icgrid

  ! 宣言文 ; Declaration statement
  ! 
  implicit none
  private

  ! 公開手続き
  ! Public procedures
  !
  public :: create_solid_rotation_field, create_planetaryVorticity_field
  public :: create_cosine_bell_field

contains

!> \~japanese
!! \brief 剛体回転場を生成し, その物理場データを指定された構造型 Field_IcGrid2D の変数に保存する.
!!
!! @param[in,out] field         生成した物理場を保存する構造型 Field_IcGrid2D の変数.
!! @param[in]     angular_speed 剛体回転軸の角速度 [rad/s].
!! @param[in]     alpha         xz 平面において, デカルト座標系の z 軸と剛体回転軸がなす角度 [rad].
!!
!! \~english
!! \brief Generates the solid body rotation, and store the physical field data into the specified variable of derived type Field_IcGrid2D.
!!
!! @param[in,out] field         The variable of derived type Field_IcGrid2D in which the created physical data will be stored.
!! @param[in]     angular_speed Angular speed [rad/s].
!! @param[in]     alpha         The angle between z axis of the cartesian coordinate system and the axis of rotation lying in the x-z plane [rad].
!!
subroutine create_solid_rotation_field(field, angular_speed, alpha)

  ! 宣言文 ; Declaration statement
  !
  type(Field_IcGrid2D), intent(inout) :: field
  real(DP), intent(in) :: angular_speed
  real(DP), intent(in) :: alpha

  ! 作業変数
  ! Work variables
  !
  type(IcGrid2D_FVM), pointer :: icgrid
  type(LocalProcRectInfo), pointer :: lcRCInfo

  integer :: rcID, i,j
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  real(DP) :: Omega_vec(3)

  ! 実行文 ; Executable statement
  !

  !
  icgrid => get_icgrid(field)
  lcRCInfo => get_LocalProcRectInfo(icgrid)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)

  Omega_vec(:) = (/ sin(alpha), 0.0d0, cos(alpha) /)

  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jEMax
      do i=iEMin, iEMax

        !
        field%val(:,i,j,rcID) = &
          & angular_speed * cross( Omega_vec(:), icgrid%rcs_AGrid(:,i,j,rcID))
        end do
    end do
  end do

end subroutine create_solid_rotation_field

!> \~japanese
!! \brief 惑星渦度場を生成し, その物理場データを指定された構造型 Field_IcGrid2D の変数に保存する.
!!
!! @param[in,out] field         生成した物理場を保存する構造型 Field_IcGrid2D の変数.
!! @param[in]     angular_speed 自転角速度 [rad/s].
!! @param[in]     alpha         xz 平面において, デカルト座標系の z 軸と剛体回転軸がなす角度 [rad].
!!
!! \~english
!! \brief Generates the planetary vorticity, and store the physical field data into the specified variable of derived type Field_IcGrid2D.
!!
!! @param[in,out] field         The variable of derived type Field_IcGrid2D in which the created physical data will be stored.
!! @param[in]     angular_speed The angular velocity [rad/s].
!! @param[in]     alpha         The angle between z axis of the cartesian coordinate system and the axis of rotation lying in the x-z plane [rad].
!!
subroutine create_planetaryVorticity_field(field, angular_speed, alpha)

  ! 宣言文 ; Declaration statement
  !
  type(Field_IcGrid2D), intent(inout) :: field
  real(DP), intent(in) :: angular_speed
  real(DP), intent(in) :: alpha

  ! 作業変数
  ! Work variables
  !
  type(IcGrid2D_FVM), pointer :: icgrid
  type(LocalProcRectInfo), pointer :: lcRCInfo
  integer :: rcID, i,j
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  real(DP) :: s_pos(3)

  ! 実行文 ; Executable statement
  !

  !
  icgrid => get_icgrid(field)
  lcRCInfo => get_LocalProcRectInfo(icgrid)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)

  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jEMax
      do i=iEMin, iEMax
         s_pos(:) = orth_to_geo_pos(icgrid%rcs_Agrid(:,i,j,rcID))

        if ( check_pole(icgrid, rcID,i,j) == NOT_POLE_FLAG ) then
          field%val(:,i,j,rcID) = 2.0d0 * angular_speed * ( &
          &    cos(s_pos(1)) * cos(s_pos(2)) * sin(alpha) + sin(s_pos(2)) * cos(alpha) &
            &   )
        else
          field%val(:,i,j,rcID) = 2.0d0 * angular_speed * sin(s_pos(2)) * cos(alpha)
        end if
      end do
    end do
  end do

end subroutine create_planetaryVorticity_field

!
!> \~japanese
!! \brief コサインベル型の物理場を生成し, その物理場データを指定された構造型 Field_IcGrid2D の変数に保存する.
!!
!! @param[in,out] field         生成した物理場を保存する構造型 Field_IcGrid2D の変数.
!! @param[in]     lambda_c コサインベルの中心点の経度.
!! @param[in]     theta_c  コサインベルの中心点の緯度.
!! @param[in]     cbR      コサインベルの半径.
!! @param[in]     h0       コサインベルの最大の高さ.
!!
!! \~english
!! \brief Generates the  cosine bell field, and store the physical field data into the specified variable of derived type Field_IcGrid2D.
!!
!! @param[in,out] field    The variable of derived type Field_IcGrid2D in which the created physical data will be stored.
!! @param[in]     lambda_c The longitude at which the center of the cosine bell is located.
!! @param[in]     theta_c  The latitude at which the center of the cosine bell is located.
!! @param[in]     cbR      The radius of cosine bell.
!! @param[in]     h0       The maximum height of the cosine bell.
!!
subroutine create_cosine_bell_field(field, lambda_c, theta_c, cbR, h0)

  ! 宣言文 ; Declaration statement
  !
  type(Field_IcGrid2D), intent(inout) :: field
  real(DP), intent(in) :: lambda_c
  real(DP), intent(in) :: theta_c
  real(DP), intent(in) :: cbR
  real(DP), intent(in) :: h0
  
  ! 作業変数
  ! Work variables
  !
  type(IcGrid2D_FVM), pointer :: icgrid
  type(LocalProcRectInfo), pointer :: lcRCInfo
  real(DP) :: ic_radius
  real(DP) :: orth_c(3)

  !> \f$( \lambda, \theta )\f$ とコサインベルの中心の間の大円距離. 
  !! The great circle distance between \f$( \lambda, \theta )\f$ and the center of cosin bell.  
  real(DP) :: r

  integer :: rcID, i,j
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax

  ! 実行文 ; Executable statement
  !
  
  !
  icgrid => get_icgrid(field)
  ic_radius = get_IcRadius(icgrid)
  lcRCInfo => get_LocalProcRectInfo(icgrid)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)
  
  orth_c = geo_to_orth_pos( (/lambda_c, theta_c, ic_radius /) )
 
  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jEMax
      do i=iEMin, iEMax

         ! 球面上のある一点とコサインベル中心点との大円距離を求める.
        r = geodesic_arc_length(ic_radius, orth_c, icgrid%rcs_AGrid(:,i,j,rcID))
        
        if( r < cbR ) then
          field%val(:,i,j,rcID) = 0.5d0 * h0 * ( 1.0d0 + cos(PI*r/cbR) )
        else
          field%val(:,i,j,rcID) = 0.0d0
        end if

      end do
    end do
  end do

end subroutine create_cosine_bell_field

end module Field_Pattern_Builder
