module numerical_filter_mod

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

#ifdef LIB_MPI
  use mpi
#endif

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

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

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

  use igmcore_math, only: &
    & PI
  
  ! 座標変換
  ! 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, get_glevel, &
    & 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, &
    & Field_IcGrid2D_Init, Field_IcGrid2D_Final, &
    & 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_field_layer, get_field3D_rank, &
    & get_field3d_name, 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, laplacian_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, dd_dVCoord

  !* DIGMO ***********************
  !*
  
  !
  !
  !
  use Vertical_Interpolation, only: &
    & interpolate_fieldr_to_z, interpolate_fieldz_to_r


use IcGrid_ncWriter_mod

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public statements
  !
  public :: numFilter_HVDiffusion 

logical :: init_writer_flag = .true.
integer :: fieldID1, fieldID2
integer :: tcount = 1
type(IcGrid_ncWriter),save :: writer

contains
subroutine numFilter_HVDiffusion( &
  & HVDiffTerm, progVar,  Rho, &
  & tfcoord, derivate, vderivate,  &
  & gamma_h, gamma_v, DelTimeLong, &
  & optr_type &
  & )
  
  type(Field_IcGrid3D), intent(inout) :: HVDiffTerm
  type(Field_IcGrid3D), intent(in) :: progVar
  type(Field_IcGrid3D), intent(in) :: Rho
  type(TFCoord_IcGrid3D), intent(in) :: tfcoord
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(in) :: vderivate
  real(DP), intent(in) :: gamma_h
  real(DP), intent(in) :: gamma_v
  real(DP), intent(in) :: DelTimeLong
  integer, intent(in) :: optr_type

  type(IcGrid3D_FVM), pointer :: icgrid3D
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(Field_IcGrid3D) :: phi, mu

  type(LocalProcRectInfo), pointer :: lpRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, KMax, rcIDMin, rcIDMax
  integer :: rcID, i, j, k, rankID
  
  real(DP) :: xi, gamma_, radius
  real(DP) :: K_H
  real(DP), allocatable :: K_V(:)
 integer :: fieldRank
  
  icgrid3D => get_icgrid3D(progVar)
  icgrid2D => get_icgrid2D(icgrid3D)
  lpRCInfo => get_LocalProcRectInfo(icgrid2D)
  fieldRank = get_field3D_rank(progVar)
  
  call get_InteriorRegion_IndexInfo( &
    & lpRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out) 

  KMax = get_vertical_level_num(icgrid3D)
  radius = get_IcRadius(icgrid2D)
  
  call Field_IcGrid3D_Init(phi, icgrid3D, 'phi', fieldRank)
  call Field_IcGrid3D_Init(mu, icgrid3D, 'mu', 1)  

  allocate(K_V(KMax))
  
  call calc_HVDiffusionCoef(icgrid3D, gamma_h, gamma_v, DelTimeLong, K_H, K_V)


  select case (optr_type)
  case (1)

    do k=1, KMax
      xi = icgrid3D%verticalGrid(k)    
      do rcID=rcIDMin, rcIDMax
        do j=jEMin, jEMax
          do i=iEMin, iEMax
            gamma_ = calc_gamma(TF_to_Z_IcGrid3D(tfCoord, xi, rcID, i, j), radius)          
            mu%field2Ds(k)%val(1,i,j,rcID) = tfcoord%G_1div2%val(1,i,j,rcID) * gamma_**2

          end do
        end do
      end do
    end do 

  case (2)
    call copy_field3D(mu, Rho)
  end select

  do k=1, KMax
    do rankID=1,fieldRank
      phi%field2Ds(k)%val(rankID,:,:,:) = &
        & progVar%field2Ds(k)%val(rankID,:,:,:) / mu%field2Ds(k)%val(1,:,:,:) 
    end do
  end do

  call calc_HVDiffTerm(HVDiffTerm, mu, phi, K_H, K_V, derivate, vderivate)

  !
  call Field_IcGrid3D_Final(phi)
  call Field_IcGrid3D_Final(mu)
    
end subroutine numFilter_HVDiffusion

subroutine calc_HVDiffTerm(HVDiffTerm, mu, phi, K_H, K_V, derivate, vderivate)
  type(Field_IcGrid3D), intent(inout) :: HVDiffTerm
  type(Field_IcGrid3D), intent(in) :: mu
  type(Field_IcGrid3D), intent(in) :: phi
  real(DP), intent(in) :: K_H
  real(DP), intent(in) :: K_V(:)
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(in) :: vderivate
  
  
  type(LocalProcRectInfo), pointer :: lpRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, KMax, rcIDMin, rcIDMax
  integer :: rcID, i, j, k, rankID
  integer :: fieldRank

  type(Field_IcGrid2D) :: hdel2
  type(Field_IcGrid3D) :: vdel_tmp1, vdel_tmp2, hdel4
  type(Field_IcGrid2D), pointer :: phi_k, mu_k, HVDiffTerm_k
  type(Field_IcGrid2D), pointer :: hdel4_k, vdel_tmp1_k, vdel_tmp2_k
 
  type(IcGrid3D_FVM), pointer :: icgrid3D
  type(IcGrid2D_FVM), pointer :: icgrid2D
  
  icgrid3D => get_icgrid3D(phi)
  icgrid2D => get_icgrid2D(icgrid3D)
  lpRCInfo => get_LocalProcRectInfo(icgrid2D)
  fieldRank = get_field3D_rank(phi)
  
  call get_InteriorRegion_IndexInfo( &
    & lpRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out) 

  KMax = get_vertical_level_num(icgrid3D)
  
  call Field_IcGrid2D_Init(hdel2, icgrid2D, 'hdel2', fieldRank)
  call Field_IcGrid3D_Init(hdel4, icgrid3D, 'hdel4', fieldRank)
  call Field_IcGrid3D_Init(vdel_tmp1, icgrid3D, 'vdel_tmp1', fieldRank)
  call Field_IcGrid3D_Init(vdel_tmp2, icgrid3D, 'vdel_tmp2', fieldRank)

  call dd_dVCoord(vderivate, phi, vdel_tmp1)
  call dd_dVCoord(vderivate, vdel_tmp1, vdel_tmp2)
  
  do k=1, KMax
    phi_k => get_field_layer(phi, k)
    mu_k => get_field_layer(mu, k)
    hdel4_k => get_field_layer(hdel4, k)
    vdel_tmp1_k => get_field_layer(vdel_tmp1, k)
    vdel_tmp2_k => get_field_layer(vdel_tmp2, k)

    call laplacian_op(derivate, phi_k, hdel2)
    
    do rankID=1,fieldRank
      hdel2%val(rankID,:,:,:) = K_H * mu_k%val(1, :,:,:) * hdel2%val(rankID,:,:,:)
      vdel_tmp1_k%val(rankID,:,:,:) = K_V(k) * mu_k%val(1, :,:,:) * vdel_tmp2_k%val(rankID,:,:,:)
    end do
    
    call laplacian_op(derivate, hdel2, hdel4_k)
  end do

  call dd_dVCoord(vderivate, vdel_tmp1, vdel_tmp2)
  
  do k=1, KMax
    vdel_tmp2_k => get_field_layer(vdel_tmp2, k)
    hdel4_k => get_field_layer(hdel4, k)
    HVDiffTerm_k => get_field_layer(HVDiffTerm, k)
  
    HVDiffTerm_k%val(:,:,:,:) = - hdel4_k%val(:,:,:,:) + vdel_tmp2_k%val(:,:,:,:)

  end do
  
  call Field_IcGrid3D_Final(vdel_tmp1)
  call Field_IcGrid3D_Final(vdel_tmp2)
  call Field_IcGrid3D_Final(hdel4)
  call Field_IcGrid2D_Final(hdel2)

end subroutine calc_HVDiffTerm

subroutine calc_HVDiffusionCoef(icgrid3D, gamma_h, gamma_v, DelTime, K_H, K_V)
  type(IcGrid3D_FVM), intent(in) :: icgrid3D
  real(DP), intent(in) :: gamma_h
  real(DP), intent(in) :: gamma_v
  real(DP), intent(in) :: DelTime
  real(DP), intent(out) :: K_H
  real(DP), intent(inout) :: K_V(:)
  
  integer :: k, vtLvNum
  real(DP) :: DelXi(size(K_V))
  reaL(DP) :: mean_grid_intrv, Ng
  type(IcGrid2D_FVM), pointer :: icgrid2D
  
  vtLvNum = get_vertical_level_num(icgrid3D)
  icgrid2D => get_icgrid2D(icgrid3D)
  
  DelXi(1) = icgrid3D%verticalGrid(2) - icgrid3D%verticalGrid(1) 
  do k=2, vtLvNum-1
    DelXi(k) = ( icgrid3D%verticalGrid(k+1) - icgrid3D%verticalGrid(k-1) ) / 2d0
  end do
  DelXi(vtLvNum) = icgrid3D%verticalGrid(vtLvNum) - icgrid3D%verticalGrid(vtLvNum-1)
  
  Ng = 10 * 4**(get_glevel(icgrid2D)) + 2
  mean_grid_intrv = sqrt(4.0d0 * PI * get_IcRadius(icgrid2D)**2 / Ng )
  K_H = gamma_h * mean_grid_intrv**4 / DelTime
  
  do k=1, vtLvNum
    K_V(k) = gamma_v * DelXi(k)**6 / DelTime
  end do

end subroutine calc_HVDiffusionCoef

end module numerical_filter_mod
