src/Williamson1992_TestCase/class_TestCase2.f90

Go to the documentation of this file.

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

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

  !
  !
  !
  use gtool_history

  ! * IGMBaseLib *
  ! *

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

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

  ! 正二十面体格子の管理クラス
  !
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_glevel, get_IcRadius, get_EffSize_Min, get_EffSize_Max, check_pole, &
    & RC_REGIONS_NUM, NOT_POLE_FLAG
  
  ! 正二十面体格子上に分布する物理場の管理クラス
  !
  !
  use Field_IcGrid2D_manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, get_icgrid
  
  !
  !
  !
  use Field_Pattern_Builder, only: &
    & create_solid_rotation_field

  ! 正二十面体格子上の物理場における誤差ノルムチェック
  !
  !
  use Field_IcGrid2D_error_norm_check, only: &
    & numerical_error_norm1, numerical_error_norm2, &
    & numerical_error_norminfinity 


  !
  !
  !
  use IcGrid_ncWriter_mod, only: &
    & IcGrid_ncWriter, IcGrid_ncWriter_Init, &
    & open_ncFile, ncdef_simulation_parameter, ncdef_GridData, ncdef_FieldData, end_ncdef, &
    & write_GridData, write_FieldData, close_ncFile, &
    & increase_recorde_counter


  ! * IGModel-SW ****
  ! *


  ! パラメータ管理
  !
  !
  use param_manager, only: &
    & alpha, output_tick, Grav, Omega, &
    & integration_time, DelTime


  ! 宣言文 ; Declaration statement
  !
  implicit none
  private
  
  ! 公開変数
  ! Public variables
  !
  
  real(DP) :: h_0

  real(DP) :: u_0

  real(DP) :: angular_speed

  ! クラス定義
  ! Class definition
  !

  !
  type, public :: TestCase2

    ! 宣言文 ; Declaration statement
    !

    ! 非公開メンバ変数
    ! Private member variable
    !

    type(Field_IcGrid2D) :: true_h
    ! f2003: type(Field_IcGrid2D), private :: true_h

    
    type(Field_IcGrid2D) :: error_h
    ! f2003: type(Field_IcGrid2D), private :: error_h
    
    type(Field_IcGrid2D) :: true_v
    ! f2003: type(Field_IcGrid2D), private :: true_v

    
    integer :: errorh_ncVarID
    ! f2003: integer, private :: errorh_ncVarID

  end type TestCase2

  !
  !
  !
  interface initialize_TestCase
    module procedure init_TestCase2
  end interface initialize_TestCase

  interface finalize_TestCase
    module procedure finalize_TestCase2
  end interface finalize_TestCase

  public :: initialize_TestCase, finalize_TestCase
  public :: set_initial_h, set_initial_hs, set_initial_v, timelevel_Updated


  ! 非公開変数
  ! Private variables
  !
  
  character(TOKEN) :: filename = 'error_norm.dat'
  
  type(IcGrid_ncWriter), save :: writer

contains

!
subroutine init_TestCase2(self, icgrid_ref )
  ! 宣言文 ; Declaration statement
  !
  type(TestCase2), intent(inout) :: self
  type(IcGrid2D_FVM), intent(in) :: icgrid_ref

  ! 作業変数
  ! Work variables
  !
  real(DP) :: ic_radius

  ! 実行文 ; Executable statement
  !

  ic_radius = get_IcRadius(icgrid_ref)
  
 
  u_0 = 2.0d0 * PI * ic_radius / ( 12.0d0 * 24.0d0 * 3600d0 )
  h_0 = 2.94e04 / Grav
  angular_speed = 2.0d0 * PI / ( 12.0d0 * 24.0d0 * 3600d0 )

  !
  call Field_IcGrid2D_Init(self%true_h, icgrid_ref, 'true_h', 1)
  call Field_IcGrid2D_Init(self%error_h, icgrid_ref, 'error_h', 1)

  !
  call Field_IcGrid2D_Init(self%true_v, icgrid_ref, 'true_v', 3)

  ! 数値誤差ノルムを書きだすための, NetCDF ファイルを用意する.  
  
  call HistoryCreate( &
    & file='error_norm.nc', title='error norm check', &
    & source=' ', institution=' ', &
    & dims=(/ 't' /), dimsizes=(/ 0 /), &
    & longnames=(/ 'time' /), units=(/ 's' /), &
    & origin=real(0d0), interval=real(output_tick) )
  
  call HistoryAddVariable( &
    & varname='norm1', dims=(/ 't' /), &
    & longname='norm1', units=' ', xtype='double' )

  call HistoryAddVariable( &
    & varname='norm2', dims=(/ 't' /), &
    & longname='norm2', units=' ', xtype='double' )

  call HistoryAddVariable( &
    & varname='norm_infinity', dims=(/ 't' /), &
    & longname='norm_infinity', units=' ', xtype='double' )

  call HistoryAddAttr('norm1', '+glevel', get_glevel(icgrid_ref))
  call HistoryAddAttr('norm1', '+alpha', alpha)
  
  !
  !
  call IcGrid_ncWriter_Init(writer, icgrid_ref)
  call open_ncFile(writer, 'h_error_check.nc')
  call ncdef_Simulation_Parameter(writer, integration_time, DelTime, output_tick)

  ! 格子点情報の定義 
  call ncdef_GridData(writer)
  self%errorh_ncVarID = ncdef_FieldData(writer, self%error_h)

  ! 定義の終了. 
  call end_ncdef(writer)

  ! 格子情報を netCDF に書きこむ. 
  !
  call write_GridData(writer)

end subroutine init_TestCase2

!
subroutine finalize_TestCase2(self)
  ! 宣言文 ; Declaration statement
  !
  type(TestCase2), intent(inout) :: self

  ! 実行文 ; Executable statement
  !

  !
  call HistoryClose()
  call close_ncFile(writer)

end subroutine finalize_TestCase2

!
subroutine set_initial_v(self, init_v)
  ! 宣言文 ; Declaration statement
  !
  type(TestCase2), intent(inout) :: self
  type(Field_IcGrid2D), intent(inout) :: init_v

  ! 実行文 ; Executable statement
  !

  write(*,*) 'init_Testcase2_v'

  call create_solid_rotation_field(self%true_v, angular_speed, alpha)
  init_v%val(:,:,:,:) =  self%true_v%val(:,:,:,:)

end subroutine set_initial_v

!
subroutine set_initial_h(self, init_h)
  ! 宣言文 ; Declaration statement
  !
  type(TestCase2), intent(inout) :: self
  type(Field_IcGrid2D), intent(inout) :: init_h

  !
  ! Work Variables
  !
  integer :: rcID, i,j
  integer :: EMin, EMax
  type(IcGrid2D_FVM), pointer :: icgrid
  real(DP) :: s_pos(3)
  real(DP) :: ic_radius
  real(DP) :: tmp

  ! 実行文　; Executable statement
  !
  
  write(*,*) 'init_Testcase2_h'

  icgrid => get_IcGrid(init_h)
  ic_radius = get_IcRadius(icgrid)
  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
          !
          tmp = dcos(s_pos(1)) * dcos(s_pos(2)) * dsin(alpha) + dsin(s_pos(2)) * dcos(alpha)
        else
          !
          tmp = dsin(s_pos(2)) * dcos(alpha)

        end if
        
        self%true_h%val(rcID,i,j, :) = &
          & h_0 - ( ic_radius * Omega * u_0 + 0.5d0 * u_0**2 ) * tmp**2 / Grav       

     end do
    end do
  end do

  init_h%val(:,:,:,:) = self%true_h%val(:,:,:,:)

end subroutine set_initial_h

!
subroutine set_initial_hs(self, init_hs)

  ! 宣言文 ; Declaration statement
  !
  type(TestCase2), intent(inout) :: self
  type(Field_IcGrid2D), intent(inout) :: init_hs

  ! 実行文 ; Executable statement
  !

  write(*,*) 'Testcase2: initialize height surface field'  
  init_hs%val(:,:,:,:) = 0.0d0

end subroutine set_initial_hs

!
subroutine timelevel_Updated(self, tstep, dt, v_n, h_n)

  ! 宣言文 ; Declaration statement
  !
  type(TestCase2), intent(inout) :: self
  integer, intent(in) :: tstep
  real(DP), intent(in) :: dt
  type(Field_IcGrid2D), intent(inout) :: v_n
  type(Field_IcGrid2D), intent(inout) :: h_n


  ! 実行文 ; Executable statement
  !
  
  if ( mod(tstep * dt, output_tick) == 0 ) then
  !
    call eval_numcal_h_solution( tstep*dt, self%true_h, h_n)
    !
    self%error_h%val(:,:,:,:) = self%true_h%val(:,:,:,:) - h_n%val(:,:,:,:)
    call write_FieldData(writer, self%errorh_ncVarID, self%error_h)
    call increase_recorde_counter(writer)

  end if

end subroutine timelevel_Updated

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 非公開手続き
! 
!!!!!!!!!!!!!!!!!!!!!!!!!

!
subroutine eval_numcal_h_solution(t, true_h, h_n)

  ! 宣言文 ; Declaration statement
  !
  real(DP), intent(in) :: t
  type(Field_IcGrid2D), intent(inout) :: true_h
  type(Field_IcGrid2D), intent(in) :: h_n
  
  ! 作業変数
  ! Work variables
  !
  real(DP) norm1Val, norm2Val, normInfityVal

  ! 実行文 ; Executable statement
  !
  
  !
  norm1Val = numerical_error_norm1(true_h, h_n)
  norm2Val = numerical_error_norm2(true_h, h_n)
  normInfityVal = numerical_error_norminfinity(true_h, h_n, maxmin_info_flag=.false.)

  call HistoryPut('norm1', norm1Val)
  call HistoryPut('norm2', norm2Val)
  call HistoryPut('norm_infinity', normInfityVal)

end subroutine eval_numcal_h_solution

end module class_TestCase2