src/motion_equation.f90

説明を見る。

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

  !* gtool ****
  !*

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

  !* IGMBaseLib ****
  !*

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

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

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


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

  ! 正二十面体格子の管理
  ! Manager of icoshedral grid
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_EffSize_Min, get_EffSize_Max, check_pole, &
    & RC_REGIONS_NUM, NOT_POLE_FLAG

  ! 正二十面体格子上の物理場に作用する微分演算を提供するクラス
  ! Class providing some differential operators acting on a physcial field on icosahedral grid.  
  !
  use Derivate_Field_IcGrid2D_Manager, only: &
    & Derivate_Field_IcGrid2D, &
    & vertical_curl_op, gradient_op

  !* IGModel-SW ****
  !*

  ! シミュレーションパラメータの管理
  ! Manager of simulation parameter
  !
  use param_manager, only: &
    & Grav, Omega, alpha

  ! OpenMP
  !
  !
  use omp_lib

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public statements
  !
  public :: motion_equation_Init, calc_motion_eq_dvdt

  ! 非公開変数
  ! Private statements
  !

  type(Field_IcGrid2D), save :: xy_Zeta

  type(Field_IcGrid2D), save :: mecha_energy

  type(Field_IcGrid2D), save :: grad_mecha_energy

contains

!
subroutine motion_equation_Init( &
    & icgrid                      &  ! (inout)
    & )
  
  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid
 
  ! 実行文 ; Executable statement 
  !

  ! 渦度場, 運動エネルギーの場, 運動エネルギーの勾配場を管理する, 構造型 Field_IcGrid2D の変数を初期化する.
  !

  call Field_IcGrid2D_Init( &
    & self=xy_Zeta, &
    & icgrid=icgrid, name='relative vorticity field', rank=1 )

  call Field_IcGrid2D_Init( &
    & self=mecha_energy, &
    & icgrid=icgrid, name='mechanical_energy_field', rank=1 )

  call Field_IcGrid2D_Init( &
    & self=grad_mecha_energy, &
    & icgrid=icgrid, name='gradient of mechanical energy field', rank=3 )
  
end subroutine motion_equation_Init

!
subroutine calc_motion_eq_dvdt( &
  & DVelDtN,                       &  ! (inout)
  & xy_VelN, xy_HeightN, xy_Coli,  &  ! (in)
  & diff_eval                      &  ! (inout)
  & )
 
  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid2D), intent(inout) :: DVelDtN
  type(Field_IcGrid2D), intent(in) :: xy_VelN
  type(Field_IcGrid2D), intent(in) :: xy_HeightN
  type(Field_IcGrid2D), intent(in) :: xy_Coli
  type(Derivate_Field_IcGrid2D), intent(inout) :: diff_eval


  ! 作業変数
  ! Work variables
  !
  integer :: EMax, EMin
  integer :: rcID,i,j
  real(DP) :: s_pos(3)
  real(DP) :: unit_k(3)
  type(IcGrid2D_FVM), pointer :: icgrid

  ! 実行文 ; Executable statements
  !

  icgrid => get_icgrid(xy_VelN)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  
  ! 速度場に回転演算子を作用させて, 相対渦度場を計算する. 
  ! Curl operator acts on the velocity field to calculate the relative vorticity field.  
  call vertical_curl_op(diff_eval, xy_VelN, xy_Zeta )
  
  ! ジオポテンシャルと運動エネルギーの和の空間勾配場を計算する.
  ! Calculates the spatial gradient geopotential and kinetic energy.
  !

  ! まず, 各格子点において力学的エネルギーを計算する.
  ! First, calculates the mechanical energy at each grid point.
  do rcID=1, RC_REGIONS_NUM
    !$omp parallel do
    do j=EMin, EMax
      do i=EMin, EMax
        mecha_energy%val(rcID,i,j,:) = &
          &   Grav * xy_HeightN%val(rcID,i,j,:) &
          & + 0.5d0 * dot( xy_VelN%val(rcID,i,j,:), xy_VelN%val(rcID,i,j,:) )
      end do
    end do
  end do

  ! 袖領域の貼りあわせ
  call paste_field2D_margin(mecha_energy)

  ! 力学的エネルギーの場に, 勾配演算子を作用させる. 
  ! Gradient operator acts on the mechanical energy field. 
  call gradient_op(diff_eval, mecha_energy, grad_mecha_energy)
  

  ! 運動方程式の半離散式の右辺を評価し, タイムレベル n における速度変化率を計算する. 
  ! Evaluates the right-hand side of the semi-discreted equation of motion,
  ! and calculates the rate of change of velocity at time level n with respect to time.
  !
  
  do rcID=1, RC_REGIONS_NUM
    !$omp parallel do private(s_pos, unit_k)
    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
          unit_k(:) = local_k(s_pos(:))
        else
          ! 両極
          if ( rcID <= 5 ) then
            unit_k(:) = (/ 0.0d0, 0.0d0, 1.0d0 /)
            s_pos(2) = PI / 2.0d0
          else
            unit_k(:) = (/ 0.0d0, 0.0d0, -1.0d0 /)
            s_pos(2) = - PI / 2.0d0
          end if
        end if

       ! 速度変化率を計算する.
       DVelDtN%val(rcID,i,j,:) = &
         & - ( xy_Zeta%val(rcID,i,j,1) + xy_Coli%val(rcID,i,j,1) ) * cross( unit_k(:), xy_VelN%val(rcID,i,j,:) ) &
         & - grad_mecha_energy%val(rcID,i,j,:)
        
      end do
    end do
  end do
  
end subroutine calc_motion_eq_dvdt

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 非公開手続き
! Private procedure
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
function local_k( &
    & s_pos        &  ! (in)
    & ) result(vec)

  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(in) :: s_pos(3)
  real(DP) :: vec(3)
  
  ! 作業変数
  ! Work variables
  !
  
  ! 地理座標系における鉛直方向の単位ベクトル.
  real(DP), parameter :: s_k(3) = (/ 0.0d0, 0.0d0, 1.0d0 /) 

  ! 実行文
  !

  ! 地理座標系における鉛直方向の単位ベクトルを, 直交座標系へと変換する. 
  vec = geo_to_orth_vec( s_k, s_pos)

end function local_k

end module motion_equation