src/Williamson1992_TestCase/class_TestCase1.f90

説明を見る。

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

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

  !
  !
  !
  use gtool_history

  ! * IGMBaseLib *
  ! *

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

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

  ! 座標変換
  ! 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, &
    & RC_REGIONS_NUM
  
  ! 正二十面体格子上に分布する物理場の管理クラス
  !
  !
  use Field_IcGrid2D_manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, get_icgrid
  
  !
  !
  !
  use Field_Pattern_Builder, only: &
    & create_solid_rotation_field, create_cosine_bell_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, integration_time, DelTime

  ! 宣言文 ; Declaration statement
  !
  implicit none
  private
  
  ! 公開変数
  ! Public variables
  !
  
  real(DP), parameter :: h_0 = 1000.0d0
  
  real(DP), parameter :: init_theta_c = 0.0d0

  real(DP), parameter :: init_lambda_c = 3.0d0 * PI / 2.0d0
  
  real(DP) :: cb_R
 
  real(DP) :: u_0
  
  ! クラス定義
  ! type definition
  !

  !
  type, public :: TestCase1

    ! 宣言文 ; 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

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

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

  end type TestCase1

  !
  !
  !
  interface initialize_TestCase
    module procedure init_TestCase1
  end interface initialize_TestCase

  interface finalize_TestCase
    module procedure finalize_TestCase1
  end interface finalize_TestCase

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

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

contains

!
subroutine init_TestCase1(self, icgrid_ref )

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

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

  ! 実行文 ; Executable statement
  !

  ic_radius = get_IcRadius(icgrid_ref)
  
  cb_R = ic_radius / 3.0d0
  u_0 = 2 * PI * ic_radius / ( 12.0d0 * 24.0d0 * 3600d0 )
  angular_speed = 2 * 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%solid_rot_field, icgrid_ref, 'solid_rotation_field', 3)
  call create_solid_rotation_field(self%solid_rot_field, angular_speed, alpha)

  t_lambdaC = init_lambda_c
  t_thetaC = init_theta_c

  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_TestCase1

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

  ! 実行文 ; Executable statement
  !

  !
  call HistoryClose()
  !
  call close_ncFile(writer)

end subroutine finalize_TestCase1

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

  ! 実行文 ; Executable statement
  !

  write(*,*) 'init_Testcase1_v'
  init_v%val(:,:,:,:) =  self%solid_rot_field%val(:,:,:,:)
 
end subroutine set_initial_v

!!
!! @param[in,out] self   TestCase1 クラスのオブジェクトの参照. 
!! @param[in,out] init_h
!!
subroutine set_initial_h(self, init_h)
  ! 宣言文 ; Declaration statement
  !
  type(TestCase1), intent(inout) :: self
  type(Field_IcGrid2D), intent(inout) :: init_h

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

  !
  call create_cosine_bell_field(init_h, init_lambda_c, init_theta_c, cb_R, h_0)
  !call eval_numcal_h_solution(0d0, self%true_h, init_h)


end subroutine set_initial_h

!
subroutine set_initial_hs(self, init_hs)

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

  ! 実行文 ; Executable statement
  !

  write(*,*) 'Testcase1 module - 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(TestCase1), 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

  ! 作業変数
  ! Work variables
  !
  real(DP) :: new_bellC(3), orth_bellC_onXY(3)
  type(IcGrid2D_FVM), pointer :: icgrid
  real(DP) :: ic_radius

  ! 実行文 ; Executable statement
  !
  
  ! 速度場を剛体回転場で上書きする. 
  v_n%val(:,:,:,:) = self%solid_rot_field%val(:,:,:,:)

  !
  !
  icgrid => get_icgrid(v_n)
  ic_radius = get_IcRadius(icgrid)

  if ( tstep /= 0 ) then
    orth_bellC_onXY(:)  = &
      & rotateZ( &
        & rotateY( geo_to_orth_pos( (/ t_lambdaC, t_thetaC, ic_radius /) ), -alpha ), &
        & 2.0d0*PI * dt / ( 12.0d0 * 24.0d0 * 3600d0 )  &
        & )

    new_bellC(:) = orth_to_geo_pos(rotateY( orth_bellC_onXY(:), alpha))
    t_lambdaC = new_bellC(1)
    t_thetaC = new_bellC(2)
  !   write(*,*) '-----', toDegrees(t_lambdaC), toDegrees(t_thetaC)
  end if

  if ( mod(tstep, int(output_tick / dt)) == 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
  !
  
  !
  call create_cosine_bell_field(true_h, t_lambdaC, t_thetaC, cb_R, h_0)
  
  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.)
!  write(*,*) 'h', true_h%val(10,33,33,:), h_n%val(10,33,33,:)

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

end subroutine eval_numcal_h_solution

end module class_TestCase1