!> Tomita et al(2004), Test case 2(Gravity wave)
!!
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
!! - Test Configuration
!!   - The basic state for thermodynamic variables is a stratified state with a constant Brunt Vaisala frequency \f$ N \f$.
!!   - In this test, there is not diffusion and rotaion of earth.
!!   - As the initial condition, the following perturbation for potential temperature is
!!      superimposed on the basic state.
!!     \f[
!!         \theta' = \Delta \theta f(\lambda, \phi) g(z)
!!     \f]
!!     where \f$\Delta \theta \f$ is the amplituide of potential temperature perturbation.
!!     And the functions \f$ f(\lambda, \theta), g(z) \f$ is same as in Test Case1.
!!
!! - Each SubTest Configuration
!!   - Case 2-1
!!     - ZT/k_max= 10[km]/21, DelTimeLong=1800[s], DelTimeLong/DelTimeShort=4
!!     - Divergence dumping coefficient as \f$ \alpha'_dH=0.05 \f$.
!!     - Brunt Vaisala frequency \f$ N=0.01 \f$ [s-1].
!!     - Number of vertical mode nv=1.
!!   - Case 2-2
!!      - Same as Case 2-1 except for \f$ N=0.02 \f$ [s-1].
!!   - Case 2-3
!!      - Same as Case 2-1 except for \f$ n_v=2 \f$.
!!   - Case 2-4
!!      - Same as Case 2-1 except for \f$ \alpha'_dH=0.05 \f.
!!
!!
module experiment_impl_template

  !* gtool ***************
  !*

  use dc_types, only: &
    & DP, TOKEN, STRING

  !* IGMBaseLib ***********
  !*

  use igmcore_math, only: &
    & PI

  use igmcore_coordinate_conversion, only: &
    & orth_to_geo_pos, geo_to_orth_pos

  use igmcore_spherical_trigonometry, only: &
    & geodesic_arc_length

  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_glevel, get_IcRadius, &
    & get_EffSize_Min, get_EffSize_Max, &
    & get_LocalProcRectInfo, &     
    & RC_REGIONS_NUM

  ! 正二十面格子における広域矩形領域のサブドメインである局所矩形領域(各 MPI プロセスと対応)の情報管理
  ! Manager of information about a local rectangle region (corresponding to each MPI process)
  !
  use LocalProcRect_Manager, only: &
    & LocalProcRectInfo, &
    & get_InteriorRegion_IndexInfo

  use IcGrid3D_FVM_Manager, only: &
    & IcGrid3D_FVM, &
    & get_IcGrid2D, &
    & get_vertical_level_num

  use Field_Grid1D_Manager, only: &
    & Field_Grid1D

  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D

  use Field_IcGrid3D_Manager, only: &
    & Field_IcGrid3D, &
    & get_IcGrid3D, get_field_layer, &
    & fill_value
    
  use TFCoord_IcGrid3D_helper, only: &
    & TFCoord_IcGrid3D, &
    & TF_to_Z_IcGrid3D, calc_gamma

  !* DIGMO ***********
  !*

  use AtmBasicState_Builder, only: &
    & create_ConstBuoFreqAtm_BasicState

  use SimParamList, only: &
    & SysParam => SystemParamList

  implicit none
  private

  public :: experiment_impl_Init, experiment_impl_Final
  public :: set_TopographicData, set_initialFieldData
  public :: LongTimeLevel_updated, ShortTimeLevel_updated

  real(DP), parameter :: N  = 0.01d0
  integer, parameter :: VERTICAL_MODE_NUM = 1

  real(DP), parameter :: T0 = 300d0
  real(DP), parameter :: P0 = 1.0d05

  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax

contains

subroutine experiment_impl_Init( &
  & experimentName &
  & )

  character(*), intent(inout) :: experimentName

end subroutine experiment_impl_Init

subroutine set_TopographicData( &
  & xy_Zs &
  & )

  type(Field_IcGrid2D), intent(inout) :: xy_Zs

end subroutine set_TopographicData

subroutine set_initialFieldData( &
  & z_RhoBasic, z_PressBasic, &
  & xyz_Vh, xyr_W, xyz_RhoPerturb, xyz_PressPerturb, xyz_e, &
  & tfcoordz, tfcoordr &
  & )

  type(Field_Grid1D), intent(inout) :: z_RhoBasic
  type(Field_Grid1D), intent(inout) :: z_PressBasic
  type(Field_IcGrid3D), intent(inout) :: xyz_Vh
  type(Field_IcGrid3D), intent(inout) :: xyr_W
  type(Field_IcGrid3D), intent(inout) :: xyz_RhoPerturb
  type(Field_IcGrid3D), intent(inout) :: xyz_PressPerturb
  type(Field_IcGrid3D), intent(inout) :: xyz_e
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordr

  type(LocalProcRectInfo), pointer :: lcRCInfo
  real(DP), pointer :: grid1D(:)
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(IcGrid3D_FVM), pointer :: icgrid3Dz

  grid1D => z_RhoBasic%grid1D
  icgrid3Dz => get_icgrid3D(xyz_e)
  icgrid2D => get_icgrid2D(icgrid3Dz)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)

  !
  !

  call create_constBuoFreqAtm_BasicState( &
    & z_RhoBasic, z_PressBasic, N, &
    & SysParam%GasConst, SysParam%Grav, SysParam%Cp, SysParam%Cv,  &
    & T0, P0 &
     & )

  call create_PressPerturbation_fields( &
    & xyz_PressPerturb, z_Pressbasic, &
    & xyz_RhoPerturb, z_RhoBasic, &
    & xyz_e, grid1D(:), &
    & icgrid2D, icgrid3Dz, tfcoordz &
    & )

end subroutine set_initialFieldData

subroutine LongTimeLevel_updated( &
  & NStepLong, xyz_RhoBasic, xyz_PressBasic, &
  & xyz_VhNl, xyr_WNl, xyz_RNl, xyz_ENl, xyz_QNl &
  & )

  integer, intent(in) :: NStepLong
  type(Field_IcGrid3D), intent(inout) :: xyz_RhoBasic
  type(Field_IcGrid3D), intent(inout) :: xyz_PressBasic
  type(Field_IcGrid3D), intent(inout) :: xyz_VhNl
  type(Field_IcGrid3D), intent(inout) :: xyr_WNl
  type(Field_IcGrid3D), intent(inout) :: xyz_RNl
  type(Field_IcGrid3D), intent(inout) :: xyz_ENl
  type(Field_IcGrid3D), intent(inout) :: xyz_QNl


end subroutine LongTimeLevel_updated

subroutine ShortTimeLevel_updated( &
  & NStepLong, NStepShort, xyz_RhoBasic, xyz_PressBasic, &
  & xyz_VhDevNs, xyr_WDevNs, xyz_RDevNs, xyz_EDevNs, &
  & xyz_QNl &
  & )

  integer, intent(in) :: NStepLong
  integer, intent(in) :: NStepShort
  type(Field_IcGrid3D), intent(inout) :: xyz_RhoBasic
  type(Field_IcGrid3D), intent(inout) :: xyz_PressBasic
  type(Field_IcGrid3D), intent(inout) :: xyz_VhDevNs
  type(Field_IcGrid3D), intent(inout) :: xyr_WDevNs
  type(Field_IcGrid3D), intent(inout) :: xyz_RDevNs
  type(Field_IcGrid3D), intent(inout) :: xyz_EDevNs
  type(Field_IcGrid3D), intent(inout) :: xyz_QNl

!  write(*,*) 'Timelevel updated experiment--', NStepLong, NStepShort

end subroutine ShortTimeLevel_updated

subroutine experiment_impl_Final()

end subroutine experiment_impl_Final

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine create_PressPerturbation_fields( &
  & xyz_PressPerturb, z_PressBasic, &
  & xyz_RhoPerturb, z_RhoBasic, & 
  & xyz_eNl, AltZ, &
  & icgrid2D, icgrid3Dz, tfcoordz &
  & )

  type(Field_IcGrid3D), intent(out) :: xyz_PressPerturb
  type(Field_Grid1D), intent(in) :: z_PressBasic
  type(Field_IcGrid3D), intent(out) :: xyz_RhoPerturb
  type(Field_Grid1D), intent(in) :: z_RhoBasic
  type(Field_IcGrid3D), intent(out) :: xyz_eNl
  real(DP), intent(in) :: AltZ(:)
  type(IcGrid2D_FVM), intent(in) :: icgrid2D
  type(IcGrid3D_FVM), intent(in) :: icgrid3Dz
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz

  type(Field_IcGrid2D), pointer :: xyz_PressPerturb_k
  type(Field_IcGrid2D), pointer :: xyz_RhoPerturb_k
  type(Field_IcGrid2D), pointer :: xyz_eNl_k

  real(DP) :: s_pos(3), orth_pos0(3)
  real(DP) :: r, Radius, SRadius
  real(DP) :: lambda0, phi0
  real(DP) :: DeltaTheta
  real(DP) :: ThetaBasic, ThetaDash
  real(DP) :: xi, z_rcij, zT
  real(DP) :: gamma_
  real(DP) :: fxy, gxi

  integer :: EMin, EMax, KzMax
  integer :: rcID, i, j, k


  EMin = get_EffSize_Min(icgrid2D)
  EMax = get_EffSize_Max(icgrid2D)
  KzMax = get_vertical_level_num(icgrid3Dz)
  SRadius = get_IcRadius(icgrid2D)

  lambda0 = 0.0d0
  phi0 = 0.0d0
  orth_pos0(:) = geo_to_orth_pos((/ lambda0, phi0, SRadius /))
  DeltaTheta = 10d0
  Radius =  SRadius / 3.0d0
  zT = icgrid3Dz%verticalGrid(KzMax)

  do k=1, KzMax
    xyz_PressPerturb_k => get_field_layer(xyz_PressPerturb, k)
    xyz_RhoPerturb_k => get_field_layer(xyz_RhoPerturb, k)
    xyz_eNl_k => get_field_layer(xyz_eNl, k)

    xi = icgrid3Dz%verticalGrid(k)
    ThetaBasic = T0 * exp(( N**2 / SysParam%Grav ) * AltZ(k))

    do rcID=rcIDMin, rcIDMax
      do j=jEMin, jEMax
        do i=iEMin, iEMax

          s_pos(:) = orth_to_geo_pos(icgrid2D%rcs_AGrid(:, i,j,rcID))
          r = geodesic_arc_length(SRadius, orth_pos0(:), icgrid2D%rcs_AGrid(:, i,j,rcID))

          if ( r < Radius ) then
            fxy = 0.5d0 * ( 1.0d0 + cos( PI * r/Radius ) )
          else
            fxy = 0.0d0
          end if

          z_rcij = TF_to_Z_IcGrid3D(tfcoordz, xi, rcID, i, j)
          gxi = sin( real(VERTICAL_MODE_NUM) * PI * z_rcij / zT )

          ThetaDash = DeltaTheta * fxy * gxi
          
          xyz_PressPerturb_k%val(1, i,j,rcID) = 0d0 
          xyz_RhoPerturb_k%val(1, i,j,rcID) = - z_RhoBasic%val(k,1) * ThetaDash / ( ThetaBasic + ThetaDash )
          
				  	xyz_eNl_k%val(1, i,j,rcID) = SysParam%Cv * ( z_PressBasic%val(k,1) + xyz_PressPerturb_k%val(1, i,j,rcID) ) &
			      	& / ( ( z_RhoBasic%val(k,1) + xyz_RhoPerturb_k%val(1, i,j,rcID) ) * SysParam%GasConst )
     
        end do
      end do
    end do
  end do
end subroutine create_PressPerturbation_fields

end module experiment_impl_template
