!> \~japanese 
!! \brief Tomita et al(2004), Test case 3(山岳波)
!! \~english  
!! \brief Tomita et al(2004), Test case 3(Mountain wave)
!!
!!
!! \~
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
!! \~japanese
!! - テストの設定
!!   - 山岳の形状は, 以下の理想的なベル型で与える. 
!!     \f[
!!       h(\lambda,\theta) = h_0 \dfrac{1}{1+(r/d)^2}
!!     \f]
!!     ここで, \f$h_0\f$ は山岳の中心の高度, \f$d\f$ はその山岳の半値幅, \f$r\f$ 中心からの距離である.   
!!   - 帯状流の初期条件として, 以下の剛体回転場を課す. 
!!     \f[
!!       \bar{u}(\phi) = u_0 \cos{\phi}
!!     \f]
!!      ここで, \f$u_0\f$ 赤道の風速である. 
!!   - 他の変数に対しては, 一定のブラントバイサラ振動数\f$N\f$ である静水圧バランスと地衡流バランスが課される. 
!!   - このテストにおいて, \f$|\Dvect{\Omega}|=0\f$とし回転を無視する. 
!!   - 上端境界における重力波の反射を抑制するために, 速度に対してレイリー減衰, 温度に対して基本場へのニュートン冷却を適用する. 
!!     その e-folding 時間は, 以下で与えられる. 
!!    \f$\xi > \xi_0\f$に対して, 
!!     \f[
!!        t^{e}(\xi) = t^{e}_{\textrm{max}} \dfrac{2}{1-\cos(\pi \dfrac{\xi-\xi_0}{z_T-\xi_0})}  
!!     \f]
!!     一方, \f$xi_0\f$の下層に対しては, \f$t^{e}\f$ を無限大とする.
!!     このテストにおいて, \f$t^{e}_{max} = 1800\f$[s], \f$z_T = 40\f$[km], および\f$\xi_0 = 25\f$[km]とする. 
!!
!! - 各サブテストの設定
!!   - Case 3-1
!!     - ZT/k_max= 40[km]/81, DelTimeLong=1800[s], DelTimeLong/DelTimeShort=4
!!     - Divergence dumping 係数 \f$ \alpha'_dH=0.20 \f$.
!!     - ブラントバイサラ係数 \f$ N=0.0187 \f$ [s-1].
!!     - 無次元水平拡散係数 \f$ \gamma_h=1.25 \times 10^{-2} \f$.  
!!     - 無次元鉛直拡散係数 \f$ \gamma_v=3.125 \times 10^{-4} \f$.  
!!   - Case 3-2
!!      - \f$ h0=4000 \f$ [m]以外 Case 3-1 と同じ.
!!
!! \~english
!! - Test Configuration
!!   - The shape of mountain is given as an ideal bell shape:
!!     \f[
!!       h(\lambda,\theta) = h_0 \dfrac{1}{1+(r/d)^2}
!!     \f]
!!     where \f$h_0\f$ is the height at the center of the mountain, \f$d\f$ is the half-width of the mountain, 
!!     and \f$r\f$ is the distance from the center.   
!!   - As the initial condition for the zonal flow, the solid body rotation field is imposed as
!!     \f[
!!       \bar{u}(\phi) = u_0 \cos{\phi}
!!     \f]
!!      where \f$u_0\f$ is the wind speed at the equator. 
!!   - For other variables, the hydrostatic balance with a constant Brunt Vaisala frequency \f$N\f$ and 
!!     the geostrophic balance is imposed. 
!!   - In this test, there is no rotaion of earth so that \f$|\Dvect{\Omega}|=0\f$. 
!!   - In order to avoid the reflection of gravity wave at the top boundary, the Rayleigh damping for the velocity
!!     and the Newtonian cooling to the basic state for the temperature are imposed. The e-folding time is given:
!!     for \f$\xi > \xi_0\f$, 
!!     \f[
!!        t^{e}(\xi) = t^{e}_{\textrm{max}} \dfrac{2}{1-\cos(\pi \dfrac{\xi-\xi_0}{z_T-\xi_0})}  
!!     \f]
!!     while for the layer below \f$xi_0\f$, \f$t^{e}\f$ is set to infinity.
!!     In this test, \f$t^{e}_{max} = 1800\f$[s], \f$z_T = 40\f$[km], and \f$\xi_0 = 25\f$[km]. 
!!
!! - Each SubTest Configuration
!!   - Case 3-1
!!     - ZT/k_max= 40[km]/81, DelTimeLong=1800[s], DelTimeLong/DelTimeShort=4
!!     - Divergence dumping coefficient as \f$ \alpha'_dH=0.20 \f$.
!!     - Brunt Vaisala frequency \f$ N=0.0187 \f$ [s-1].
!!     - Nondimensional horizontal diffusion coefficient \f$ \gamma_h=1.25 \times 10^{-2} \f$.  
!!     - Nondimensional vertical diffusion coefficient \f$ \gamma_v=3.125 \times 10^{-4} \f$.  
!!   - Case 3-2
!!      - Same as Case 3-1 except for \f$ h0=4000 \f$ [m].
!!
!!
module experiment_impl_template

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

  use igmcore_math, only: &
    & PI

  use igmcore_coordinate_conversion, only: &
    & orth_to_geo_pos, geo_to_orth_pos, &
    & geo_to_orth_vec

  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, check_pole, &
    & NORTH_POLE_FLAG, SOUTH_POLE_FLAG, &
    & NOT_POLE_FLAG, &
    & 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, &
    & Field_IcGrid2D_Init, Field_IcGrid2D_Final, &
    & get_icgrid, global_mean, get_field2d_name

  use Field_IcGrid3D_Manager, only: &
    & Field_IcGrid3D, &
    & Field_IcGrid3D_Init, Field_IcGrid3D_Final, &
    & get_IcGrid3D, get_field_layer, get_field3d_name, &
    & fill_value, copy_field3D

  use TFCoord_IcGrid3D_helper, only: &
    & TFCoord_IcGrid3D, &
    & TF_to_Z_IcGrid3D, interpolate_FieldZ1D_into_FieldTF3D, &
    & calc_gamma, global_mean_Field_TFIcGrid3Dz

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

  use AtmBasicState_Builder, only: &
    & create_ConstBuoFreqAtm_BasicState

  use SimParamList, only: &
    & SysParam => SystemParamList, &
    & TimeIntegParam => TimeIntegParamList

  ! 宣言文
  ! Declaration statements
  implicit none
  private

  ! 公開手続き
  ! Public procedures
  !
  public :: experiment_impl_Init, experiment_impl_Final
  public :: set_TopographicData, set_initialFieldData
  public :: LongTimeLevel_RKStep_Updated, LongTimeLevel_Component_Updated, ShortTimeLevel_updated

  ! 非公開変数
  ! Private variables
  !
  
  !** 大気の基本場を特徴づけるのパラメータ
  !* Parameters which characterize the basic atmosphereic state. 
  !* 
  
  !> \~japanese 一定である浮力振動数の値 [s-1]. 
  !! \~english The value of constant buoyancy frequency [s-1].  
  real(DP), parameter :: N  = 0.0187d0
  
  !> \~japanese 等温大気の温度 [K]. 
  !! この温度は, 既知である一定な浮力振動数の値と一貫性があるように決められる, すなわち \f$ T_0 = g^2 / (N^2 * Cp) \f$. 
  !!
  !> \~english The temperature of isothermal atmosphere [K]. 
  !! This temperture will be determined such that it is consistent with the known value of constant buoyancy frequency. 
  !! i.e, \f$ T_0 = g^2 / (N^2 * Cp) \f$. 
  real(DP) :: T0

  !> \~japanese 基準圧力(極における圧力) [Pa]. 
  !! \~english  The reference pressure which is pressure at the equator [Pa]. 
  real(DP), parameter :: P_pole = 1.0d05
  
  !> \~japanese 赤道における剛体回転速度 [m/s]. 
  !! \~english The velocity of solid body rotaion at the equator [m/s]. 
  real(DP), parameter :: u0 = 40d0
  
  !** 地形の高度場を特徴づけるパラメータ
  !* The parameters which characterize the height field of mountain.  
  !*
  
  real(DP), parameter :: h0 = 1000d0
  
  real(DP), parameter :: d  = 1250d03
  
  
  real(DP) :: lambda0 = - PI / 2d0
  real(DP) :: phi0    =   0d0 

  !
  !
    
  real(DP), pointer :: dumpCoefDist_z(:)
  real(DP), pointer :: dumpCoefDist_r(:)  
  real(DP) :: te_max = 1800d0
  real(DP) :: dump_criticalLevel = 25d03
  
  !
  !

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

  type(TFCoord_IcGrid3D), pointer :: tfcoord_r, tfcoord_z
  type(Field_IcGrid3D), pointer :: xyz_ENl_ref, xyz_RhoBasic_ref, & 
    &                              xyz_RNl_ref, xyz_VhNl_ref, xyr_WNl_ref, &
    &                              xyz_QNl_ref
  
contains

subroutine experiment_impl_Init( &
  & experimentName &
  & )

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

end subroutine experiment_impl_Init

!
! 経度 lambda0, 緯度 phi0 に中心を持つローレンツ型の山を, 下部地形として設定する. 
!
subroutine set_TopographicData( &
  & xy_Zs &
  & )

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

  type(LocalProcRectInfo), pointer :: lcRCInfo
  type(IcGrid2D_FVM), pointer :: icgrid2D
  
  integer :: rcID, i, j
  real(DP) :: orth_pos0(3), ic_radius
  real(DP) :: r

  icgrid2D => get_icgrid(xy_Zs)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)
  
  ic_radius = get_IcRadius(icgrid2D) 
  orth_pos0(:) = geo_to_orth_pos((/ lambda0, phi0, ic_radius /))

  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jEMax
      do i=iEMin, iEMax
      
        r = geodesic_arc_length(ic_radius, orth_pos0(:), icgrid2D%rcs_AGrid(:, i,j,rcID))
        xy_Zs%val(1, i,j,rcID) = h0 / ( 1d0 + ( r / d )**2 )
        
      end do
    end do
  end do
  
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), target, intent(in) :: tfcoordz
  type(TFCoord_IcGrid3D), target, intent(in) :: tfcoordr

  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(IcGrid3D_FVM), pointer :: icgrid3Dz, icgrid3Dr

  icgrid3Dz => get_icgrid3D(xyz_e)
  icgrid3Dr => get_icgrid3D(xyr_W)  
  icgrid2D => get_icgrid2D(icgrid3Dz)
  tfcoord_r => tfcoordr
  tfcoord_z => tfcoordz
  
  allocate(dumpCoefDist_z(size(icgrid3Dz%verticalGrid)))
  allocate(dumpCoefDist_r(size(icgrid3Dr%verticalGrid)))
    
  call calc_dumpCoefDist( &
    & dumpCoefDist_z, icgrid3Dz%verticalGrid &
    & )

  call calc_dumpCoefDist( &
    & dumpCoefDist_r, icgrid3Dr%verticalGrid &
    & )

  ! 初期場を設定する. 
  ! Configure the intial field. 
  !
  
  ! (鉛直)等温大気の温度として, 与えられた一定の浮力振動数 N と一貫性のある温度を選ぶ. 
  T0 =  SysParam%Grav * ( SysParam%Grav ) / ( N**2 * SysParam%Cp ) 
  
  ! 浮力振動数 N が一定である基本場を生成し, 基本密度場 z_RhoBasic と 基本圧力場 z_PressBasic を設定する. 
  !
  call create_constBuoFreqAtm_BasicState( &
    & z_RhoBasic, z_PressBasic, N, &
    & SysParam%GasConst, SysParam%Grav, SysParam%Cp, SysParam%Cv,  &
    & T0, P_pole &
     & )

  ! 剛体回転場を生成し, 初期の速度場 xyz_vh を設定する.  
  ! また, 剛体回転場と地衡風平衡にある圧力場を求め, 圧力と密度の摂動場 xyz_PressPerturb, xyz_RhoPerturb を設定する. 
  call create_ZonalFlow_PressPerturbation_fields( &
    & xyz_PressPerturb, z_Pressbasic, &
    & xyz_RhoPerturb, z_RhoBasic, &
    & xyz_e, xyz_vh, &
    & icgrid2D, icgrid3Dz, tfcoordz &
    & )

end subroutine set_initialFieldData

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


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

  xyz_VhNl_ref => xyz_VhNl
  xyr_WNl_ref => xyr_WNl
  xyz_RNl_ref => xyz_RNl
  xyz_RhoBasic_ref => xyz_RhoBasic
  xyz_ENl_ref => xyz_ENl
  xyz_QNl_ref => xyz_QNl
     
end subroutine LongTimeLevel_RKStep_Updated

subroutine LongTimeLevel_Component_Updated( &
  & NStepLong, RKStep,  &
  & xyz_GRNl, xyz_GVhNl, xyr_GWNl, xyz_GEaNl &
  & )

  integer, intent(in) :: NStepLong
  integer, intent(in) :: RKStep
  type(Field_IcGrid3D), intent(inout) :: xyz_GRNl
  type(Field_IcGrid3D), intent(inout) :: xyz_GVhNl
  type(Field_IcGrid3D), intent(inout) :: xyr_GWNl
  type(Field_IcGrid3D), intent(inout) :: xyz_GEaNl
  

  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(Field_IcGrid2D), pointer :: xyz_VhNl_k, xyr_WNl_k, xyz_ENl_k, &
    &                              xyz_RNl_k, xyz_RhoBasic_k

  integer :: rcID, i, j, k
  integer :: KzMax, KrMax
  real(DP) :: rhoNl_rijk, z_rcij
  real(DP) :: s_pos(3), a
  
  real(DP) :: u_rijk(3), T_phi_rijk, DelE, gm_Q(1)
  real(DP) :: VhNl_ori(3)

  icgrid3Dz => get_icgrid3D(xyz_VhNl_ref)    
  KzMax = get_vertical_level_num(icgrid3Dz)
  KrMax = get_vertical_level_num(get_icgrid3D(xyr_WNl_ref))
  icgrid2D => get_icgrid2D(icgrid3Dz)
  a = get_IcRadius(icgrid2D)
  
  do k=2, KrMax-1
    xyr_WNl_k => get_field_layer(xyr_WNl_ref, k) 

    xyr_GWNl%field2Ds(k)%val(:,:,:,:) =   xyr_GWNl%field2Ds(k)%val(:,:,:,:) &
      &                                 - dumpCoefDist_r(k) * xyr_WNl_k%val(:,:,:,:)     
  end do

  do k=1, KzMax
    xyz_RhoBasic_k => get_field_layer(xyz_RhoBasic_ref, k)
    xyz_RNl_k => get_field_layer(xyz_RNl_ref, k)
    xyz_VhNl_k => get_field_layer(xyz_VhNl_ref, k)
    xyz_ENl_k => get_field_layer(xyz_ENl_ref, k)
    
    do rcID=rcIDMin, rcIDMAx
      do j=jEMin, jEMax
        do i=iEMin, iEMax
          
          z_rcij = TF_to_Z_IcGrid3D(tfcoord_z, icgrid3Dz%verticalGrid(k), rcID,i,j)
          s_pos(:) = orth_to_geo_pos(icgrid2D%rcs_AGrid(:, i,j,rcID))
          s_pos(3) = a + z_rcij
            
          u_rijk(:) = geo_to_orth_vec((/ calc_u_phi(s_pos(2)), 0d0, 0d0 /), s_pos(:))
          if ( check_pole(icgrid2D, rcID, i, j) /= NOT_POLE_FLAG ) then
            u_rijk(:) = 0d0
          end if

          rhoNl_rijk = xyz_RhoBasic_k%val(1, i,j,rcID) + xyz_RNl_k%val(1, i,j,rcID)

          xyz_GVhNl%field2Ds(k)%val(1:3, i,j,rcID) =  & 
              &    xyz_GVhNl%field2Ds(k)%val(1:3, i,j,rcID) &
              &  - dumpCoefDist_z(k) * ( xyz_VhNl_k%val(1:3, i,j,rcID) - rhoNl_rijk * u_rijk(:) ) 

                        
          T_phi_rijk = calc_T_phi(s_pos(2), calc_u_phi(s_pos(2)), &
            &            a, SysParam%Grav, SysParam%Omega, SysParam%GasConst/SysParam%Cp, SysParam%GasConst)
                      
          DelE = ( xyz_ENl_k%val(1, i,j,rcID) - rhoNl_rijk * SysParam%Cv * T_phi_rijk )
          xyz_GEaNl%field2Ds(k)%val(1, i,j,rcID) =  xyz_GEaNl%field2Ds(k)%val(1, i,j,rcID) & 
            &                           - dumpCoefDist_z(k) * DelE 

         
         xyz_QNl_ref%field2Ds(k)%val(1,i,j,rcID) =   - dumpCoefDist_z(k) * DelE
            
          if (k>60 .and. rcID==1.and.i==1.and.j==1)then
        write(*,*)  k,xyz_ENl_k%val(1, i,j,rcID)/(rhoNl_rijk * SysParam%Cv),&
        & T_phi_rijk, xyz_RNl_k%val(1, i,j,rcID), xyz_QNl_ref%field2Ds(k)%val(1,i,j,rcID)
           end if
                      
        end do
      end do
    end do
    
  end do

end subroutine LongTimeLevel_Component_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


end subroutine ShortTimeLevel_updated

subroutine experiment_impl_Final()

  deallocate(dumpCoefDist_z, dumpCoefDist_r)
  
end subroutine experiment_impl_Final

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

subroutine calc_dumpCoefDist( &
  & dumpCoefDist, Xi &
  & )
  
  real(DP), intent(inout) :: dumpCoefDist(:)
  real(DP), intent(in) :: Xi(:)
  
  integer :: k, kMax
  real(DP) :: Xi0, Z_T
  
  kMax = size(Xi)
  Z_T = SysParam%domain_top
  Xi0 = dump_criticalLevel
  
  do k=1, kMax
    
    if( Xi(k) > Xi0 ) then
      dumpCoefDist(k) = 0.5d0 * ( 1d0 - cos(PI * ( Xi(k) - Xi0 )/( Z_T - Xi0 )) ) / te_max
    else
      dumpCoefDist(k) = 0d0
    end if
    
  end do
  
end subroutine calc_dumpCoefDist

subroutine create_ZonalFlow_PressPerturbation_fields( &
  & xyz_PressPerturb, z_PressBasic, &
  & xyz_RhoPerturb, z_RhoBasic, & 
  & xyz_eNl, xyz_vhNl, &
  & 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
  type(Field_IcGrid3D), intent(out) :: xyz_vhNl
  type(IcGrid2D_FVM), intent(in) :: icgrid2D
  type(IcGrid3D_FVM), intent(in) :: icgrid3Dz
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz

  type(Field_IcGrid3D) :: xyz_RhoBasic, xyz_PressBasic
  type(Field_IcGrid2D), pointer :: xyz_PressBasic_k, xyz_PressPerturb_k, &
    &                              xyz_RhoBasic_k, xyz_RhoPerturb_k
  type(Field_IcGrid2D), pointer :: xyz_eNl_k, xyz_vhNl_k

  real(DP) :: s_pos(3)
  real(DP) :: xi, z_rcij, zT

  real(DP) :: a, Omega, g
  real(DP) :: Rd, kapper

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

  real(DP) :: p_rijk, p_phi_rijk, rho_rijk, T_phi_rijk, u_rijk
  real(DP) :: factorz

  KzMax = get_vertical_level_num(icgrid3Dz)

  a = get_IcRadius(icgrid2D)
  Omega = SysParam%Omega
  g = SysParam%Grav
  Rd = SysParam%GasConst
  kapper = Rd / SysParam%Cp

  call Field_IcGrid3D_Init(xyz_RhoBasic, icgrid3Dz, 'xyz_RhoBasic', 1)
  call Field_IcGrid3D_Init(xyz_PressBasic, icgrid3Dz, 'xyz_PressBasic', 1)

  call interpolate_FieldZ1D_into_FieldTF3D( &
    & tfcoordz, z_RhoBasic, xyz_RhoBasic &
    & )

  call interpolate_FieldZ1D_into_FieldTF3D( &
    & tfcoordz, z_PressBasic, xyz_PressBasic &
    & )

  do k=1, KzMax

    xyz_PressBasic_k => get_field_layer(xyz_PressBasic, k)
    xyz_PressPerturb_k => get_field_layer(xyz_PressPerturb, k)
    xyz_RhoBasic_k => get_field_layer(xyz_RhoBasic, k)
    xyz_RhoPerturb_k => get_field_layer(xyz_RhoPerturb, k)
    xyz_vhNl_k => get_field_layer(xyz_vhNl, k)
    xyz_eNl_k => get_field_layer(xyz_eNl, k)

    xi = icgrid3Dz%verticalGrid(k)

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

          z_rcij = TF_to_Z_IcGrid3D(tfcoordz, xi, rcID, i, j)
          s_pos(:) = orth_to_geo_pos(icgrid2D%rcs_AGrid(:, i,j,rcID))
          s_pos(3) = a + z_rcij

          u_rijk = calc_u_phi(s_pos(2))
          p_phi_rijk = calc_p_phi(s_pos(2), a, g, Omega, kapper)
          T_phi_rijk = calc_T_phi(s_pos(2), u_rijk, a, g, Omega, kapper, Rd)

          factorz = exp( - N**2 / ( g * kapper ) * z_rcij )
          p_rijk = p_phi_rijk * factorz
          rho_rijk = ( p_phi_rijk / ( Rd * T_phi_rijk ) ) * factorz

          xyz_PressPerturb_k%val(1, i,j,rcID) = p_rijk - xyz_PressBasic_k%val(1, i,j,rcID)
          
           ! 厳密には, 
          xyz_RhoPerturb_k%val(1, i,j,rcID) = rho_rijk -  xyz_RhoBasic_k%val(1, i,j,rcID) 
				   ! だが, 上では密度の(スプライン)補間値を使うために, 以下の近似式の方が精度が良くなるので以下を使う. 
				   
				   !xyz_RhoPerturb_k%val(1, i,j,rcID) = xyz_PressPerturb_k%val(1, i,j,rcID) / ( Rd * T_phi_rijk ) 
				   
				   xyz_eNl_k%val(1, i,j,rcID) = SysParam%Cv * T_phi_rijk

				   xyz_vhNl_k%val(:, i,j,rcID) = geo_to_orth_vec((/ u_rijk, 0d0, 0d0 /), s_pos(:))

          if ( check_pole(icgrid2D, rcID, i, j) /= NOT_POLE_FLAG ) then
            xyz_vhNl_k%val(:, i,j,rcID) = 0d0
          end if

        end do
      end do
    end do
  end do

  call Field_IcGrid3D_Final(xyz_RhoBasic)
  call Field_IcGrid3D_Final(xyz_PressBasic)

end subroutine create_ZonalFlow_PressPerturbation_fields

function calc_u_phi(phi) result(u)
  real(DP), intent(in) :: phi
  real(DP) :: u

  u = u0 * cos(phi)

end function calc_u_phi

function calc_T_phi(phi, u, a, g, Omega, kapper, Rd) result(T_phi)
  real(DP), intent(in) :: phi
  real(DP), intent(in) :: u
  real(DP), intent(in) :: a
  real(DP), intent(in) :: g
  real(DP), intent(in) :: Omega
  real(DP), intent(in) :: kapper
  real(DP), intent(in) :: Rd
  real(DP) :: T_phi

  real(DP) :: tmp

  tmp = g - u**2 / a - 2d0 * Omega * u * cos(phi)
  T_phi = g * tmp * kapper / ( N**2 * Rd )

end function calc_T_phi

function calc_p_phi(phi, a, g, Omega, kapper) result(p_phi)
  real(DP), intent(in) :: phi
  real(DP), intent(in) :: a
  real(DP), intent(in) :: g
  real(DP), intent(in) :: Omega
  real(DP), intent(in) :: kapper
  real(DP) :: p_phi

  real(DP) :: tmp_G_phi, tmp_G_pole

  tmp_G_pole = calc_g1(0d0, a, g, Omega, kapper) / calc_g2(0d0, a, Omega)
  tmp_G_phi = calc_g1(phi, a, g, Omega, kapper) / calc_g2(phi, a, Omega)

  
  p_phi = P_pole * ( tmp_G_pole / tmp_G_phi )**( N**2 * a / ( 4d0 * g * kapper ) )

end function calc_p_phi

function calc_g1(phi, a, g, Omega, kapper) result(g1)
  real(DP), intent(in) :: phi
  real(DP), intent(in) :: a
  real(DP), intent(in) :: g
  real(DP), intent(in) :: Omega
  real(DP), intent(in) :: kapper
  real(DP) :: g1

  real(DP) :: factor1, factor2

  factor1 = ( 3d0 + 4d0 * cos(2d0*phi) + cos(4d0*phi) ) * u0**2
  factor2 = ( 1d0 + cos(2d0*phi) ) * u0 * a * g

  g1 =   2d0 * factor1 * ( u0**2 + 4d0 * a * Omega * (  u0 + a * Omega ) ) &
    &  - 16d0 * factor2 * ( u0 + 2d0 * a * Omega ) + 16d0 * a**2 * g**2

end function calc_g1

function calc_g2(phi, a, Omega) result(g2)
  real(DP), intent(in) :: phi
  real(DP), intent(in) :: a
  real(DP), intent(in) :: Omega
  real(DP) :: g2

  g2 = u0**2 * ( u0**2 + 4d0 * a * Omega * ( u0 + a * Omega ) )

end function calc_g2

end module experiment_impl_template
