src/conservation_analysis.f90

説明を見る。

module conservation_analysis

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

  !* gtool **
  !*

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

  !* IGMBaseLib **
  !*

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

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

  ! 正二十面格子上の物理場データの管理
  ! Manager of the physical field data on icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, Field_IcGrid2D_Final, &
    & get_IcGrid, global_mean

  ! 正二十面体格子上の物理場に作用する微分演算
  ! Some differntial operators acting on a physical field on icosahedral grid
  !
  use Derivate_Field_IcGrid2D_Manager, only: &
    & Derivate_Field_IcGrid2D, &
    & vertical_curl_op

  !* IGModel-SW **
  !*

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

  ! 宣言文 ; Declaration statement
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedures
  !
  public :: init_conservation_analysis, calc_total_energy, calc_potential_enstrophy

  ! 非公開変数
  ! Private variables

  real(DP) :: Ep0

contains

!
subroutine init_conservation_analysis(init_h, hs)
  ! 宣言文 ; Declaration statement
  !
  type(Field_IcGrid2D), intent(in) :: init_h
  type(Field_IcGrid2D), intent(in) :: hs

  ! 作業変数
  ! Work variables
  !
  type(Field_IcGrid2D) :: tmp_Ep_field
  type(IcGrid2D_FVM), pointer :: icgrid
  integer :: rcID, i, j
  integer :: EMax, EMin

  ! 実行文 ; Executable statements
  !
  icgrid => get_IcGrid(init_h)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  call Field_IcGrid2D_Init(tmp_Ep_field, icgrid, 'tmp_PE_field', 1)

  do rcID=1, RC_REGIONS_NUM
    do j= EMin, EMax
      do i= EMin, EMax
        tmp_Ep_field%val(rcID,i,j,1) = &
          & 0.5d0 * Grav * &
          &   ( init_h%val(rcID,i,j,1)**2 - hs%val(rcID,i,j,1)**2 )

      end do
    end do
  end do

  Ep0 = global_mean(tmp_Ep_field)
  !write(*,*) 'Ep0=',Ep0
  call Field_IcGrid2D_Final(tmp_Ep_field)

end subroutine init_conservation_analysis

!
function calc_total_energy(h, hs, v) result(total_energy)
  ! 宣言文 ; Declaration statement
  !
  type(Field_IcGrid2D), intent(in) :: h
  type(Field_IcGrid2D), intent(in) :: hs
  type(Field_IcGrid2D), intent(in) :: v
  real(DP) total_energy

  ! 作業変数
  ! Work variables
  !
  type(Field_IcGrid2D) :: tmp_TE_field
  type(IcGrid2D_FVM), pointer :: icgrid
  integer :: rcID, i, j
  integer :: EMax, EMin

  ! 実行文
  !
  icgrid => get_IcGrid(h)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  call Field_IcGrid2D_Init(tmp_TE_field, icgrid, 'tmp_TE_field', 1)

  do rcID=1, RC_REGIONS_NUM
    do j= EMin, EMax
      do i= EMin, EMax
        tmp_TE_field%val(rcID,i,j,1) = &
          &   0.5d0 * ( h%val(rcID,i,j,1) - hs%val(rcID,i,j,1) ) * dot(v%val(rcID,i,j,:), v%val(rcID,i,j,:)) &
          & + 0.5d0 * Grav * ( h%val(rcID,i,j,1)**2 - hs%val(rcID,i,j,1)**2 ) &
          & - Ep0
      end do
    end do
  end do

  total_energy = global_mean(tmp_TE_field)
  call Field_IcGrid2D_Final(tmp_TE_field)

end function calc_total_energy

!
function calc_potential_enstrophy(h, hs, v, f, diff_eval) result(potential_enstrophy)
  ! 宣言文 ; Declaration statement
  !
  type(Field_IcGrid2D), intent(in) :: h
  type(Field_IcGrid2D), intent(in) :: hs
  type(Field_IcGrid2D), intent(in) :: v
  type(Field_IcGrid2D), intent(in) :: f
  type(Derivate_Field_IcGrid2D), intent(inout) :: diff_eval
  real(DP) potential_enstrophy

  ! 作業変数
  ! Work variables
  !
  type(Field_IcGrid2D) :: tmp_PE_field
  type(Field_IcGrid2D) :: vorticity
  type(IcGrid2D_FVM), pointer :: icgrid
  integer :: rcID, i, j
  integer :: EMax, EMin

  ! 実行文
  !
  icgrid => get_IcGrid(h)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  call Field_IcGrid2D_Init(tmp_PE_field, icgrid, 'tmp_PE_field', 1)
  call Field_IcGrid2D_Init(vorticity, icgrid, 'vorticity', 1)

  ! 速度場に回転演算子を作用させて, 渦度場を計算する.
  !
  call vertical_curl_op(diff_eval, v, vorticity )

  do rcID=1, RC_REGIONS_NUM
    do j= EMin, EMax
      do i= EMin, EMax
        tmp_PE_field%val(rcID,i,j,1) = &
          & 0.5d0 * ( vorticity%val(rcID,i,j,1) + f%val(rcID,i,j,1) )**2 &
          &   / ( h%val(rcID,i,j,1) - hs%val(rcID,i,j,1) )
      end do
    end do
  end do

  potential_enstrophy = global_mean(tmp_PE_field)
  call Field_IcGrid2D_Final(tmp_PE_field)

end function calc_potential_enstrophy

end module conservation_analysis