src/Williamson1992_TestCase/class_TestCase5.f90

説明を見る。

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

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

  !
  !
  !
  use gtool_history

  ! * IGMBaseLib *
  ! *

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

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

  ! 正二十面体格子の管理クラス
  !
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_EffSize_Min, get_EffSize_Max, get_glevel, get_IcRadius, check_pole, &
    & RC_REGIONS_NUM, NOT_POLE_FLAG, NORTH_POLE_FLAG, SOUTH_POLE_FLAG

  ! 正二十面体格子上に分布する物理場の管理クラス
  !
  !
  use Field_IcGrid2D_manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, Field_IcGrid2D_Final, &
    & get_icgrid, paste_field2D_margin


  !
  !
  !
  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

  !
  !
  !
  use field_manager, only: &
    & diff_eval, xy_Coli

  !
  !
  !
  use conservation_analysis, only: &
    & init_conservation_analysis, &
    & calc_total_energy, calc_potential_enstrophy

  ! 宣言文 ; Declaration statement
  !
  implicit none
  private
  
  ! 公開変数
  ! Public variables
  !

  real(DP) :: h_0 

  real(DP), parameter :: hs_0 = 2000.0d0

  real(DP), parameter :: theta_c = PI / 6.0d0

  real(DP), parameter :: lambda_c = - PI / 2.0d0
  
  real(DP), parameter :: cb_R = PI / 9.0d0
 
  real(DP), parameter :: u_0 = 20.0d0

  real(DP) :: angular_speed

  ! クラス定義
  ! Class definition
  !

  !
  type, public :: TestCase5

    ! 宣言文 ; Declaration statement
    !
    integer :: dummy
  end type TestCase5

  !
  !
  !
  interface initialize_TestCase
    module procedure init_TestCase5
  end interface initialize_TestCase

  interface finalize_TestCase
    module procedure finalize_TestCase5
  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) :: init_TE

  real(DP) :: init_PE

  real(DP) :: TE

  real(DP) :: PE

  type(Field_IcGrid2D), pointer :: ini_v

  type(Field_IcGrid2D), pointer :: ini_h

  type(Field_IcGrid2D), pointer :: ini_hs

contains

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

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

  ! 実行文 ; Executable statement
  !

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

  !
  allocate(ini_v)
  allocate(ini_h)
  allocate(ini_hs)

  call Field_IcGrid2D_Init(ini_v, icgrid_ref, 'init_v', 3)
  call Field_IcGrid2D_Init(ini_h, icgrid_ref, 'init_h', 1)
  call Field_IcGrid2D_Init(ini_hs, icgrid_ref, 'init_hs', 1)

  call calc_init_v()
  call paste_field2D_margin(ini_v)
  call calc_init_h()
  call calc_init_hs()

  call init_conservation_analysis(ini_h, ini_hs)
  init_TE = calc_total_energy(ini_h, ini_hs, ini_v)
  init_PE = calc_potential_enstrophy(ini_h, ini_hs, ini_v, xy_Coli, diff_eval)

  write(*,*) 'initTE', init_TE, init_PE
  open(10, file='TE_PE.dat',status='unknown')

  !
  !

  !
  call HistoryCreate( &
    & file='dissipation_rate.nc', title='disspation rates', &
    & source=' ', institution=' ', &
    & dims=(/ 't' /), dimsizes=(/ 0 /), &
    & longnames=(/ 'time' /), units=(/ 's' /), &
    & origin=real(0d0), interval=real(output_tick) )

  call HistoryAddVariable( &
    & varname='TE', dims=(/ 't' /), &
    & longname='total energy', units='m m2 s-2', xtype='double' )

  call HistoryAddVariable( &
    & varname='PE', dims=(/ 't' /), &
    & longname='potential enstrophy', units='s-2 m-1', xtype='double' )

  call HistoryAddVariable( &
    & varname='TE_DISSIP_RATE', dims=(/ 't' /), &
    & longname='dissipation rate of total energy', units='1', xtype='double' )

  call HistoryAddVariable( &
    & varname='PE_DISSIP_RATE', dims=(/ 't' /), &
    & longname='dissipation rate of potential enstrophy', units='1', xtype='double' )

  call HistoryAddAttr('TE', '+glevel', get_glevel(icgrid_ref))
  call HistoryAddAttr('TE', '+alpha', alpha)

end subroutine init_TestCase5

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

  ! 実行文 ; Executable statement
  !
  call HistoryClose()

end subroutine finalize_TestCase5

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

  ! 実行文 ; Executable statement
  !

  write(*,*) 'init_Testcase5_v'

  init_v%val = ini_v%val

  call Field_IcGrid2D_Final(ini_v)
  deallocate(ini_v)

end subroutine set_initial_v

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

  ! 実行文 ; Executable statement
  !
  init_h%val = ini_h%val

  call Field_IcGrid2D_Final(ini_h)
  deallocate(ini_h)

end subroutine set_initial_h

!
subroutine set_initial_hs(self, init_hs)

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

  ! 実行文　; Executable statement
  !
  init_hs%val = ini_hs%val

end subroutine set_initial_hs

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

  ! 宣言文 ; Declaration statement
  !
  type(TestCase5), 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) :: TE
  real(DP) :: PE

  ! 実行文 ; Executable statement
  !
  

  !write(*,*) TE, PE
  if ( mod( tstep, int(output_tick / dt) ) == 0 ) then
    TE = calc_total_energy(h_n, ini_hs, v_n)
    PE = calc_potential_enstrophy(h_n, ini_hs, v_n, xy_Coli, diff_eval)

    call HistoryPut('TE', TE)
    call HistoryPut('PE', PE)

    call HistoryPut('TE_DISSIP_RATE', abs( TE - init_TE ) / init_TE )
    call HistoryPut('PE_DISSIP_RATE', abs( PE - init_PE ) / init_PE )

  end if

end subroutine timelevel_Updated

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

subroutine calc_init_h()
  ! 宣言文 ; Declaration statement
  !

  !
  ! 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
  !

  icgrid => get_IcGrid(ini_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 = cos(s_pos(1)) * cos(s_pos(2)) * sin(alpha) + sin(s_pos(2)) * cos(alpha)
        else
          !
          tmp = sin(s_pos(2)) * cos(alpha)
        end if

        ini_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

end subroutine calc_init_h

!
!
!
subroutine calc_init_hs()
  ! 宣言文 ; Declaration statement
  !

  ! 作業変数
  ! 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) :: r

  ! 実行文　; Executable statement
  !

  write(*,*) 'init_Testcase5_hs'

  icgrid => get_IcGrid(ini_hs)
  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
          s_pos(1) = 0.0d0
        end if

        r = sqrt(min(cb_R**2, (s_pos(1) - lambda_c)**2 + (s_pos(2) - theta_c)**2 ))
!        write(*,*) cb_R**2, (s_pos(1) - lambda_c)**2 + (s_pos(2) - theta_c)**2

        ini_hs%val(rcID,i,j, :) = hs_0 * ( 1.0d0 - r / cb_R )
        if(rcID==1 .and. i == 1 .and. j==1 ) write(*,*) 'hs:', ini_hs%val(rcID,i,j,:), s_pos
     end do
    end do
  end do
end subroutine calc_init_hs

!
!
!
subroutine calc_init_v()

  ! 実行文 ; Executable statement
  !

  write(*,*) 'init_Testcase5_v'

  call create_solid_rotation_field(ini_v, angular_speed, alpha)

end subroutine calc_init_v

end module class_TestCase5