src/Williamson1992_TestCase/class_TestCase3.f90

説明を見る。

module class_TestCase3
  
  ! モジュール引用 ; 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, geo_to_orth_vec

  ! 正二十面体格子の管理クラス
  !
  !
  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, 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

  ! OpenMP
  !
  !
  use omp_lib

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

  real(DP) :: u_0

  real(DP) :: angular_speed

  ! クラス定義
  ! Class definition
  !

  !
  type, public :: TestCase3

    ! 宣言文 ; Declaration statement
    !

    ! 非公開メンバ変数
    ! Private member variable
    !
    type(Field_IcGrid2D) :: true_h
    type(Field_IcGrid2D) :: error_h
    type(Field_IcGrid2D) :: true_v
    type(Field_IcGrid2D) :: error_v
    integer :: errorh_ncVarID
    integer :: errorV_ncVarID

    ! f2003 code
    !type(Field_IcGrid2D), private :: true_h
    !type(Field_IcGrid2D), private :: error_h
    !type(Field_IcGrid2D), private :: true_v
    !type(Field_IcGrid2D), private :: error_v
    !integer, private :: errorh_ncVarID
    !integer, private :: errorV_ncVarID

  end type TestCase3

  !
  !
  !
  interface initialize_TestCase
    module procedure init_TestCase3
  end interface initialize_TestCase

  interface finalize_TestCase
    module procedure finalize_TestCase3
  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_TestCase3(self, icgrid_ref )
  ! 宣言文 ; Declaration statement
  !
  type(TestCase3), 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)
  call Field_IcGrid2D_Init(self%error_v, icgrid_ref, 'error_v', 3)

  !
  call create_fields(self%true_h, self%true_v)

  !
  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='v_norm1', dims=(/ 't' /), &
    & longname='v_norm1', units=' ', xtype='double' )

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

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

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

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

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


  call HistoryAddAttr('v_norm1', '+glevel', get_glevel(icgrid_ref))
  call HistoryAddAttr('v_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)
  self%errorv_ncVarID = ncdef_FieldData(writer, self%error_v)

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

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

end subroutine init_TestCase3

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

  ! 実行文 ; Executable statement
  !

  !
  call HistoryClose()
  call close_ncFile(writer)

end subroutine finalize_TestCase3

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

  ! 実行文 ; Executable statement
  !

  write(*,*) 'init_TestCase3_v'

  init_v%val(:,:,:,:) =  self%true_v%val(:,:,:,:)

end subroutine set_initial_v

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

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

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

end subroutine set_initial_h

!
subroutine set_initial_hs(self, init_hs)

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

  ! 実行文 ; Executable statement
  !

  write(*,*) 'TestCase3: 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(TestCase3), 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, self%true_v, h_n, v_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, true_v, &
  & h_n, v_n &
  & )

  ! 宣言文 ; Declaration statement
  !
  real(DP), intent(in) :: t
  type(Field_IcGrid2D), intent(in) :: true_h
  type(Field_IcGrid2D), intent(in) :: true_v
  type(Field_IcGrid2D), intent(in) :: h_n
  type(Field_IcGrid2D), intent(in) :: v_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('h_norm1', norm1Val)
  call HistoryPut('h_norm2', norm2Val)
  call HistoryPut('h_norm_infinity', normInfityVal)

  norm1Val = numerical_error_norm1(true_v, v_n)
  norm2Val = numerical_error_norm2(true_v, v_n)
  normInfityVal = numerical_error_norminfinity(true_v, v_n, maxmin_info_flag=.false.)

  call HistoryPut('v_norm1', norm1Val)
  call HistoryPut('v_norm2', norm2Val)
  call HistoryPut('v_norm_infinity', normInfityVal)

end subroutine eval_numcal_h_solution

subroutine create_fields(h_field, v_field)
  ! モジュール引用 ; Use statements
  !

  use igmcore_math, only: toDegrees
  use igmcore_linear_algebra, only: rotateY

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid2D), intent(inout) :: h_field
  type(Field_IcGrid2D), intent(inout) :: v_field

  ! 作業変数
  ! Work variables
  !
  type(IcGrid2D_FVM), pointer :: icgrid
  integer :: rcID, i, j
  integer :: EMin, EMax
  integer :: pole_flag
  integer :: ndiv
  
  real(DP) :: u, v, u_ds
  real(DP) :: ic_radius, s_pos(3), s_pos_dash(3)
  integer, parameter :: MAX_GLEVEL = 9
  real(DP) :: integ_dtheta

  !
 
  ! 実行文 ; Executable statements
  !
  
  !
  icgrid => get_icgrid(v_field)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)
  ic_radius = get_IcRadius(icgrid)
  s_pos_dash(3) = ic_radius

  !
  integ_dtheta = 2.0d0 * PI / ( 10.0d0 * 2**(MAX_GLEVEL - 1) ) / 100.0

  !
  write(*,*) 'Create the fields ..'

  do rcID=1, RC_REGIONS_NUM
    !$omp parallel do private(s_pos, s_pos_dash, u_ds, u, v, ndiv)
    do j=EMin, EMax
      do i=EMin, EMax
        pole_flag = check_pole(icgrid, rcID,i,j)
        if(  pole_flag == NORTH_POLE_FLAG ) then
          s_pos(:) = (/ PI, PI / 2.0d0, ic_radius /)
        else if( pole_flag == SOUTH_POLE_FLAG ) then
          s_pos(:) = (/ 0.0d0, - PI / 2.0d0, ic_radius /)
        else
          s_pos(:) = orth_to_geo_pos(icgrid%rcs_AGrid(rcID,i,j,:))
        end if

        !
        s_pos_dash(2) = &
          & dasin( &
          &   dsin(s_pos(2)) * dcos(alpha) + dsin(alpha) * dcos(s_pos(1)) * dcos(s_pos(2)) &
          & )

        s_pos_dash(1) = dasin( dsin(s_pos(1)) * dcos(s_pos(2)) / dcos(s_pos_dash(2)) )

        !
        if ( dcos(alpha) * dcos(s_pos(1)) * dcos(s_pos(2)) - dsin(alpha) * dsin(s_pos(2)) < 0.0d0 ) then
          s_pos_dash(1) = PI - s_pos_dash(1)         
        end if

        !
        u_ds = calc_u_dash(s_pos_dash(2)) 
        v = u_ds * dsin(alpha) * dsin(s_pos_dash(1)) / dcos(s_pos(2))
        u = v * dsin(s_pos(2)) * dtan(s_pos(1)) + u_ds * dcos(s_pos_dash(1)) / dcos(s_pos(1))

        if ( (abs( s_pos(2)) - PI / 2.0)**2 < 1.0e-20 ) then
          v_field%val(rcID,i,j,:) = geo_to_orth_vec((/ u, 0.0d0, 0.0d0 /), &
            &                        (/ PI, s_pos_dash(2), ic_radius /) )
        else
          !
          v_field%val(rcID,i,j, :) = geo_to_orth_vec((/ u, v, 0.0d0 /), s_pos)
        end if

        ndiv = int( ( s_pos_dash(2) - (- 0.5d0 * PI ) ) / integ_dtheta )

        !
        h_field%val(rcID,i,j,:) = calc_h(s_pos_dash(2), ic_radius, ndiv)
      end do
    end do
  end do

  write(*,*) 'Finish creating the fields'

end subroutine create_fields

function calc_h(theta_dash, ic_radius, n) result(val)
  real(DP), intent(in) :: theta_dash
  real(DP), intent(in) :: ic_radius
  integer, intent(in) :: n
  real(DP) val

  !
  !
  !
  real(DP) :: integ_ret
  real(DP) :: tau_lbound, tau
  real(DP) :: h, h2
  real(DP) :: e_sum, o_sum
  integer :: i

  !
  !
  if( n == 0 ) then
    val = h_0
    return
  end if

  tau_lbound = - PI / 2.0d0
  h = ( theta_dash - tau_lbound ) / dble(n)

  e_sum = 0.0d0
  do i=2, n-2, 2
    tau = tau_lbound + i * h
    e_sum = e_sum + h_integrand(tau, ic_radius)
  end do

  o_sum = 0.0d0
  do i=1, n-1, 2
    tau = tau_lbound + i * h
    o_sum = e_sum + h_integrand(tau, ic_radius)
  end do

  integ_ret =   h * ( h_integrand(tau_lbound, ic_radius) &
    &         + h_integrand(theta_dash, ic_radius) + 2.0d0 * e_sum + 4.0d0 * o_sum ) / 3.0d0

  val = h_0 - ic_radius * integ_ret / Grav
       
end function calc_h

function h_integrand(tau, a) result(val)
  !
  !
  real(DP), intent(in) :: tau
  real(DP), intent(in) :: a
  real(DP) val

  !
  !
  real(DP) :: u_ds

  !
  !
  u_ds = calc_u_dash(tau)

  val = u_ds * ( 2.0d0 * Omega * dsin(tau) + u_ds * dtan(tau) / a )
  
end function h_integrand

function calc_u_dash(theta_dash) result(val)
  real(DP), intent(in) :: theta_dash
  real(DP) val

  !
  !
  !
  real(DP), parameter :: theta_b = - PI / 6.0d0
  real(DP), parameter :: theta_e = PI / 2.0d0
  real(DP), parameter :: x_e = 0.3d0
  real(DP) :: x

  !
  !
  x = x_e * ( theta_dash - theta_b ) / ( theta_e - theta_b )
   
  val = u_0 * b(x) * b(x_e - x) * dexp(4.0d0 / x_e)
       
end function calc_u_dash

function b(x) result(val)
  !
  !
  real(DP), intent(in) :: x
  real(DP) val

  !
  !

  if( x <= 0.0d0 ) then
    val = 0.0d0
  else
    val = dexp(-1.0d0/x) 
  end if

end function b

end module class_TestCase3