src/core/grid/STDGrid_Builder.f90

Go to the documentation of this file.

module STDGrid_Builder

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

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

  ! 数学
  ! Mathematic
  !
  use igmcore_math, only: PI     ! 円周率. PI.

  ! 線形代数
  ! Linear algebra
  !
  use igmcore_linear_algebra, only: &
    & vec_length, rotateZ

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

  ! 球面三角関数
  ! Spherical trigonometry
  !
  use igmcore_spherical_trigonometry, only: &
    & spherical_penta_area, &
    & spherical_hex_area

  ! 正二十面体格子データの管理
  ! Manager of Icosahedral grid data
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & paste_margin_width, calc_icosahedron_vertex, set_rcregion_bounds, &
    & get_IdMin, get_EffSize_Min, get_EffSize_Max, get_IcRadius, &
    & generate_CV5_GPindex, generate_CV6_GPindex, &
    & ICOSAHEDRON_VERTEX_NUM, RC_REGIONS_NUM

  ! 宣言文 ; Declaration statements
  !
  implicit none 
  private

  ! 公開手続き
  ! Public procedures
  !
  public :: construct_STDGrid
  public :: calc_STDGrid_GP_vertex, calc_CV_vertex, calc_CV_Area
  public :: project_on_sphere

contains 

!
subroutine construct_STDGrid(icgrid)

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid

  ! 実行文 ; Executable statements
  !

  call calc_STDGrid_GP_vertex(icgrid)
  call calc_CV_vertex(icgrid)
  call calc_CV_Area(icgrid)

end subroutine construct_STDGrid

!
subroutine calc_STDGrid_GP_vertex( &
  & icgrid &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid

  ! 作業変数
  ! Work variables
  !
  real(DP) :: orth_icvertex(ICOSAHEDRON_VERTEX_NUM, 3)
  integer :: ringID
  integer :: rcID

  integer :: EMin, EMax
  integer :: upper_IDs(6) = (/ 2, 3, 4, 5, 6, 2 /)
  integer :: lower_IDs(6) = (/ 7, 8, 9, 10, 11, 7 /)

  ! 実行文 ; Executable statements
  !

  !write(*,*) 'calc STDGP_vertex'
  ! (単位)正二十面体の 12 個の頂点の座標を計算する.
  call calc_icosahedron_vertex(orth_icvertex)

  ! 球の半径分だけ拡大する
  orth_icvertex(:,:) = get_IcRadius(icgrid) * orth_icvertex(:,:)

  ! 各矩形領域の端にあたる格子点が, 正二十面体の頂点に重なるように, 座標を設定する.
  do ringID=1,5
    ! 上側の矩形領域に対して
    call set_rcregion_bounds( &
      & self=icgrid, &                                    ! (inout)
      & rcID=ringID, &                                      ! (in)
      & top=orth_icvertex(1,:), &                           ! (in)
      & left=orth_icvertex(upper_IDs(ringID),:), &          ! (in)
      & bottom=orth_icvertex(lower_IDs(ringID), :), &       ! (in)
      & right=orth_icvertex(upper_IDs(ringID+1),:) )        ! (in)

    ! 下側の矩形領域に対して
    call set_rcregion_bounds( &
      & self=icgrid, &                                    ! (inout)
      & rcID=ringID+5, &                                    ! (in)
      & top=orth_icvertex(upper_IDs(ringID+1),:), &         ! (in)
      & left=orth_icvertex(lower_IDs(ringID),:), &          ! (in)
      & bottom=orth_icvertex(12,:), &                       ! (in)
      & right=orth_icvertex(lower_IDs(ringID+1),:) )        ! (in)

  end do

  ! 配列インデックスの境界に関する情報を取得
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  ! 各矩形領域の端の格子点座標が決まったので, それをもとに再帰的に各矩形領域を分割していく.
  do rcID=1, RC_REGIONS_NUM
    call split_rectangle_region( &
      & rc_AGrid=icgrid%rcs_AGrid(rcID, EMin:EMax,EMin:EMax,:), & ! (inout)
      & ic_radius=get_IcRadius(icgrid), &                         ! (in)
      & row1=EMin, row2=EMax, &                                        ! (in)
      & col1=EMin, col2=EMax )                                         ! (in)
  end do

  ! のりしろを貼り合わせる
  call paste_margin_width(icgrid)

end subroutine calc_STDGrid_GP_vertex

!
subroutine calc_CV_vertex(icgrid)

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid

  ! 作業変数
  ! Work variable
  !
  integer :: CV6_GPindx(6, 3, 2)
  integer :: CV5_GPindx(5, 3, 2)
  integer :: rcID, i, j, cvID
  integer :: EMin, EMax, idMin
  real(DP) :: ic_radius

  ! 実行文 ; Executable statemnets
  !

!  write(*,*) 'calc STDCV_vertex'

  ! 配列インデックスの境界に関する情報を取得
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)
  idMin = get_IdMin(icgrid)
  ! 半径の取得
  ic_radius = get_IcRadius(icgrid)

  ! 各格子点について, そのコントロールボリュームの頂点座標を計算する.
  do rcID=1, RC_REGIONS_NUM
    do j=EMin, EMax
      do i=EMin, EMax

        if ( ( rcID > 5 .and. i==EMin .and. j==EMin  ) .or. ( i==EMin .and. j== EMax ) &
          & .or. ( i==EMax .and. j== EMin ) .or. ( rcID <= 5 .and. i==EMax .and. j== EMax ) ) then
          ! 矩形領域の左右端の特異格子点のコントロールボリュームの頂点(5個)座標を求める.

          CV5_GPindx(:,:,:) = generate_CV5_GPindex(icgrid, GP_i=i, GP_j=j, rcID=rcID )
          do cvID=1,5
            icgrid%rcs_CV(rcID,i,j,cvID, 1:3) = &
              & calc_each_CV_vertex( icgrid%rcs_AGrid(rcID,:,:,:), CV5_GPindx(cvID,:,:), ic_radius, idMin )
          end do

        else if( rcID <= 5 .and. i==EMin .and. j==EMin  ) then
          ! 北極にある格子点が持つコントロールボリュームの頂点(5個)座標を求める.
          do cvID=1,5
             icgrid%rcs_CV(rcID,EMin,EMin, cvID, :) = &
               & (  icgrid%rcs_AGrid(cvID,EMin,EMin, :)    &
               &  + icgrid%rcs_AGrid(cvID,EMin+1,EMin, :) &
               &  + icgrid%rcs_AGrid(cvID,EMin,EMin+1, :) &
               &  ) / 3.0d0

             call project_on_sphere( ic_radius, icgrid%rcs_CV(rcID,i,j,cvID,:) )
           end do

        else if( rcID > 5 .and. i==EMax .and. j== EMax ) then
          ! 南極にある格子点が持つコントロールボリュームの頂点(5個)座標を求める.
          do cvID=1,5
             icgrid%rcs_CV(rcID,EMax,EMax, cvID, :) = &
               & (   icgrid%rcs_AGrid(11-cvID,EMax,EMax, :) &
               &   + icgrid%rcs_AGrid(11-cvID,EMax-1,EMax, :) &
               &   + icgrid%rcs_AGrid(11-cvID,EMax,EMax-1, :) &
               &  ) / 3.0d0

             call project_on_sphere(ic_radius, icgrid%rcs_CV(rcID,i,j,cvID,:))
           end do

        else
          ! 特異点以外の格子点のコントロールボリュームの頂点(6個)座標を求める.
          CV6_GPindx(:,:,:) =  generate_CV6_GPindex(icgrid, GP_i=i, GP_j=j, rcID=rcID )

          do cvID=1,6
            icgrid%rcs_CV(rcID,i,j,cvID, 1:3) = &
              & calc_each_CV_vertex(icgrid%rcs_AGrid(rcID,:,:,:), CV6_GPindx(cvID,:,:), ic_radius, idMin)
          end do

        end if
      end do
    end do
  end do

end subroutine calc_CV_vertex

!
subroutine calc_CV_Area( &
  & icgrid &
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid


  ! 作業変数
  ! Work variables
  !
  integer :: rcID, i, j
  integer :: EMin, EMax, idMin, k
  real(DP) :: ic_radius
  real(DP) :: CV5(5, 3), CV6(6, 3)

  ! 実行文 ; Executable statements
  !

!  write(*,*) 'calc CV_Area'

  ! 配列インデックスの境界に関する情報を取得
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)
  ! 半径の取得
  ic_radius = get_IcRadius(icgrid)

  ! それぞれの格子点について, それに付随するコントールボリュームの面積を求める.
  do rcID=1,RC_REGIONS_NUM
    do j=EMin, EMax
      do i=EMin, EMax
        if ( ( i==EMin .and. j==EMin ) .or. ( i==EMax .and. j== EMin ) .or. &
          &  ( i==EMin .and. j== EMax ) .or. ( i==EMax .and. j== EMax ) ) then
          ! 特異格子点のコントールボリューム(頂点数 5 個)に対して

          CV5(:,:) = icgrid%rcs_CV(rcID,i,j,1:5, :)

          icgrid%rcs_CVArea(rcID,i,j) = &
            & spherical_penta_area( ic_radius, CV5(1,:), CV5(2,:), CV5(3,:),CV5(4,:),CV5(5,:) )

        else
          ! 通常の格子点のコントールボリューム(頂点数 6 個)に対して

          CV6(:,:) = icgrid%rcs_CV(rcID,i,j,1:6, :)

          icgrid%rcs_CVArea(rcID,i,j) = &
            & spherical_hex_area( ic_radius, CV6(1,:), CV6(2,:), CV6(3,:),CV6(4,:),CV6(5,:), CV6(6,:) )

        end if
      end do
    end do
  end do

end subroutine calc_CV_Area


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! モジュール外非公開手続き
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
recursive subroutine split_rectangle_region( &
  & rc_AGrid,                         &  ! (inout)
  & ic_radius, row1, row2, col1, col2 &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(inout) :: rc_AGrid(:,:,:)
  real(DP), intent(in) :: ic_radius
  integer, intent(in) :: row1
  integer, intent(in) :: row2
  integer, intent(in) :: col1
  integer, intent(in) :: col2

  ! 作業変数
  ! Work variable
  !
  integer mid_row, mid_col

  ! 実行文 ; Executable statements
  !

  if( mod(row1+row2, 2) == 1 ) then
    ! 分割予定の格子点が存在しないので, 再帰処理を終了する.
    return
  end if

  ! 対象の矩形領域の各辺の中点と矩形領域内の中心点の座標を計算する.
  mid_row = (row1 + row2) / 2
  mid_col = (col1 + col2) / 2

  rc_AGrid(row1, mid_col, :) = 0.5d0 * ( rc_AGrid(row1, col1, :) + rc_AGrid(row1, col2, :) )
  rc_AGrid(row2, mid_col, :) = 0.5d0 * ( rc_AGrid(row2, col1, :) + rc_AGrid(row2, col2, :) )
  rc_AGrid(mid_row, col1, :) = 0.5d0 * ( rc_AGrid(row1, col1, :) + rc_AGrid(row2, col1, :) )
  rc_AGrid(mid_row, col2, :) = 0.5d0 * ( rc_AGrid(row1, col2, :) + rc_AGrid(row2, col2, :) )
  rc_AGrid(mid_row, mid_col, :) = 0.5d0 * ( rc_AGrid(row2, col1, :) + rc_AGrid(row1, col2, :) )

  ! 平面図形として矩形平面上に存在する 5 点を, 原点を中心として球面上に射影する.
  call project_on_sphere(ic_radius, rc_AGrid(row1, mid_col, :))
  call project_on_sphere(ic_radius, rc_AGrid(row2, mid_col, :))
  call project_on_sphere(ic_radius, rc_AGrid(mid_row, col1, :))
  call project_on_sphere(ic_radius, rc_AGrid(mid_row, col2, :))
  call project_on_sphere(ic_radius, rc_AGrid(mid_row, mid_col, :))

  ! 頂点の座標が決まった４つの矩形を, それぞれさらに分割する.
  call split_rectangle_region(rc_AGrid, ic_radius, row1, mid_row, col1, mid_col)
  call split_rectangle_region(rc_AGrid, ic_radius, row1, mid_row, mid_col, col2)
  call split_rectangle_region(rc_AGrid, ic_radius, mid_row, row2, col1, mid_col)
  call split_rectangle_region(rc_AGrid, ic_radius, mid_row, row2, mid_col, col2)

end subroutine split_rectangle_region

!
subroutine project_on_sphere( &
  & radius, &  ! (in)
  & vec     &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  real(8), intent(in) :: radius    ! 正二十面体格子を内包する球の半径
  real(8), intent(inout) :: vec(3) ! 射影したい点の位置ベクトル

  ! 実行文 ; Executable statements
  !

  vec = radius / vec_length(vec) * vec
end subroutine project_on_sphere

!
function calc_each_CV_vertex( &
  & rc_AGrid, GP3_index, ic_radius, idMin &  ! (in)
  & ) result(cv_pt)

  ! 宣言文 ; Declaration statements
  !
  integer, intent(in) :: idMin
  real(DP), intent(in) :: rc_AGrid(idMin:,idMin:,:)
  integer, intent(in) :: GP3_index(:,:)
  real(DP), intent(in) :: ic_radius
  real(DP) cv_pt(3)

  ! 作業変数
  ! Work variables
  !
  integer  k

  ! 実行文 ; Executable statements
  !

  cv_pt(:) = 0.0d0

  do k=1,3
    cv_pt(:) = cv_pt(:) + rc_AGrid(GP3_index(k, 1), GP3_index(k, 2), :)
  end do

  cv_pt(:) = cv_pt(:) / 3.0d0
  call project_on_sphere(ic_radius, cv_pt(:))

end function calc_each_CV_vertex

end module STDGrid_Builder