!> \~japanese
!! \brief
!! \english
!! \brief
!!
!!
module small_time_step_mod

  ! モジュール引用 ; Use statements
  !

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

  ! 種類型パラメータ
  ! Kind type parameters
  !
  use dc_types, only: &
    & DP

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

use igmcore_math

  ! 座標変換
  ! Coordinate conversion
  !
  use igmcore_coordinate_conversion, only: &
    & orth_to_geo_pos, geo_to_orth_vec

  ! 正二十面体格子データの管理
  ! Manager of icosahedral grid data
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_EffSize_Min, get_EffSize_Max, get_IcRadius, &
    & 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

  ! 3 次元正二十面体格子データの管理
  ! Manager of 3D icosahedral grid data
  !
  use IcGrid3D_FVM_Manager, only: &
    & IcGrid3D_FVM, &
    & get_vertical_level_num, get_icgrid2D

  !
  !
  !
  use IcosahedralLorentzGrid_Manager, only: &
    & IcosahedralLorentzGrid, &
    & get_TerrainFollowingCoordinate, &
    & LORENTZGRID_HARF_LEVEL, LORENTZGRID_FULL_LEVEL
  
  ! 2 次元正二十面格子上の物理場データの管理
  ! Manager of the physical field on 2D icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & paste_field2D_margin

  ! 3 次元正二十面格子上の物理場データの管理
  ! Manager of the physical field on 3D icosahedral grid
  !
  use Field_IcGrid3D_Manager, only: &
    & Field_IcGrid3D, &
    & Field_IcGrid3D_Init, Field_IcGrid3D_Final, &
    & get_field_layer, get_icgrid3D, &
    & fill_value

  ! 正二十面格子と伴に地形に沿う座標系を利用するための補助モジュール
  ! Helper module for the use of terrain following coordinate with icosahedral grid.
  !
  use TFCoord_IcGrid3D_helper, only: &
    & TFCoord_IcGrid3D, &
    & TF_to_Z_IcGrid3D, calc_gamma, global_mean_Field_TFIcGrid3Dz

  ! 正二十面格子上の物理場データに対する水平空間微分演算の提供
  ! Providing horizontal spatial diffrential operations for data of physical field
  ! on icosahedral grid
  !
  use Derivate_Field_IcGrid2D_Manager, only: &
    & Derivate_Field_IcGrid2D

  ! 3 次元正二十面格子上の物理場データに対する鉛直微分演算の提供
  ! Provide vertical differential operations for data of physical field on 3D icosahedral grid.
  !
  use VDerivate_Field_IcGrid3D_Manager, only: &
    & VDerivate_Field_IcGrid3D, d_dVCoord

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

  ! 物理場データの鉛直補間
  ! Vertical interpolation of physical field data.
  !
  use Vertical_Interpolation, only: &
    & interpolate_fieldr_to_z, interpolate_fieldz_to_r

  !
  !
  !
  use TF_HVDerivateTool, only: &
    & TF_HGradient, TF_HDivergence, &
    & TF_VDerivate, TF_VDerivateVFlux

  use BoundaryCondition_Configurator, only: &
    & apply_UpperLowerBoundary_SlipCondition

  ! 球面座標系上でのデータ操作のためのユティリティ
  !
  !
  use Utility_GeoCoord, only: &
    & local_k

  ! 鉛直運動量を予測するための一次元ヘルムホルツ方程式のソルバ
  ! Solver for one-dimensional helmholtz equation to predict future vertical momentum.
  !
  use helmholtzEq_ForWDevAs_solver_mod, only: &
    & solve_helmholtzEq_forWDevAs
  
  ! 宣言文 ; Declaration statements
  !
  implicit none
  private
  
  ! 公開手続き
  ! Public procedures
  !
  public :: calc_VhDevAs, calc_WDevAs, calc_RDevAs, calc_EDevAs
  public :: calc_SR, calc_SW, calc_SP
  public :: calc_smallTStepDevVar, calc_PDev_From_EDev

contains


subroutine calc_PDev_From_EDev( &
  & xyz_PDev,       &  ! (inout)
  & xyz_EDev,       &  ! (in)
  & GasConst, Cv    &  ! (in)
  & )
  
  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: xyz_PDev
  type(Field_IcGrid3D), intent(in) :: xyz_EDev
  real(DP), intent(in) :: GasConst
  real(DP), intent(in) :: Cv

  ! 作業変数
  ! Local variables
  !
  type(Field_IcGrid2D), pointer :: xyz_PDev_k
  type(Field_IcGrid2D), pointer :: xyz_EDev_k

  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  integer :: k, kzMax

  ! 実行文 ; Executable statements
  !

  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 
  
  icgrid3Dz => get_icgrid3D(xyz_EDev)
  kzMax = get_vertical_level_num(icgrid3Dz)

  do k=1, kzMax
    xyz_PDev_k => get_field_layer(xyz_PDev, k)
    xyz_EDev_k => get_field_layer(xyz_EDev, k)

    xyz_PDev_k%val(:,:,:,:) = ( GasConst / Cv ) * xyz_EDev_k%val(:,:,:,:)

  end do

end subroutine calc_PDev_From_EDev

subroutine calc_VhDevAs( &
  & xyz_VhDevAs,                                   &  ! (inout)
  & xyz_VhDevNs, xyr_WDevNs, xyz_PDevNs,           &  ! (in)
  & xyz_RhoNl, xyz_VhNl, xyr_WNl, xyz_HGradPNl, xyz_GVhNl,    &  ! (in)
  & tfcoordz, derivate, vzderivate, &
  & DelTimeShort, alpha_dH          &
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: xyz_VhDevAs
  type(Field_IcGrid3D), intent(in) :: xyz_VhDevNs
  type(Field_IcGrid3D), intent(in) :: xyr_WDevNs
  type(Field_IcGrid3D), intent(in) :: xyz_PDevNs
  type(Field_IcGrid3D), intent(in) :: xyz_RhoNl    
  type(Field_IcGrid3D), intent(in) :: xyz_VhNl  
  type(Field_IcGrid3D), intent(in) :: xyr_WNl
  type(Field_IcGrid3D), intent(in) :: xyz_HGradPNl
  type(Field_IcGrid3D), intent(in) :: xyz_GVhNl
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(in) :: vzderivate
  real(DP), intent(in) :: DelTimeShort
  real(DP), intent(in) :: alpha_dH


  ! 作業変数
  ! Local variables
  !
  type(Field_IcGrid3D) :: xyz_VhNs, xyr_WNs
  type(Field_IcGrid3D) :: xyz_TFHDivVhNs
  type(Field_IcGrid3D) :: xyz_TFVDelWNs
  type(Field_IcGrid3D) :: xyz_PDevModNs
  type(Field_IcGrid3D) :: xyz_TFHGradPDevNs
  type(Field_IcGrid3D) :: xyz_UniCoef

  type(Field_IcGrid2D), pointer :: xyz_VhDevAs_k, xyz_VhDevNs_k
  type(Field_IcGrid2D), pointer :: xyz_VhNs_k, xyz_WNs_k
  type(Field_IcGrid2D), pointer :: xyz_TFHDivVhNs_k, xyz_TFVDelWNs_k
  type(Field_IcGrid2D), pointer :: PNl_k, xyz_PDevNs_k, xyz_PDevModNs_k, xyz_TFHGradPNl_k
  type(Field_IcGrid2D), pointer :: xyz_TFHGradPDevNs_k
  type(Field_IcGrid2D), pointer :: xyz_GVhNl_k

  integer :: rcID, i,j,k
  type(LocalProcRectInfo), pointer :: lpRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, KzMax, rcIDMin, rcIDMax

  type(IcGrid3D_FVM), pointer :: icgrid3Dz, icgrid3Dr
  type(IcGrid2D_FVM), pointer :: icgrid2D
  real(DP) :: s_pos(3), xi, dumpcoef, u_rijk(3)
  real(DP) :: univec_k(3)

  ! 実行文 ; Executable statements
  !

  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 
  
  icgrid3Dz => get_icgrid3D(xyz_VhDevNs)
  icgrid3Dr => get_icgrid3D(xyr_WDevNs)
  icgrid2D => get_icgrid2D(icgrid3Dz)

  lpRCInfo => get_LocalProcRectInfo(icgrid2D)  
  call get_InteriorRegion_IndexInfo( &
    & lpRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out) 

  KzMax = get_vertical_level_num(icgrid3Dz)

	! 3 次元正二十面格子上の物理場データを保持する作業変数を初期化する. 
	! Initialize all local variables which contain the data of field on 3D icosahedral grid. 
	!
  
  call Field_IcGrid3D_Init(xyz_VhNs, icgrid3Dz, 'xyz_VhNs', 3)
  call Field_IcGrid3D_Init(xyr_WNs, icgrid3Dr, 'xyr_WNs', 1)
  call Field_IcGrid3D_Init(xyz_TFHDivVhNs, icgrid3Dz, 'xyz_TFHDivVhDevNs', 1)
  call Field_IcGrid3D_Init(xyz_TFVDelWNs, icgrid3Dz, 'xyz_VDerivateWNs', 1)
  call Field_IcGrid3D_Init(xyz_PDevModNs, icgrid3Dz, 'xyz_PDevModNs', 1)
  call Field_IcGrid3D_Init(xyz_TFHGradPDevNs, icgrid3Dz, 'HGrad_PDevNs', 3)
  call Field_IcGrid3D_Init(xyz_UniCoef, icgrid3Dz, 'xyz_UniCoef', 1)

	call fill_value(xyz_UniCoef, 1.0d0)

  ! flux division method によって時間的に分割されている変数を結合して, もとに戻す.
  !
  ! phi(t) + phi*(m*tau) => phi(t+m*tau)
	call combine_LTSVarAndSTSDevVar(xyz_VhNs, xyz_VhNl, xyz_VhDevNs)
	call combine_LTSVarAndSTSDevVar(xyr_WNs, xyr_WNl, xyr_WDevNs)
	
  call TF_HDivergence( &
    & xyz_VhNs, xyz_TFHDivVhNs, &
    & tfcoordz, derivate, vzderivate )

  call TF_VDerivateVFlux( &
    & xyr_WNs, xyz_UniCoef, &
    & xyz_TFVDelWNs, tfcoordz)

  do k=1, kzMax
    xyz_TFHDivVhNs_k => get_field_layer(xyz_TFHDivVhNs, k)
    xyz_TFVDelWNs_k => get_field_layer(xyz_TFVDelWNs, k)
    xyz_PDevNs_k => get_field_layer(xyz_PDevNs, k)
    xyz_PDevModNs_k => get_field_layer(xyz_PDevModNs, k)

    xyz_PDevModNs_k%val(1, :,:,:) =   xyz_PDevNs_k%val(1, :,:,:) &
      &                            - alpha_dH * ( xyz_TFHDivVhNs_k%val(1, :,:,:) + xyz_TFVDelWNs_k%val(1, :,:,:) )
  
  end do

  call TF_HGradient( &
    & xyz_PDevModNs, xyz_TFHGradPDevNs, &
    & tfcoordz, derivate, vzderivate )

  !$omp parallel do private(xyz_TFHGradPDevNs_k, xyz_TFHGradPNl_k, AtildeNl_h_k, CtildeNl_h_k, xyz_VhDevNs_k, xyz_VhDevAs_k)
  do k=1, KzMax

    xyz_TFHGradPNl_k => get_field_layer(xyz_HGradPNl, k)
    xyz_TFHGradPDevNs_k => get_field_layer(xyz_TFHGradPDevNs, k)
    xyz_GVhNl_k => get_field_layer(xyz_GVhNl, k)
    xyz_VhDevNs_k => get_field_layer(xyz_VhDevNs, k)
    xyz_VhDevAs_k => get_field_layer(xyz_VhDevAs, k)

    xyz_VhDevAs_k%val(:,:,:,:) =   xyz_VhDevNs_k%val(:,:,:,:)                &
      &                          + DelTimeShort * (                          &
      &                            - xyz_TFHGradPDevNs_k%val(:,:,:,:)        &
      &                            - xyz_TFHGradPNl_k%val(:,:,:,:)           &
      &                            + xyz_GVhNl_k%val(:,:,:,:)                &
      &                               )



!xi = icgrid3Dz%VerticalGrid(k)
!if ( xi > 25d03 ) then
!DumpCoef = 0.5d0 * (1d0-cos(PI *(xi-25d03)/(40d03-25d03))) / 1800d0
!else
!DumpCoef = 0d0
!end if
!
!do rcID=rcIDMin,rcIDMax
!do j=jEMin, jEMax
!do i=iEMin, iEMax
!s_pos(:) = orth_to_geo_pos(icgrid2D%rcs_AGrid(:,i,j,rcID))
!
!if (abs(s_pos(2)) == PI/2d0) then
!u_rijk(:) = 0d0
!else
!u_rijk(:)=geo_to_orth_vec((/ 40d0*cos(s_pos(2)),0d0,0d0 /), s_pos(:))
!end if
!
!xyz_VhDevAs_k%val(:,i,j,rcID) = xyz_VhDevAs_k%val(:,i,j,rcID) - DumpCoef*DelTimeShort* &
!  & (  xyz_VhNs%field2Ds(k)%val(:,i,j,rcID) &
!  &  - xyz_RhoNl%field2Ds(k)%val(1,i,j,rcID) * u_rijk(:) )
!
!end do
!end do
!end do
    
  end do

	! 3 次元正二十面格子上の物理場データを保持するために用意した作業変数を最終化する. 
	! Finalize local variables prepared for storing data of field on 3D icosahedral grid. 
	!
	
  call Field_IcGrid3D_Final(xyz_VhNs)
  call Field_IcGrid3D_Final(xyr_WNs)
  call Field_IcGrid3D_Final(xyz_TFHGradPDevNs)
  call Field_IcGrid3D_Final(xyz_TFHDivVhNs)
  call Field_IcGrid3D_Final(xyz_TFVDelWNs)
  call Field_IcGrid3D_Final(xyz_PDevModNs)
  call Field_IcGrid3D_Final(xyz_UniCoef)

end subroutine calc_VhDevAs

subroutine calc_SR( &
  & SR, &
  & VhDevAs, VhNl, WNl, xyz_GRNl, &
  & tfcoordz, derivate, vzderivate &
  & )
  
  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: SR
  type(Field_IcGrid3D), intent(in) :: VhDevAs
  type(Field_IcGrid3D), intent(in) :: VhNl
  type(Field_IcGrid3D), intent(in) :: WNl
  type(Field_IcGrid3D), intent(in) :: xyz_GRNl
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(in) :: vzderivate

  ! 作業変数  
  ! Local variables
  !
  type(Field_IcGrid3D) :: tmp_TFHDiv_VhDevAsplusVhNl
  type(Field_IcGrid3D) :: TFHDiv_VhDevAsplusVhNl
  type(Field_IcGrid3D) :: xyz_UniCoef
  type(Field_IcGrid3D) :: xyz_TFVDelWNl

  type(Field_IcGrid2D), pointer :: VhDevAs_k
  type(Field_IcGrid2D), pointer :: VhNl_k
  type(Field_IcGrid2D), pointer :: tmp_TFHDiv_VhDevAsplusVhNl_k
  type(Field_IcGrid2D), pointer :: TFHDiv_VhDevAsplusVhNl_k
  type(Field_IcGrid2D), pointer :: TFVdel_WzNl_k
  type(Field_IcGrid2D), pointer :: SR_k, xyz_GRNl_k

  integer :: rcID, i,j,k
  type(LocalProcRectInfo), pointer :: lpRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, KzMax, rcIDMin, rcIDMax

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

  ! 実行文 ; Executable statements
  !
  
  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 

  icgrid3Dz => get_icgrid3D(VhNl)
  icgrid2D => get_icgrid2D(icgrid3Dz)

  lpRCInfo => get_LocalProcRectInfo(icgrid2D)  
  call get_InteriorRegion_IndexInfo( &
    & lpRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out) 

  kzMax = get_vertical_level_num(icgrid3Dz)

	! 3 次元正二十面格子上の物理場データを保持する作業変数を初期化する. 
	! Initialize all Local variables which contain the data of field on 3D icosahedral grid. 
	!

  call Field_IcGrid3D_Init(tmp_TFHDiv_VhDevAsplusVhNl, icgrid3Dz, 'VhDevAs+VhNl', 3)
  call Field_IcGrid3D_Init(TFHDiv_VhDevAsplusVhNl, icgrid3Dz, 'HDiv(VhDevAs+VhNl)', 1)
  call Field_IcGrid3D_Init(xyz_TFVDelWNl, icgrid3Dz, 'TFVdel_WzNl', 1)
  call Field_IcGrid3D_Init(xyz_UniCoef, icgrid3Dz, 'xyz_UniConf', 1)

  call fill_value(xyz_UniCoef, 1.0d0)
  
  call combine_LTSVarAndSTSDevVar(tmp_TFHDiv_VhDevAsplusVhNl, VhNl, VhDevAs)

  call TF_HDivergence( &
    & tmp_TFHDiv_VhDevAsplusVhNl, TFHDiv_VhDevAsplusVhNl, &
    & tfcoordz, derivate, vzderivate )

  call TF_VDerivateVFlux( &
    & WNl, xyz_UniCoef, &
    & xyz_TFVDelWNl, tfcoordz)


  !$omp parallel do private(TFHDiv_VhDevAsplusVhNl_k, TFVdel_WzNl_k, SR_k)
  do k=1, KzMax

    TFHDiv_VhDevAsplusVhNl_k => get_field_layer(TFHDiv_VhDevAsplusVhNl, k)
    TFVdel_WzNl_k => get_field_layer(xyz_TFVDelWNl, k)
    SR_k => get_field_layer(SR, k)
    xyz_GRNl_k => get_field_layer(xyz_GRNl, k)

    SR_k%val(1, :,:,:) = - ( TFHDiv_VhDevAsplusVhNl_k%val(1, :,:,:) + TFVdel_WzNl_k%val(1, :,:,:) ) &
      &                  + xyz_GRNl_k%val(1,:,:,:)                                                 

  end do

	! 3 次元正二十面格子上の物理場データを保持するために用意した作業変数を最終化する. 
	! Finalize local variables prepared for storing data of field on 3D icosahedral grid. 
	!
	
  call Field_IcGrid3D_Final(tmp_TFHDiv_VhDevAsplusVhNl)
  call Field_IcGrid3D_Final(TFHDiv_VhDevAsplusVhNl)
  call Field_IcGrid3D_Final(xyz_TFVDelWNl)
  call Field_IcGrid3D_Final(xyz_UniCoef)

end subroutine calc_SR

subroutine calc_SW( &
  & SW, &
  & buoyancyNl, GWNl, &
  & tfcoordz  &
  &  )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: SW
  type(Field_IcGrid3D), intent(in) :: buoyancyNl
  type(Field_IcGrid3D), intent(in) :: GWNl
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz

  ! 作業変数
  ! Local variables
  !
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(IcGrid3D_FVM), pointer :: icgrid3Dr

  type(Field_IcGrid3D) :: xyr_buoyancyNl

  type(Field_IcGrid2D), pointer :: SW_k, xyr_buoyancyNl_k, GWNl_k

  real(DP) :: s_pos(3)
  real(DP) :: univec_k(3)

  integer :: rcID, i, j, k
  type(LocalProcRectInfo), pointer :: lpRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, KrMax, rcIDMin, rcIDMax

  ! 実行文 ; Executable statements
  !

  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 

  icgrid3Dr => get_icgrid3D(SW)
  icgrid2D => get_icgrid2D(icgrid3Dr)

  lpRCInfo => get_LocalProcRectInfo(icgrid2D)  
  call get_InteriorRegion_IndexInfo( &
    & lpRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out) 

  KrMax = get_vertical_level_num(icgrid3Dr)

	! 3 次元正二十面格子上の物理場データを保持する作業変数を初期化する. 
	! Initialize all Local variables which contain the data of field on 3D icosahedral grid. 
	!

  call Field_IcGrid3D_Init(xyr_buoyancyNl, icgrid3Dr,  'xyr_bouyancyNl', 1)
  call interpolate_fieldz_to_r(buoyancyNl, xyr_buoyancyNl)
  
  do k=2, KrMax
    SW_k => get_field_layer(SW, k)
    xyr_buoyancyNl_k => get_field_layer(xyr_buoyancyNl, k)
    GWNl_k => get_field_layer(GWNl, k)    
  
    SW_k%val(1,:,:,:) = xyr_buoyancyNl_k%val(1,:,:,:) + GWNl_k%val(1,:,:,:)
  end do

  SW%field2Ds(1)%val(1,:,:,:) = 0d0
  
	! 3 次元正二十面格子上の物理場データを保持するために用意した作業変数を最終化する. 
	! Finalize local variables prepared for storing data of field on 3D icosahedral grid. 
	!
	
  call Field_IcGrid3D_Final(xyr_buoyancyNl)

end subroutine calc_SW

subroutine calc_SP( &
  & SP, &
  & VhDevAs, VhNl, WNl, &
  & xyz_RhoNl, xyz_HGradPNl, xyz_hNl, xyz_GEaNl, &
  & gtildeNl, &
  & tfcoordz, derivate, vzderivate, &
  & Rd, Cp, Cv &
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: SP
  type(Field_IcGrid3D), intent(in) :: VhDevAs
  type(Field_IcGrid3D), intent(in) :: VhNl
  type(Field_IcGrid3D), intent(in) :: WNl
  type(Field_IcGrid3D), intent(in) :: xyz_RhoNl
  type(Field_IcGrid3D), intent(in) :: xyz_HGradPNl
  type(Field_IcGrid3D), intent(in) :: xyz_hNl
  type(Field_IcGrid3D), intent(in) :: xyz_GEaNl
  type(Field_IcGrid3D), intent(in) :: gtildeNl
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(in) :: vzderivate
  real(DP), intent(in) :: Rd
  real(DP), intent(in) :: Cp
  real(DP), intent(in) :: Cv

  ! 作業変数
  ! Local variables
  !
  type(Field_IcGrid3D) :: xyz_WNl
  type(Field_IcGrid3D) :: EnthalpyHFluxNl
  type(Field_IcGrid3D) :: HDiv_EnthalpyHFluxNl
  type(Field_IcGrid3D) :: VDel_EnthalpyVFluxNl

  type(Field_IcGrid2D), pointer :: SP_k
  type(Field_IcGrid2D), pointer :: xyz_RhoNl_k
  type(Field_IcGrid2D), pointer :: xyz_hNl_k
  type(Field_IcGrid2D), pointer :: VhDevAs_k
  type(Field_IcGrid2D), pointer :: VhNl_k
  type(Field_IcGrid2D), pointer :: xyz_WNl_k
  type(Field_IcGrid2D), pointer :: xyz_GEaNl_k
  type(Field_IcGrid2D), pointer :: gtildeNl_k

  type(Field_IcGrid2D), pointer :: xyz_HGradPNl_k
  type(Field_IcGrid2D), pointer :: EnthalpyHFlux_AsPlusNl_k
  type(Field_IcGrid2D), pointer :: HDiv_EnthalpyHFlux_AsPlusNl_k
  type(Field_IcGrid2D), pointer :: VDel_EnthalpyVFluxNl_k

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

  real(DP) :: h_rijk

  integer :: rcID, i, j, k, rankID
  type(LocalProcRectInfo), pointer :: lpRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, KzMax, rcIDMin, rcIDMax

  ! 実行文 ; Executable statements
  !

  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 

  icgrid3Dz => get_icgrid3D(VhNl)
  icgrid2D => get_icgrid2D(icgrid3Dz)

  lpRCInfo => get_LocalProcRectInfo(icgrid2D)  
  call get_InteriorRegion_IndexInfo( &
    & lpRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out) 

  KzMax = get_vertical_level_num(icgrid3Dz)

	! 3 次元正二十面格子上の物理場データを保持する作業変数を初期化する. 
	! Initialize all Local variables which contain the data of field on 3D icosahedral grid. 
	!

  call Field_IcGrid3D_Init(xyz_WNl, icgrid3Dz, 'xyz_WzNl', 1)
  call Field_IcGrid3D_Init(EnthalpyHFluxNl, icgrid3Dz, 'EnthalpyHFlux_AsPlusNl', 3)
  call Field_IcGrid3D_Init(HDiv_EnthalpyHFluxNl, icgrid3Dz, 'HDiv_EnthalpyHFlux_AsPlusNl', 1)
  call Field_IcGrid3D_Init(VDel_EnthalpyVFluxNl, icgrid3Dz, 'VDel_EnthalpyVFluxNl', 1)

  call interpolate_fieldr_to_z(WNl, xyz_WNl)

  !$omp parallel do private(VhDevAs_k, VhNl_k, xyz_hNl_k, xyz_WNl_k, EnthalpyHFlux_AsPlusNl_k)
  do k=1, KzMax

    VhDevAs_k => get_field_layer(VhDevAs, k)
    VhNl_k => get_field_layer(VhNl, k)
    xyz_hNl_k => get_field_layer(xyz_hNl, k)
    xyz_WNl_k => get_field_layer(xyz_WNl, k)
    EnthalpyHFlux_AsPlusNl_k => get_field_layer(EnthalpyHFluxNl, k)
    
    do rankID=1, 3
      EnthalpyHFlux_AsPlusNl_k%val(rankID, :,:,:) = &
            & xyz_hNl_k%val(1, :,:,:) * ( &
            &  VhDevAs_k%val(rankID, :,:,:) + VhNl_k%val(rankID, :,:,:) &
            & )
     end do
  end do

  call TF_HDivergence( &
    & EnthalpyHFluxNl, HDiv_EnthalpyHFluxNl, &
    & tfcoordz, derivate, vzderivate )

  call TF_VDerivateVFlux( &
    & WNl, xyz_hNl, &
    & VDel_EnthalpyVFluxNl, tfcoordz)

  do k=1, KzMax
    SP_k => get_field_layer(SP, k)
    xyz_RhoNl_k => get_field_layer(xyz_RhoNl, k)
    VhNl_k => get_field_layer(VhNl, k)
    xyz_WNl_k => get_field_layer(xyz_WNl, k)
    xyz_GEaNl_k => get_field_layer(xyz_GEaNl, k)

    gtildeNl_k => get_field_layer(gtildeNl, k)
    HDiv_EnthalpyHFlux_AsPlusNl_k => get_field_layer(HDiv_EnthalpyHFluxNl, k)
    VDel_EnthalpyVFluxNl_k => get_field_layer(VDel_EnthalpyVFluxNl, k)
    xyz_HGradPNl_k => get_field_layer(xyz_HGradPNl, k)

    do rcID=rcIDMin, rcIDMax
      !$omp parallel do
      do j=jEMin, jEMax
        do i=iEMin, iEMax

        SP_k%val(1, i,j,rcID) = ( Rd / Cv ) * &
            & ( &
            &     xyz_GEaNl_k%val(1, i,j,rcID)        &
            &   - HDiv_EnthalpyHFlux_AsPlusNl_k%val(1, i,j,rcID)                   &
            &   - VDel_EnthalpyVFluxNl_k%val(1, i,j,rcID)                   &
            &   + dot_product(VhNl_k%val(1:3, i,j,rcID), xyz_HGradPNl_k%val(1:3, i,j,rcID))          &
            &     / xyz_RhoNl_k%val(1, i,j,rcID) &
            &   - gtildeNl_k%val(1, i,j,rcID) * xyz_WNl_k%val(1, i,j,rcID)          &
            & )

        end do
      end do
    end do
  end do

	! 3 次元正二十面格子上の物理場データを保持するために用意した作業変数を最終化する. 
	! Finalize local variables prepared for storing data of field on 3D icosahedral grid. 
	!
	
  call Field_IcGrid3D_Final(xyz_WNl)
  call Field_IcGrid3D_Final(EnthalpyHFluxNl)
  call Field_IcGrid3D_Final(HDiv_EnthalpyHFluxNl)
  call Field_IcGrid3D_Final(VDel_EnthalpyVFluxNl)

end subroutine calc_SP

subroutine calc_WDevAs( &
  & xyr_WDevAs, xyr_WDevNs, &
  & xyz_RDevNs, xyz_PDevNs, xyz_VhDevAs, &
  & xyz_VhNl, xyr_WNl, xyz_RhoNl, xyz_hNl, xyz_gtildeNl, &
  & xyz_SR, xyz_SP, xyr_SW, &
  & IcosaLorentzGrid, tfcoordz, tfcoord, vrderivate, &
  & Grav, Rd, Cp, Cv, DelTimeShort, &
  & hydrostatic_opt &
  & )

use igmcore_math

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: xyr_WDevAs
  type(Field_IcGrid3D), intent(in) :: xyr_WDevNs
  type(Field_IcGrid3D), intent(in) :: xyz_RDevNs
  type(Field_IcGrid3D), intent(in) :: xyz_PDevNs
  type(Field_IcGrid3D), intent(in) :: xyz_VhDevAs
  type(Field_IcGrid3D), intent(in) :: xyz_VhNl
  type(Field_IcGrid3D), intent(in) :: xyr_WNl
  type(Field_IcGrid3D), intent(in) :: xyz_RhoNl
  type(Field_IcGrid3D), intent(in) :: xyz_hNl
  type(Field_IcGrid3D), intent(in) :: xyz_gtildeNl
  type(Field_IcGrid3D), intent(in) :: xyz_SR
  type(Field_IcGrid3D), intent(in) :: xyz_SP
  type(Field_IcGrid3D), intent(in) :: xyr_SW
  type(IcosahedralLorentzGrid), intent(in) :: icosaLorentzGrid
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  type(TFCoord_IcGrid3D), intent(in) :: tfcoord
  type(VDerivate_Field_IcGrid3D), intent(in) :: vrderivate
  real(DP), intent(in) :: Grav
  real(DP), intent(in) :: Rd
  real(DP), intent(in) :: Cp
  real(DP), intent(in) :: Cv
  real(DP), intent(in) :: DelTimeShort
  logical, intent(in) :: hydrostatic_opt

  ! 作業変数
  ! Local variables
  !
  type(IcGrid3D_FVM), pointer :: icgrid3Dz, icgrid3Dr
  type(Field_IcGrid3D) :: xyz_VhAs
  integer :: k
  
  type(Field_IcGrid3D) :: xyr_RDevNs
  type(Field_IcGrid3D) :: xyr_PDevNs
  type(Field_IcGrid3D) :: xyr_hNl
  type(Field_IcGrid3D) :: xyr_gtildeNl
  type(Field_IcGrid3D) :: xyr_SR
  type(Field_IcGrid3D) :: xyr_SP

  integer :: rcID, i, j, KrMax
  type(LocalProcRectInfo), pointer :: lpRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, KzMax, rcIDMin, rcIDMax

real(DP) :: DumpCoef, xi

  ! 実行文 ; Executable statements
  !

  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 
  
  icgrid3Dz => get_icgrid3D(xyz_VhNl)
  icgrid3Dr => get_icgrid3D(xyr_WDevNs)

  lpRCInfo => get_LocalProcRectInfo(get_icgrid2D(icgrid3Dz))  
  call get_InteriorRegion_IndexInfo( &
    & lpRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out) 
  
	! 3 次元正二十面格子上の物理場データを保持する作業変数を初期化する. 
	! Initialize all Local variables which contain the data of field on 3D icosahedral grid. 
	!

  call Field_IcGrid3D_Init(xyr_RDevNs, icgrid3Dr, 'xyr_RDevNs', 1)
  call Field_IcGrid3D_Init(xyr_PDevNs, icgrid3Dr, 'xyr_PDevNs', 1)
  call Field_IcGrid3D_Init(xyr_hNl, icgrid3Dr, 'xyr_hNl', 1)
  call Field_IcGrid3D_Init(xyr_gtildeNl, icgrid3Dr, 'xyr_gtildeNl', 1)
  call Field_IcGrid3D_Init(xyr_SR, icgrid3Dr, 'xyr_SR', 1)
  call Field_IcGrid3D_Init(xyr_SP, icgrid3Dr, 'xyr_SP', 1)

  call interpolate_fieldz_to_r(xyz_RDevNs, xyr_RDevNs)
  call interpolate_fieldz_to_r(xyz_PDevNs, xyr_PDevNs)
  call interpolate_fieldz_to_r(xyz_hNl, xyr_hNl)
  call interpolate_fieldz_to_r(xyz_gtildeNl, xyr_gtildeNl)
  call interpolate_fieldz_to_r(xyz_SR, xyr_SR)
  call interpolate_fieldz_to_r(xyz_SP, xyr_SP)

  ! 鉛直運動量 W に対する一次元ヘルムホルツ方程式を陰的に解いて,
  ! 未来の W を予測する.
  !

  
  !
  !
  !

  call solve_helmholtzEq_forWDevAs( &
    & xyr_WDevAs, xyr_WDevNs,              &
    & xyr_RDevNs, xyr_PDevNs, xyz_VhDevAs, &
    & xyr_hNl, xyr_gtildeNl, &
    & xyr_WNl,xyr_SR, xyr_SP, xyr_SW, &
    & IcosaLorentzGrid, tfcoord, vrderivate, &
    & Grav, Rd, Cp, Cv, DelTimeShort, &
    & hydrostatic_opt &
  & )


	! 3 次元正二十面格子上の物理場データを保持するために用意した作業変数を最終化する. 
	! Finalize local variables prepared for storing data of field on 3D icosahedral grid. 
	!
	
  call Field_IcGrid3D_Final(xyr_RDevNs)
  call Field_IcGrid3D_Final(xyr_PDevNs)
  call Field_IcGrid3D_Final(xyr_hNl)
  call Field_IcGrid3D_Final(xyr_gtildeNl)
  call Field_IcGrid3D_Final(xyr_SR)
  call Field_IcGrid3D_Final(xyr_SP)

end subroutine calc_WDevAs

subroutine calc_RDevAs( &
  & xyz_RDevAs, xyz_RDevNs, &
  & xyr_WDevAs, xyz_SR, &
  & tfcoordz, vzderivate, &
  & DelTimeShort &
  & )
  
  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: xyz_RDevAs
  type(Field_IcGrid3D), intent(in) :: xyz_RDevNs
  type(Field_IcGrid3D), intent(in) :: xyr_WDevAs
  type(Field_IcGrid3D), intent(in) :: xyz_SR
  type(TFCoord_IcGrid3D), intent(in)  :: tfcoordz
  type(VDerivate_Field_IcGrid3D), intent(in) :: vzderivate
  real(DP), intent(in) :: DelTimeShort

  ! 作業変数
  ! Local variables
  !
  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  type(IcGrid2D_FVM), pointer :: icgrid2D

  type(Field_IcGrid3D) :: xyz_WDevAs, xyz_UniCoef
  type(Field_IcGrid3D) :: xyz_TFVDelWDevAs

  type(Field_IcGrid2D), pointer :: xyz_RDevAs_k, xyz_RDevNs_k
  type(Field_IcGrid2D), pointer :: xyz_WDevAs_k
  type(Field_IcGrid2D), pointer :: VDel_xyz_WDevAs_k, xyz_SR_k
  integer :: kzMax
  integer :: k

  ! 実行文 ; Executable statements
  !

   ! 格子に関する情報を取得する.
  ! Get the information about grid.

  icgrid3Dz => get_icgrid3D(xyz_RDevAs)
  kzMax = get_vertical_level_num(icgrid3Dz)

	! 3 次元正二十面格子上の物理場データを保持する作業変数を初期化する. 
	! Initialize all Local variables which contain the data of field on 3D icosahedral grid. 
	!

  call Field_IcGrid3D_Init(xyz_UniCoef, icgrid3Dz, 'xyz_UniConf', 1)
  call Field_IcGrid3D_Init(xyz_TFVDelWDevAs, icgrid3Dz, 'xyz_TFVDelWDevAs', 1)

  call fill_value(xyz_UniCoef, 1.0d0)

  call TF_VDerivateVFlux( &
    & xyr_WDevAs, xyz_UniCoef, &
    & xyz_TFVDelWDevAs, tfcoordz)

  !$omp parallel do private(xyz_RDevAs_k, xyz_RDevNs_k, xyz_SR_k, VDel_xyz_WDevAs_k)
  do k=1, kzMax
    xyz_RDevAs_k => get_field_layer(xyz_RDevAs, k)
    xyz_RDevNs_k => get_field_layer(xyz_RDevNs, k)
    xyz_SR_k => get_field_layer(xyz_SR, k)
    VDel_xyz_WDevAs_k => get_field_layer(xyz_TFVDelWDevAs, k)

    xyz_RDevAs_k%val(1, :,:,:) = &
       &   xyz_RDevNs_k%val(1, :,:,:) &
       & + DelTimeShort * ( - VDel_xyz_WDevAs_k%val(1, :,:,:) + xyz_SR_k%val(1, :,:,:) )

  end do

	! 3 次元正二十面格子上の物理場データを保持するために用意した作業変数を最終化する. 
	! Finalize local variables prepared for storing data of field on 3D icosahedral grid. 
	!
	
  call Field_IcGrid3D_Final(xyz_UniCoef)  
  call Field_IcGrid3D_Final(xyz_TFVDelWDevAs)

end subroutine calc_RDevAs


subroutine calc_EDevAs( &
  & EDevAs, EDevNs, ENl, &
  & RDevAs, VhDevAs, WDevAs, &
  & RDevNs, VhDevNs, WDevNs, &
  & xyz_RhoNl, VhNl, WNl, xyz_hNl, &
  & xyz_GEaNl, TEAs, TENs, &
  & tfcoordz, derivate, vzderivate, &
  & Grav, Cp, Cv, DelTimeShort )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: EDevAs
  type(Field_IcGrid3D), intent(inout) :: EDevNs
  type(Field_IcGrid3D), intent(in) :: ENl
  type(Field_IcGrid3D), intent(in) :: RDevAs
  type(Field_IcGrid3D), intent(in) :: VhDevAs
  type(Field_IcGrid3D), intent(in) :: WDevAs
  type(Field_IcGrid3D), intent(in) :: RDevNs
  type(Field_IcGrid3D), intent(in) :: VhDevNs
  type(Field_IcGrid3D), intent(in) :: WDevNs
  type(Field_IcGrid3D), intent(in) :: xyz_RhoNl
  type(Field_IcGrid3D), intent(in) :: VhNl
  type(Field_IcGrid3D), intent(in) :: WNl
  type(Field_IcGrid3D), intent(in) :: xyz_hNl
  type(Field_IcGrid3D), intent(in) :: xyz_GEaNl
  type(Field_IcGrid3D), intent(in) :: TEAs
  type(Field_IcGrid3D), intent(in) :: TENs
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(in) :: vzderivate
  real(DP), intent(in) :: Grav
  real(DP), intent(in) :: Cp
  real(DP), intent(in) :: Cv
  real(DP), intent(in) :: DelTimeShort

  ! 作業変数
  ! Local variables
  !
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  type(Field_IcGrid3D) :: xyr_WAs, xyz_WNl, xyz_WDevNs, xyz_WDevAs
  type(Field_IcGrid3D) :: TEHFlux, HDiv_TEHFlux
  type(Field_IcGrid3D) :: TEVFlux, VDel_TEVFlux

  type(Field_IcGrid2D), pointer :: EDevNs_k, EDevAs_k, ENl_k
  type(Field_IcGrid2D), pointer :: xyz_hNl_k
  type(Field_IcGrid2D), pointer :: TEAs_k, TENs_k
  type(Field_IcGrid2D), pointer :: RDevAs_k, RDevNs_k, xyz_RhoNl_k
  type(Field_IcGrid2D), pointer :: VhNl_k, VhDevNs_k, VhDevAs_k
  type(Field_IcGrid2D), pointer :: WzNl_k, WDevzNs_k, WDevzAs_k
  type(Field_IcGrid2D), pointer :: xyz_GEaNl_k
  type(Field_IcGrid2D), pointer :: TEHFlux_k, HDiv_TEHFlux_k
  type(Field_IcGrid2D), pointer :: TEVFlux_k, VDel_TEVFlux_k

  real(DP) :: z_rijk, G_1div2_rijk
  real(DP) :: h_rijk, Phi_rijk, ke_rijk
  real(DP) :: tmp_TEFlux
  real(DP) :: Vh_rijk(3), W_rijk, Rho_rijk
  real(DP) :: s_pos(3)
  real(DP) :: univec_k(3)
  real(DP) :: gamma_
  real(DP) :: radius

  integer :: rcID, i, j, k, valRankID
  type(LocalProcRectInfo), pointer :: lpRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, KzMax, rcIDMin, rcIDMax

  ! 実行文 ; Executable statements
  !

  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 
  
  icgrid3Dz => get_icgrid3D(VhNl)
  icgrid2D => get_icgrid2D(icgrid3Dz)

  lpRCInfo => get_LocalProcRectInfo(icgrid2D)  
  call get_InteriorRegion_IndexInfo( &
    & lpRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out) 

  KzMax = get_vertical_level_num(icgrid3Dz)
  radius = get_IcRadius(icgrid2D)

	! 3 次元正二十面格子上の物理場データを保持する作業変数を初期化する. 
	! Initialize all Local variables which contain the data of field on 3D icosahedral grid. 
	!

  call Field_IcGrid3D_Init(xyz_WDevNs, icgrid3Dz, 'WDevzNs', 1)
  call Field_IcGrid3D_Init(xyz_WDevAs, icgrid3Dz, 'WDevzAs', 1)
  call Field_IcGrid3D_Init(xyr_WAs, get_icgrid3D(WDevNs), 'xyr_WAs', 1)
  call Field_IcGrid3D_Init(xyz_WNl, icgrid3Dz, 'WzNl', 1)
  call Field_IcGrid3D_Init(TEHFlux, icgrid3Dz, 'TEHFlux', 3)
  call Field_IcGrid3D_Init(TEVFlux, icgrid3Dz, 'TEVFlux', 1)
  call Field_IcGrid3D_Init(HDiv_TEHFlux, icgrid3Dz, 'HDiv_TEHFlux', 1)
  call Field_IcGrid3D_Init(VDel_TEVFlux, icgrid3Dz, 'VDel_TEVFlux', 1)

  call interpolate_fieldr_to_z(WDevNs, xyz_WDevNs)
  call interpolate_fieldr_to_z(WDevAs, xyz_WDevAs)
  call interpolate_fieldr_to_z(WNl, xyz_WNl)
    
  call combine_LTSVarAndSTSDevVar(xyr_WAs, WNl, WDevAs)
  
  do k=1, kzMax
    TEHFlux_k => get_field_layer(TEHFlux, k)
    TEVFlux_k => get_field_layer(TEVFlux, k)
    TENs_k => get_field_layer(TENs, k)

    ENl_k => get_field_layer(ENl, k)
    EDevNs_k => get_field_layer(EDevNs, k)

    xyz_RhoNl_k => get_field_layer(xyz_RhoNl, k)
    RDevNs_k => get_field_layer(RDevNs, k)

    VhNl_k => get_field_layer(VhNl, k)
    VhDevNs_k => get_field_layer(VhDevNs, k)
    VhDevAs_k => get_field_layer(VhDevAs, k)
    
    WzNl_k => get_field_layer(xyz_WNl, k)
    WDevzNs_k => get_field_layer(xyz_WDevNs, k)
    WDevzAs_k => get_field_layer(xyz_WDevAs, k)

    xyz_hNl_k => get_field_layer(xyz_hNl, k)
!
    do rcID=rcIDMin, rcIDMax
      !$omp parallel do private(i, z_rijk, Rho_rijk, Vh_rijk, W_rijk, Phi_rijk, ke_rijk, tmp_TEFlux)
      do j=jEMin, jEMax
        do i=iEMin, iEMax
 

          !
          !
          z_rijk = TF_to_Z_IcGrid3D(tfCoordz, icgrid3Dz%verticalGrid(k), rcID, i, j)

          Rho_rijk = xyz_RhoNl_k%val(1, i,j,rcID)
          Vh_rijk(1:3) = VhNl_k%val(1:3, i,j,rcID)
          W_rijk = WzNl_k%val(1, i,j,rcID)

          Phi_rijk = Grav * z_rijk
          ke_rijk = 0.5d0 * ( dot_product(Vh_rijk(:), Vh_rijk(:)) + W_rijk**2 ) / Rho_rijk**2

          tmp_TEFlux = xyz_hNl_k%val(1, i,j,rcID) + Phi_rijk + ke_rijk

          TEHFlux_k%val(1:3, i,j,rcID) = &
            & tmp_TEFlux * ( Vh_rijk(1:3) + VhDevAs_k%val(1:3, i,j,rcID) )

          TEVFlux_k%val(1, i,j,rcID) = tmp_TEFlux

          !
          !
          Rho_rijk = Rho_rijk + RDevNs_k%val(1, i,j,rcID)
          Vh_rijk(1:3) = Vh_rijk(1:3) + VhDevNs_k%val(1:3, i,j,rcID)
          W_rijk = W_rijk + WDevzNs_k%val(1, i,j,rcID)

          ke_rijk = 0.5d0 * ( dot_product(Vh_rijk(:), Vh_rijk(:)) + W_rijk**2 ) / Rho_rijk**2

          TENs_k%val(1, i,j,rcID) = &
            &    ( ENl_k%val(1, i,j,rcID) + EDevNs_k%val(1, i,j,rcID) ) &
            &  + Rho_rijk * ( Phi_rijk + ke_rijk )

        end do
      end do
    end do
  end do

  call TF_HDivergence( &
    & TEHFlux, HDiv_TEHFlux, &
    & tfcoordz, derivate, vzderivate )
    
  call TF_VDerivateVFlux( &
    & xyr_WAs, TEVFlux, &
    & VDel_TEVFlux, tfcoordz)

  !$omp parallel do private(TEAs_k, TENs_k, HDiv_TEHFlux_k, VDel_TEVFlux_k, xyz_GEaNl_k)
  do k=1, KzMax

    TEAs_k => get_field_layer(TEAs, k)
    TENs_k => get_field_layer(TENs, k)
    HDiv_TEHFlux_k => get_field_layer(HDiv_TEHFlux, k)
    VDel_TEVFlux_k => get_field_layer(VDel_TEVFlux, k)
    xyz_GEaNl_k => get_field_layer(xyz_GEaNl, k)

    TEAs_k%val(1, :,:,:) = TENs_k%val(1, :,:,:) +  &
                        & DelTimeShort * ( &
                        &   - HDiv_TEHFlux_k%val(1, :,:,:) &
                        &   - VDel_TEVFlux_k%val(1, :,:,:) &
                        &   + xyz_GEaNl_k%val(1, :,:,:)     &
                            &  )
  end do

  do k=1, KzMax

    TEAs_k => get_field_layer(TEAs, k)

    ENl_k => get_field_layer(ENl, k)
    EDevAs_k => get_field_layer(EDevAs, k)

    VhNl_k => get_field_layer(VhNl, k)
    VhDevAs_k => get_field_layer(VhDevAs, k)
    WzNl_k => get_field_layer(xyz_WNl, k)
    WDevzAs_k => get_field_layer(xyz_WDevAs, k)

    xyz_RhoNl_k => get_field_layer(xyz_RhoNl, k)
    RDevAs_k => get_field_layer(RDevAs, k)

    do rcID=rcIDMin, rcIDMax
      !$omp parallel do private(z_rijk, Rho_rijk, Vh_rijk, W_rijk, Phi_rijk, ke_rijk)
      do j=jEMin, jEMax
        do i=iEMin, iEMax

          z_rijk = TF_to_Z_IcGrid3D(tfCoordz, icgrid3Dz%verticalGrid(k), rcID, i, j)
          Rho_rijk = xyz_RhoNl_k%val(1, i,j,rcID) + RDevAs_k%val(1, i,j,rcID)
          Vh_rijk(1:3) = VhNl_k%val(1:3, i,j,rcID) + VhDevAs_k%val(1:3, i,j,rcID)
          W_rijk = WzNl_k%val(1, i,j,rcID) + WDevzAs_k%val(1, i,j,rcID)

          Phi_rijk = Grav * z_rijk
          ke_rijk = 0.5d0 * ( dot_product(Vh_rijk(:), Vh_rijk(:)) + W_rijk**2 ) / Rho_rijk**2

          EDevAs_k%val(1, i,j,rcID) = &
            &    TEAs_k%val(1, i,j,rcID)            &
            &  - Rho_rijk * ( Phi_rijk + ke_rijk ) &
            &  - ENl_k%val(1, i,j,rcID)
          
        end do
      end do
    end do

  end do
  
  ! 3 次元正二十面格子上の物理場データを保持するために用意した作業変数を最終化する. 
	! Finalize local variables prepared for storing data of field on 3D icosahedral grid. 
	!
  
  call Field_IcGrid3D_Final(xyz_WDevNs)
  call Field_IcGrid3D_Final(xyz_WDevAs)
  call Field_IcGrid3D_Final(xyr_WAs)
  call Field_IcGrid3D_Final(xyz_WNl)
  call Field_IcGrid3D_Final(TEHFlux)
  call Field_IcGrid3D_Final(TEVFlux)
  call Field_IcGrid3D_Final(HDiv_TEHFlux)
  call Field_IcGrid3D_Final(VDel_TEVFlux)

end subroutine calc_EDevAs

subroutine combine_LTSVarAndSTSDevVar(combinedVar, currentLTSVar, devVar)
	type(Field_IcGrid3D), intent(inout) :: combinedVar
	type(Field_IcGrid3D), intent(in) :: devVar
	type(Field_IcGrid3D), intent(in) :: currentLTSVar
  
  ! 作業変数
  ! Local variables
  !
  integer :: k
  integer :: vtLvNum
  type(IcGrid3D_FVM), pointer :: icgrid3D
  type(Field_IcGrid2D), pointer :: combinedVar_k,devVar_k, currentLTSVar_k

  ! 実行文 ; Executable statements
  !

  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 
  
  icgrid3D => get_icgrid3D(devVar)
  vtLvNum = get_vertical_level_num(icgrid3D)

	!$omp parallel do private(combinedVar_k, devVar_k, currentLTSVar_k)
  do k=1, vtLvNum
    combinedVar_k => get_field_layer(combinedVar, k)
    devVar_k => get_field_layer(devVar, k)
    currentLTSVar_k => get_field_layer(currentLTSVar, k)
    
    combinedVar_k%val(:,:,:,:) = devVar_k%val(:,:,:,:) + currentLTSVar_k%val(:,:,:,:)
    
  end do
  		
end subroutine combine_LTSVarAndSTSDevVar

subroutine calc_smallTStepDevVar(devVar, var, currentLTVar)

  ! 宣言文 ; Declaration statements
  !

  type(Field_IcGrid3D), intent(inout) :: devVar
  type(Field_IcGrid3D), intent(in) :: var
  type(Field_IcGrid3D), intent(in) :: currentLTVar
  
  ! 作業変数
  ! Local variables
  !
  integer :: k
  integer :: vtLvNum
  type(IcGrid3D_FVM), pointer :: icgrid3D
  type(Field_IcGrid2D), pointer :: devVar_k, var_k, currentLTVar_k

  ! 実行文 ; Executable statements
  !

  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 
  
  icgrid3D => get_icgrid3D(var)
  vtLvNum = get_vertical_level_num(icgrid3D)

  do k=1, vtLvNum
    devVar_k => get_field_layer(devVar, k)
    var_k => get_field_layer(var, k)
    currentLTVar_k => get_field_layer(currentLTVar, k)
    
    devVar_k%val(:,:,:,:) = var_k%val(:,:,:,:) - currentLTVar_k%val(:,:,:,:)
    
  end do

end subroutine calc_smallTStepDevVar

end module small_time_step_mod
