!> Williamson, et al (1992) のテストケース 5 を実施するためのクラスを提供するモジュール.  
!
! Authors::   Yuta Kawai
! Version::   $Id: class_TestCase5.f90,v 1.1 2011-07-15 15:31:55 ykawai Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
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 

  !> 同心円型をした孤立山岳の最大の高度. 
  !! The maximum height of an isolated mpountain.   
  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
  !

  !
  !> 球面浅水モデルの標準テスト(Williamson et al, 1992)の TestCase5 を行うクラス. 
  !! 
  !! 
  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

!
!> TestCase5 クラスのオブジェクトを初期化する. 
!!
!! @param[in,out] self       TestCase5 クラスのオブジェクトの参照. 
!! @param[in]     icgrid_ref IcGrid2D_FVM クラスのオブジェクトの参照.
!!
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

!
!> TestCase5 クラスのオブジェクトを最終化する. 
!!
!! @param[in,out] self TestCase5 クラスのオブジェクトの参照. 
!! 
subroutine finalize_TestCase5(self)
  ! 宣言文 ; Declaration statement
  !
  type(TestCase5), intent(inout) :: self

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

end subroutine finalize_TestCase5

!
!> 初期の速度場として, 剛体回転場を設定する(TestCase2 と同じ). 
!!
!! @param[in,out] self   TestCase5 クラスのオブジェクトの参照. 
!! @param[in,out] init_v
!!
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

!
!> 初期の表面高度場を, 初期の速度場である剛体回転場に対して地衡風平衡を満たすように設定する(TestCase2 と同じ). 
!!
!! @param[in,out] self   TestCase5 クラスのオブジェクトの参照. 
!! @param[in,out] init_h
!!
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

!
!> 下部境界の地形の高度場として, 孤立した山岳地形を設定する. 
!!
!! 孤立した山岳は以下の式(Williamson, 1992 の式(134)) で与えられる. 
!! 
!! 
!! @param[in,out] self   TestCase5 クラスのオブジェクトの参照. 
!! @param[in,out] init_hs
!!
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

!
!> タイムレベルの更新された後の処理を行う. 
!!
!! 時間積分ループにおいて, 物理場管理モジュール(field_manager)の update_timeLevel が呼び出された後に, 
!! このサブルーチンは呼び出される. 
!! 
!! @param[in,out] self   TestCase5 クラスのオブジェクトの参照. 
!! @param[in]     tstep
!! @param[in]     dt
!! @param[in,out] v_n
!! @param[in,out] h_n
!!
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
