src/util/field/Field_Pattern_Builder.f90

Go to the documentation of this file.

module Field_Pattern_Builder

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

  ! 種類型パラメータ
  ! Kind type parameter 
  !
  use dc_types, only: DP  ! 倍精度実数型. Double precision. 
  
  ! 数学
  ! 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, &
    & RC_REGIONS_NUM, NORTH_POLE_FLAG, SOUTH_POLE_FLAG, NOT_POLE_FLAG

  ! 正二十面体格子上の物理場データの管理
  ! 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

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
  integer :: rcID, i,j
  integer :: EMin, EMax
  real(DP) :: Omega_vec(3)

  ! 実行文 ; Executable statement
  !

  !
  icgrid => get_icgrid(field)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)
  Omega_vec(:) = (/ sin(alpha), 0.0d0, cos(alpha) /)

  do rcID=1, RC_REGIONS_NUM
    do j=EMin, EMax
      do i=EMin, EMax

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

end subroutine create_solid_rotation_field

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
  integer :: rcID, i,j
  integer :: EMin, EMax
  real(DP) :: s_pos(3)

  ! 実行文 ; Executable statement
  !

  !
  icgrid => get_icgrid(field)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  do rcID=1, RC_REGIONS_NUM
    do j=EMin, EMax
      do i=EMin, EMax
         s_pos(:) = orth_to_geo_pos(icgrid%rcs_Agrid(rcID,i,j,:))

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

end subroutine create_planetaryVorticity_field

!
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
  real(DP) :: ic_radius
  real(DP) :: orth_c(3)

  real(DP) :: r

  integer :: rcID, i,j
  integer :: EMin, EMax

  ! 実行文 ; Executable statement
  !
  
  !
  icgrid => get_icgrid(field)
  ic_radius = get_IcRadius(icgrid)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)
  
  orth_c = geo_to_orth_pos( (/lambda_c, theta_c, ic_radius /) )
 
  do rcID=1, RC_REGIONS_NUM
    do j=EMin, EMax
      do i=EMin, EMax
        r = geodesic_arc_length(ic_radius, orth_c, icgrid%rcs_AGrid(rcID,i,j,:))
        
        if( r < cbR ) then
          field%val(rcID,i,j, :) = 0.5d0 * h0 * ( 1.0d0 + cos(PI*r/cbR) )
        else
          field%val(rcID,i,j, :) = 0.0d0
        end if

      end do
    end do
  end do

end subroutine create_cosine_bell_field

end module Field_Pattern_Builder