src/Williamson1992_TestCase/class_TestCase6.f90

説明を見る。

module class_TestCase6
  
  ! モジュール引用 ; 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, &
    & 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, Field_IcGrid2D_Final, &
    & get_icgrid, paste_field2D_margin
  

  !
  !
  !
  use Field_Pattern_Builder, only: &
    & create_solid_rotation_field


  !
  !
  !
  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 :: ome = 7.848d-06

  real(DP), parameter :: K =  7.848d-06

  real(DP), parameter :: R = 4.0d0
 
  ! クラス定義
  ! Class definition
  !

  !
  type, public :: TestCase6

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

  !
  !
  !
  interface initialize_TestCase
    module procedure init_TestCase6
  end interface initialize_TestCase

  interface finalize_TestCase
    module procedure finalize_TestCase6
  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_TestCase6(self, icgrid_ref )
  ! 宣言文 ; Declaration statement
  !
  type(TestCase6), intent(inout) :: self
  type(IcGrid2D_FVM), intent(in) :: icgrid_ref

  ! 作業変数
  ! Work variables
  !

  ! 実行文 ; Executable statement
  !
  h_0 = 8.0e03

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

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

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

end subroutine finalize_TestCase6

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

  ! 実行文 ; Executable statement
  !

  write(*,*) 'init_TestCase6_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(TestCase6), 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(TestCase6), 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(TestCase6), 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_v()

  ! 実行文 ; Executable statement
  !

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

  ! 実行文　; Executable statement
  !

  write(*,*) 'init_TestCase6_v'

  icgrid => get_IcGrid(ini_h)
  a = 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_u =   a * ome * cos(s_pos(2)) &
             &   + a * K * cos(s_pos(2))**(R-1) * (R*sin(s_pos(2))**2 - cos(s_pos(2))**2) * cos(R*s_pos(1))

           s_v = - a * K * R * cos(s_pos(2))**(R-1) * sin(s_pos(2)) * sin(R * s_pos(1))

           ini_v%val(rcID,i,j,:) = geo_to_orth_vec((/ s_u, s_v, 0.0d0 /), s_pos(:))

        else
          ini_v%val(rcID,i,j,:) = 0.0d0
        end if

      end do
    end do
  end do

end subroutine calc_init_v

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 =   Grav * h_0 + ic_radius**2 * A(s_pos(2)) &
            &   + ic_radius**2 * B(s_pos(2)) * cos(R * s_pos(1)) &
            &   + ic_radius**2 * C(s_pos(2)) * cos(2.0d0 * R * s_pos(1))
        else
          !
          tmp = Grav * h_0
        end if

        ini_h%val(rcID,i,j, :) = tmp / 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_TestCase6_hs'

  ini_hs%val = 0.0d0

end subroutine calc_init_hs

function A(theta) result(val)

  real(DP), intent(in) :: theta
  real(DP) val

  val =    0.5d0 * ome * (2 * Omega + ome) * cos(theta)**2 &
    &    + 0.25d0 * K**2 * cos(theta) **(2*R) * ( &
    &       (R+1.0d0) * cos(theta)**2 + ( 2.0d0 * R**2 - R - 2.0d0 ) - 2.0d0 * R**2 * cos(theta)**(-2) &
    &        )

end function A

function B(theta) result(val)

  real(DP), intent(in) :: theta
  real(DP) val

  val =  2.0d0 * (Omega + ome) * K / ( ( R + 1.0d0 ) * ( R + 2.0d0 ) ) * cos(theta)**R * &
    &       ( ( R**2 + 2.0d0 * R + 2.0d0 ) - ( R + 1.0d0 )**2 * cos(theta)**2 )

end function B

function C(theta) result(val)

  real(DP), intent(in) :: theta
  real(DP) val

  val =  0.25d0 * K**2 * cos(theta)**(2*R) * ( ( R + 1.0d0 ) * cos(theta)**2 - ( R + 2.0d0 ) )

end function C

end module class_TestCase6