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

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

#ifdef LIB_MPI
  use mpi
#endif

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

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

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

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


  ! 正二十面体格子データの管理
  ! 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

  ! 正二十面格子と伴に地形に沿う座標系を利用するための補助モジュール
  ! 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

  ! 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_icgrid3D, get_field3d_name, get_field_layer, &
    & copy_field3D

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

  ! 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 Utility_GeoCoord, only: &
    & local_k

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

  !
  !
  !  
  use numerical_filter_mod, only: &
    & numFilter_HVDiffusion
  
  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public statements
  !
  public :: calc_RhoNl, calc_hNl
  public :: calc_AtildeNl, calc_CtildeNl, calc_gtildeNl
  public :: calc_buoyancyNl, calc_HGradPNl
  public :: calc_GR_GVh_GW
  public :: apply_numerical_filter
  public :: update_variable

contains

subroutine calc_RhoNl( &
  & xyz_RhoNl,             &  ! (inout)
  & xyz_RhoBasic, xyz_RNl  &  ! (in)
  & )
  
  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: xyz_RhoNl
  type(Field_IcGrid3D), intent(in) :: xyz_RhoBasic
  type(Field_IcGrid3D), intent(in) :: xyz_RNl

  ! 作業変数
  ! Local variables
  !
  integer :: k
  integer :: KzMax
  type(IcGrid3D_FVM), pointer :: icgrid3D
  type(Field_IcGrid2D), pointer :: xyz_RhoNl_k
  type(Field_IcGrid2D), pointer :: xyz_RhoBasic_k
  type(Field_IcGrid2D), pointer :: xyz_RNl_k

  ! 実行文 ; Executable statements
  !

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

  do k=1, KzMax
    xyz_RhoNl_k => get_field_layer(xyz_RhoNl, k)
    xyz_RhoBasic_k => get_field_layer(xyz_RhoBasic, k)
    xyz_RNl_k => get_field_layer(xyz_RNl, k)

    xyz_RhoNl_k%val(:,:,:,:) = xyz_RhoBasic_k%val(:,:,:,:) + xyz_RNl_k%val(:,:,:,:)

  end do

end subroutine calc_RhoNl

subroutine calc_hNl( &
  & xyz_hNl,            &  ! (inout)
  & xyz_ENl, xyz_RhoNl, &  ! (in)
  & Cp, Cv              &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: xyz_hNl
  type(Field_IcGrid3D), intent(in) :: xyz_ENl
  type(Field_IcGrid3D), intent(in) :: xyz_RhoNl
  real(DP), intent(in) :: Cp
  real(DP), intent(in) :: Cv
  
  ! 作業変数
  ! Local variables
  !
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(IcGrid3D_FVM), pointer :: icgrid3Dz
 
  type(Field_icGrid2D), pointer :: xyz_hNl_k
  type(Field_icGrid2D), pointer :: xyz_ENl_k
  type(Field_icGrid2D), pointer :: xyz_RhoNl_k
  
  integer :: k, KzMax
  
  ! 実行文 ; Executable statements
  !

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

  do k=1, KzMax
    xyz_hNl_k => get_field_layer(xyz_hNl, k)
    xyz_ENl_k => get_field_layer(xyz_ENl, k)
    xyz_RhoNl_k => get_field_layer(xyz_RhoNl, k)

    xyz_hNl_k%val(1, :,:,:) = ( Cp / Cv ) * xyz_ENl_k%val(1, :,:,:) &
              &                            / xyz_RhoNl_k%val(1, :,:,:)
  end do

end subroutine calc_hNl

subroutine calc_AtildeNl( &
  & AtildeNl,             &
  & xyz_RhoNl,            &
  & xyz_VhNl, xyr_WNl,    &
  & tfcoordz, derivate, vzderivate &
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: AtildeNl
  type(Field_IcGrid3D), intent(in) :: xyz_RhoNl
  type(Field_IcGrid3D), intent(in) :: xyz_VhNl
  type(Field_IcGrid3D), intent(in) :: xyr_WNl
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(in) :: vzderivate

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

  type(Field_IcGrid3D) :: AtildeNl_ei(3)
  type(Field_IcGrid3D) :: xyz_WNl

  type(Field_IcGrid2D), pointer :: AtildeNl_k
  type(Field_IcGrid2D), pointer :: AtildeNl_e1_k, AtildeNl_e2_k, AtildeNl_e3_k

  real(DP) :: e1(3) = (/ 1.0d0, 0.0d0, 0.0d0 /)
  real(DP) :: e2(3) = (/ 0.0d0, 1.0d0, 0.0d0 /)
  real(DP) :: e3(3) = (/ 0.0d0, 0.0d0, 1.0d0 /)

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

  ! 実行文 ; Executable statements
  !
 
  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 
  !
  
  icgrid3Dz => get_icgrid3D(xyz_VhNl)
  icgrid3Dr => get_icgrid3D(xyr_WNl)
  lpRCInfo => get_LocalProcRectInfo(get_icgrid2D(icgrid3Dz))
  
  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(AtildeNl_ei(1), icgrid3Dz, "AtildeNl_e1", 1)
  call Field_IcGrid3D_Init(AtildeNl_ei(2), icgrid3Dz, "AtildeNl_e2", 1)
  call Field_IcGrid3D_Init(AtildeNl_ei(3), icgrid3Dz, "AtildeNl_e3", 1)
  call Field_IcGrid3D_Init(xyz_WNl, icgrid3Dz, "xyz_WNl", 1)

  call interpolate_fieldr_to_z(xyr_WNl, xyz_WNl)

  do eiID=1, 3
    call calc_tmpAtildeNl_i( &
      & AtildeNl_ei(eiID), eiID, &
      & xyz_VhNl, xyr_WNl, xyz_WNl, xyz_RhoNl, &
      & icgrid3Dz, icgrid3Dr, tfcoordz, &
      & derivate, vzderivate )
  end do

  do kz=1, KzMax
    AtildeNl_k => get_field_layer(AtildeNl, kz)
    AtildeNl_e1_k => get_field_layer(AtildeNl_ei(1), kz)
    AtildeNl_e2_k => get_field_layer(AtildeNl_ei(2), kz)
    AtildeNl_e3_k => get_field_layer(AtildeNl_ei(3), kz)

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

          AtildeNl_k%val(1:3, i,j,rcID) = &
            &   AtildeNl_e1_k%val(1, i,j,rcID) * e1(1:3) &
            & + AtildeNl_e2_k%val(1, i,j,rcID) * e2(1:3) &
            & + AtildeNl_e3_k%val(1, i,j,rcID) * e3(1:3)

        end do
      end do
    end do
  end do

  do eiID=1, 3
    call Field_IcGrid3D_Final(AtildeNl_ei(eiID))
  end do

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

end subroutine calc_AtildeNl

subroutine calc_CtildeNl( &
  & CtildeNl,          &   ! (inout)
  & xyz_VhNl, xyr_WNl, &   ! (in)
  & Omega              &   ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: CtildeNl
  type(Field_IcGrid3D), intent(in) :: xyz_VhNl
  type(Field_IcGrid3D), intent(in) :: xyr_WNl
  real(DP), intent(in) :: Omega

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

  type(Field_icGrid2D), pointer :: CtildeNl_k
  type(Field_icGrid2D), pointer :: xyz_VhNl_k
  type(Field_icGrid2D), pointer :: xyz_WNl_k

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

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

  ! 実行文 ; Executable statements
  !
  
  ! 格子に関する情報を取得する. 
  ! Get the information about grid. 
  !
  
  icgrid3Dz => get_icgrid3D(xyz_VhNl)
  icgrid2D => get_icgrid2D(icgrid3Dz)
  lpRCInfo => get_LocalProcRectInfo(get_icgrid2D(icgrid3Dz))
  
  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_WNl", 1)

  call interpolate_fieldr_to_z(xyr_WNl, xyz_WNl)

  do k=1, KzMax
    xyz_VhNl_k => get_field_layer(xyz_VhNl, k)
    xyz_WNl_k => get_field_layer(xyz_WNl, k)
    CtildeNl_k => get_field_layer(CtildeNl, 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))
          univec_k(:) = local_k(s_pos(1), s_pos(2))
          Vei(1:3) = xyz_VhNl_k%val(1:3, i,j,rcID) + xyz_WNl_k%val(1, i,j,rcID) * univec_k(1:3)

          CtildeNl_k%val(1:3, i,j,rcID) = 2.0d0 * Omega * (/ - Vei(2), Vei(1), 0.0d0 /)

        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)

end subroutine calc_CtildeNl

subroutine calc_GR_GVh_GW( &
  & G_R, G_Vh, G_W, &
  & xyz_AtildeNl, xyz_CtildeNl, &
  & hydrostatic_option &
  & )

  type(Field_IcGrid3D), intent(inout) :: G_R
  type(Field_IcGrid3D), intent(inout) :: G_Vh  
  type(Field_IcGrid3D), intent(inout) :: G_W
  type(Field_IcGrid3D), intent(in) :: xyz_AtildeNl
  type(Field_IcGrid3D), intent(in) :: xyz_CtildeNl
  logical, intent(in) :: hydrostatic_option
  
  ! 作業変数
  ! Local variables
  !
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(IcGrid3D_FVM), pointer :: icgrid3Dz

  type(Field_IcGrid3D) :: xyz_GW

  type(Field_IcGrid2D), pointer :: xyz_GVh_k, xyz_GW_k
  type(Field_IcGrid2D), pointer :: AtildeNl_k
  type(Field_IcGrid2D), pointer :: CtildeNl_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, KzMax, rcIDMin, rcIDMax

  real(DP) :: AtildeCtildeNlh_rcijk(3), AtildeCtildeNlz_rcijk
  
  ! 実行文 ; Executable statements
  !

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

  icgrid3Dz => get_icgrid3D(xyz_AtildeNl)
  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_GW, icgrid3Dz,  'xyz_GW', 1)

  do k=1, KzMax
    xyz_GVh_k => get_field_layer(G_Vh, k)
    xyz_GW_k => get_field_layer(xyz_GW, k)
    AtildeNl_k => get_field_layer(xyz_AtildeNl, k)
    CtildeNl_k => get_field_layer(xyz_CtildeNl, k)

    do rcID=rcIDMin, rcIDMax
      !$omp parallel do private(i, s_pos, univec_k)
      do j=jEMin, jEMax
        do i=iEMin, iEMax

          s_pos(1:3) = orth_to_geo_pos(icgrid2D%rcs_AGrid(1:3, i,j,rcID))
          univec_k(1:3) = local_k(s_pos(1), s_pos(2))

          AtildeCtildeNlz_rcijk    = dot_product( AtildeNl_k%val(1:3, i,j,rcID) + CtildeNl_k%val(1:3, i,j,rcID), univec_k(1:3))
          AtildeCtildeNlh_rcijk(:) =   ( AtildeNl_k%val(1:3, i,j,rcID) + CtildeNl_k%val(1:3, i,j,rcID) ) &
              &                        - AtildeCtildeNlz_rcijk * univec_k(1:3)

          xyz_GVh_k%val(1:3, i,j,rcID) = - AtildeCtildeNlh_rcijk(:)

          if ( .not. hydrostatic_option ) then   
            xyz_GW_k%val(1, i,j,rcID) = - AtildeCtildeNlz_rcijk
          else
            xyz_GW_k%val(1, i,j,rcID) = 0d0  
          end if

        end do
      end do
    end do

  end do

  call interpolate_fieldz_to_r(xyz_GW, G_W)

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



subroutine apply_numerical_filter(& 
  & xyz_RAl, xyz_VhAl, xyr_WAl, xyz_EAl, xyz_RhoAl, &
  & xyz_HVDiff_R, xyz_HVDiff_Vh, xyr_HVDiff_W, xyz_HVDiff_Ea, &
  & tfcoordz, tfcoordr, &
  & derivate, vderivatez, vderivater, &
  & gamma_h, gamma_v, DelTimeLong & 
  & )

  type(Field_IcGrid3D), intent(in) :: xyz_RAl
  type(Field_IcGrid3D), intent(in) :: xyz_VhAl
  type(Field_IcGrid3D), intent(in) :: xyr_WAl
  type(Field_IcGrid3D), intent(in) :: xyz_EAl
  type(Field_IcGrid3D), intent(in) :: xyz_RhoAl    
  type(Field_IcGrid3D), intent(inout) :: xyz_HVDiff_R
  type(Field_IcGrid3D), intent(inout) :: xyz_HVDiff_Vh
  type(Field_IcGrid3D), intent(inout) :: xyr_HVDiff_W
  type(Field_IcGrid3D), intent(inout) :: xyz_HVDiff_Ea
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordr
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(inout) :: vderivatez
  type(VDerivate_Field_IcGrid3D), intent(inout) :: vderivater
  real(DP), intent(in) :: gamma_h
  real(DP), intent(in) :: gamma_v
  real(DP), intent(in) :: DelTimeLong
  
  type(Field_IcGrid3D) :: xyz_HVDiffTerm_sca, xyz_HVDiffTerm_vec, xyr_HVDiffTerm_sca
  type(Field_IcGrid3D) :: xyr_Rho
  type(IcGrid3D_FVM), pointer :: icgrid3Dz, icgrid3Dr
  integer :: k, KrMax, KzMax
  
  icgrid3Dz => get_icgrid3D(xyz_VhAl)
  icgrid3Dr => get_icgrid3D(xyr_WAl)
  
  KzMax = get_vertical_level_num(icgrid3Dz)
  KrMax = get_vertical_level_num(icgrid3Dr)
  
  call Field_IcGrid3D_Init(xyz_HVDiffTerm_sca, icgrid3Dz, 'xyz_HVDiffTerm', 1)
  call Field_IcGrid3D_Init(xyz_HVDiffTerm_vec, icgrid3Dz, 'xyz_HVDiffTerm', 3)
  call Field_IcGrid3D_Init(xyr_HVDiffTerm_sca, icgrid3Dr, 'xyz_HVDiffTerm', 1)

  !
  !
  call numFilter_HVDiffusion( &
  & xyz_HVDiffTerm_sca, xyz_RhoAl,  xyz_RhoAl, &
  & tfcoordz, derivate, vderivatez,  &
  & gamma_h, gamma_v, DelTimeLong, 1 &
  & )
        
  do k=1, KzMax
    xyz_HVDiff_R%field2Ds(k)%val(:,:,:,:) = &
      &  ( xyz_HVDiff_R%field2Ds(k)%val(:,:,:,:) + xyz_HVDiffTerm_sca%field2Ds(k)%val(:,:,:,:) ) 
  end do
  
    
  !
  !
  call numFilter_HVDiffusion( &
  & xyz_HVDiffTerm_vec, xyz_VhAl,  xyz_RhoAl, &
  & tfcoordz, derivate, vderivatez,  &
  & gamma_h, gamma_v, DelTimeLong, 2 &
  & )
  
  do k=1, KzMax
    xyz_HVDiff_Vh%field2Ds(k)%val(:,:,:,:) = &
      &  ( xyz_HVDiff_Vh%field2Ds(k)%val(:,:,:,:) + xyz_HVDiffTerm_vec%field2Ds(k)%val(:,:,:,:) ) 
  end do
  
  !
  !
  call numFilter_HVDiffusion( &
  & xyz_HVDiffTerm_sca, xyz_EAl,  xyz_RhoAl, &
  & tfcoordz, derivate, vderivatez,  &
  & gamma_h, gamma_v, DelTimeLong, 2 &
  & )

  do k=1, KzMax
    xyz_HVDiff_Ea%field2Ds(k)%val(:,:,:,:) = &
      &  ( xyz_HVDiff_Ea%field2Ds(k)%val(:,:,:,:) + xyz_HVDiffTerm_sca%field2Ds(k)%val(:,:,:,:) ) 
  end do

  
  !
  !
  call Field_IcGrid3D_Init(xyr_Rho, icgrid3Dr, 'xyr_Rho', 1)
  call interpolate_fieldz_to_r(xyz_RhoAl, xyr_Rho)

  call numFilter_HVDiffusion( &
    & xyr_HVDiffTerm_sca, xyr_WAl,  xyr_Rho, &
   & tfcoordr, derivate, vderivater,  &
   & gamma_h, gamma_v, DelTimeLong, 2 &
   & )

  do k=1, KrMax-1
    xyr_HVDiff_W%field2Ds(k)%val(:,:,:,:) = &
      & ( xyr_HVDiff_W%field2Ds(k)%val(:,:,:,:) + xyr_HVDiffTerm_sca%field2Ds(k)%val(:,:,:,:) ) 
  end do  
  
  call Field_IcGrid3D_Final(xyr_Rho)
  
  !
  call Field_IcGrid3D_Final(xyz_HVDiffTerm_sca)
  call Field_IcGrid3D_Final(xyz_HVDiffTerm_vec)
  call Field_IcGrid3D_Final(xyr_HVDiffTerm_sca)

  
end subroutine apply_numerical_filter

  
!> \~japanese
!! \brief \f$ \tilde{g}^t \f$ の場を計算する. 
!! \~english
!! \brief Calculate the field \f$ \tilde{g}^t \f$. 
!!
!!
!! \~japanese
!! @param[inout] gtildeNl
!! @param[in]    xyz_RhoNl
!! @param[in]    bouyancyNl
!! @param[in]    Grav        重力加速度. [m/s^2]
!!
!! \~english
!! @param[inout] gtildeNl
!! @param[in]    xyz_RhoNl
!! @param[in]    bouyancyNl
!! @param[in]    Grav        Gravitational acceleration. [m/s^2]
!!
subroutine calc_gtildeNl( &
  & gtildeNl,              &  ! (inout)
  & xyz_RhoNl, buoyancyNl, &  ! (in)
  & Grav                   &  ! (in)
  & )
  
  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: gtildeNl
  type(Field_IcGrid3D), intent(in) :: xyz_RhoNl
  type(Field_IcGrid3D), intent(in) :: buoyancyNl
  real(DP), intent(in) :: Grav

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

  type(Field_IcGrid2D), pointer :: gtildeNl_k
  type(Field_IcGrid2D), pointer :: buoyancyNl_k
  type(Field_IcGrid2D), pointer :: xyz_RhoNl_k

  real(DP) :: gamma_
  real(DP) :: radius

  integer :: KzMax
  integer :: k

  ! 実行文 ; Executable statements
  !

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

  icgrid3Dz => get_icgrid3D(buoyancyNl)
  icgrid2D => get_icgrid2D(icgrid3Dz)
  KzMax = get_vertical_level_num(icgrid3Dz)
  radius = get_IcRadius(icgrid2D)

  do k=1, KzMax
    gtildeNl_k => get_field_layer(gtildeNl, k)
    buoyancyNl_k => get_field_layer(buoyancyNl, k)
    xyz_RhoNl_k => get_field_layer(xyz_RhoNl, k)

    gtildeNl_k%val(:,:,:,:) =  &
            &    Grav &
            &  + buoyancyNl_k%val(:,:,:,:) / xyz_RhoNl_k%val(:,:,:,:) 

  end do

end subroutine calc_gtildeNl


subroutine calc_buoyancyNl( &
  & buoyancyNl,                &  ! (inout)
  & xyz_PNl, xyz_RNl,          &  ! (in)
  & tfcoordz, vzderivate, Grav &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: buoyancyNl
  type(Field_IcGrid3D), intent(in) :: xyz_PNl
  type(Field_IcGrid3D), intent(in) :: xyz_RNl
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  type(VDerivate_Field_IcGrid3D), intent(in) :: vzderivate
  real(DP), intent(in) :: Grav

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

  type(Field_IcGrid2D), pointer :: buoyancyNl_k
  type(Field_IcGrid2D), pointer :: xyz_PNl_k
  type(Field_IcGrid2D), pointer :: xyz_RNl_k

  type(Field_IcGrid3D) :: tmp_Vdel_field3D
  type(Field_IcGrid2D), pointer :: tmp_Vdel_field3D_k
  type(Field_IcGrid3D) :: Vdel_field3D
  type(Field_IcGrid2D), pointer :: Vdel_field3D_k

  real(DP) :: gamma_
  real(DP) :: radius

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

  ! 実行文 ; Executable statements
  !

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

  icgrid3Dz => get_icgrid3D(xyz_PNl)
  icgrid2D => get_icgrid2D(icgrid3Dz)
  lpRCInfo => get_LocalProcRectInfo(get_icgrid2D(icgrid3Dz))
  
  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(Vdel_field3D, icgrid3Dz,  'Vdel_field3D', 1)
  call Field_IcGrid3D_Init(tmp_Vdel_field3D, icgrid3Dz,  'tmp_Vdel_field3D', 1)

  do k=1, KzMax
    xyz_PNl_k => get_field_layer(xyz_PNl, k)
    tmp_Vdel_field3D_k => get_field_layer(tmp_Vdel_field3D, k)

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

          gamma_ = calc_gamma(TF_to_Z_IcGrid3D(tfCoordz, icgrid3Dz%verticalGrid(k), rcID, i, j), radius)

          tmp_Vdel_field3D_k%val(1, i,j,rcID) = xyz_PNl_k%val(1, i,j,rcID) / ( tfcoordz%G_1div2%val(1, i,j,rcID) * gamma_**2 )

        end do
      end do
    end do
  end do

  call d_dVCoord(vzderivate, tmp_Vdel_field3D, Vdel_field3D)

  do k=1, KzMax
    buoyancyNl_k => get_field_layer(buoyancyNl, k)
    xyz_RNl_k => get_field_layer(xyz_RNl, k)
    Vdel_field3D_k => get_field_layer(Vdel_field3D, k)

    do rcID=rcIDMin, rcIDMax
      !$omp parallel do private(i, gamma_)
      do j=jEMin, jEMax
        do i=iEMin, iEMax
 
          gamma_ = calc_gamma(TF_to_Z_IcGrid3D(tfCoordz, icgrid3Dz%verticalGrid(k), rcID, i, j), radius)

          buoyancyNl_k%val(1, i,j,rcID) =  - ( gamma_**2 * Vdel_field3D_k%val(1, i,j,rcID) &
                                  &        + xyz_RNl_k%val(1, i,j,rcID) * Grav )
        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(Vdel_field3D)
  call Field_IcGrid3D_Final(tmp_Vdel_field3D)

end subroutine calc_buoyancyNl

subroutine calc_HGradPNl( &
  & xyz_HGradPNl, &
  & xyz_PNl,      &
  & tfcoordz, derivate, vzderivate &
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: xyz_HGradPNl
  type(Field_IcGrid3D), intent(in) :: xyz_PNl
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(in) :: vzderivate

  ! 実行文 ; Executable statements
  !
  
  call TF_HGradient( &
    & xyz_PNl, xyz_HGradPNl, &
    & tfcoordz, derivate, vzderivate )

end subroutine calc_HGradPNl

subroutine update_variable( &
  & variableAl,             &
  & variableNl, variableNs  &
  & )

  ! 宣言文 ; Declaration statements
  !  
  type(Field_IcGrid3D), intent(inout) :: variableAl
  type(Field_IcGrid3D), intent(in) :: variableNl
  type(Field_IcGrid3D), intent(in) :: variableNs

  ! 作業変数
  ! Local variables
  !
  integer :: k
  integer :: vtLvNum
  type(IcGrid3D_FVM), pointer :: icgrid3D
  type(Field_IcGrid2D), pointer :: variableAl_k, variableNl_k, variableNs_k

  ! 実行文 ; Executable statements
  !

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

  !$omp parallel do private(variableAl_k, variableNl_k, variableNs_k)
  do k=1, vtLvNum
    variableAl_k => get_field_layer(variableAl, k)
    variableNl_k => get_field_layer(variableNl, k)
    variableNs_k => get_field_layer(variableNs, k)

    variableAl_k%val(:,:,:,:) = variableNl_k%val(:,:,:,:) + variableNs_k%val(:,:,:,:)

  end do

end subroutine update_variable

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! モジュール外非公開手続き
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine calc_tmpAtildeNl_i( &
  & AtildeNl_i, &
  & ei_ID,                                  &
  & xyz_VhNl, xyr_WNl,  xyz_WNl, xyz_RhoNl, &
  & icgrid3Dz, icgrid3Dr, tfCoordz,         &
  & derivate, vzderivate &
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: AtildeNl_i
  integer, intent(in) :: ei_ID
  type(Field_IcGrid3D), intent(in) :: xyz_VhNl
  type(Field_IcGrid3D), intent(in) :: xyr_WNl
  type(Field_IcGrid3D), intent(in) :: xyz_WNl
  type(Field_IcGrid3D), intent(in) :: xyz_RhoNl
  type(IcGrid3D_FVM), intent(in) :: icgrid3Dz
  type(IcGrid3D_FVM), intent(in) :: icgrid3Dr
  type(TFCoord_IcGrid3D), intent(in) :: tfCoordz
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(in) :: vzderivate

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

  type(Field_IcGrid3D) :: hdiv_field
  type(Field_IcGrid2D), pointer :: hdiv_field_k
  type(Field_IcGrid3D) :: tmp_hdiv_field
  type(Field_IcGrid2D), pointer :: tmp_hdiv_field_k

  type(Field_IcGrid3D) :: xyz_Vei, xyz_VDelVeiWFlux
  type(Field_IcGrid2D), pointer :: xyz_Vei_k, xyz_VDelVeiWFlux_k

  type(Field_IcGrid2D), pointer :: AtildeNl_i_k

  type(Field_IcGrid2D), pointer :: xyz_VhNl_k
  type(Field_IcGrid2D), pointer :: xyz_WNl_k
  type(Field_IcGrid2D), pointer :: xyz_RhoNl_k

  real(DP) :: v_ei_rij, VhNl_rijk(3)
  real(DP) :: xi
  real(DP) :: gamma_
  real(DP) :: radius
  real(DP) :: s_pos(3)
  real(DP) :: uvec_k(3)

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

  ! 実行文 ; Executable statements
  !

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

  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(hdiv_field, icgrid3Dz, "hdiv_field", 1)
  call Field_IcGrid3D_Init(tmp_hdiv_field, icgrid3Dz, "tmp_hdiv_field", 3)
  call Field_IcGrid3D_Init(xyz_Vei, icgrid3Dz, "xyz_Vei", 1)
  call Field_IcGrid3D_Init(xyz_VDelVeiWFlux, icgrid3Dz, "xyz_VDelVeiWFlux", 1)

  do k=1, KzMax
    xyz_VhNl_k => get_field_layer(xyz_VhNl, k)
    xyz_WNl_k => get_field_layer(xyz_WNl, k)
    xyz_RhoNl_k => get_field_layer(xyz_RhoNl, k)

    tmp_hdiv_field_k => get_field_layer(tmp_hdiv_field, k)
    xyz_Vei_k => get_field_layer(xyz_Vei, k)
    
    xi = icgrid3Dz%verticalGrid(k)
    
    do rcID=rcIDMin, rcIDMax
      !$omp parallel do private(i, gamma_, s_pos, uvec_k, v_ei_rij, VhNl_rijk)
      do j=jEMin, jEMax
        do i=iEMin, iEMax

          gamma_ = calc_gamma(TF_to_Z_IcGrid3D(tfCoordz, xi, rcID, i, j), radius)

          s_pos(:) = orth_to_geo_pos(icgrid2D%rcs_AGrid(:, i,j,rcID))
          uvec_k(:) = local_k(s_pos(1), s_pos(2))

          VhNl_rijk(1:3) = xyz_VhNl_k%val(1:3, i,j,rcID)
          v_ei_rij =  ( VhNl_rijk(ei_ID) + xyz_WNl_k%val(1, i,j,rcID) * uvec_k(ei_ID)) &
                  &      / xyz_RhoNl_k%val(1, i,j,rcID)

          tmp_hdiv_field_k%val(1:3, i,j,rcID) = v_ei_rij * VhNl_rijk(1:3)

          xyz_Vei_k%val(1, i,j,rcID) = v_ei_rij

        end do
      end do
    end do

  end do

  call TF_HDivergence( &
    & tmp_hdiv_field, hdiv_field, &
    & tfcoordz, derivate, vzderivate )

  call TF_VDerivateVFlux( &
    & xyr_WNl, xyz_Vei, &
    & xyz_VDelVeiWFlux, tfcoordz)

 
  !$omp parallel do private (hdiv_field_k, xyz_VDelVeiWFlux_k, AtildeNl_i_k)
  do k=1, KzMax
    hdiv_field_k => get_field_layer(hdiv_field, k)
    xyz_VDelVeiWFlux_k => get_field_layer(xyz_VDelVeiWFlux, k)
    AtildeNl_i_k => get_field_layer(AtildeNl_i, k)

    AtildeNl_i_k%val(1,:,:,:) =   hdiv_field_k%val(1,:,:,:)      &
                      &         + xyz_VDelVeiWFlux_k%val(1,:,:,:) 
  end do

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

end subroutine calc_tmpAtildeNl_i

end module large_time_step_RK2_mod
