!
!
! \f[
!   G^{1/2} \nabala_{h0} {\gamma \phi} =   \tilde{\nabla}_{h0} ( G^{1/2} \gamma \phi )
!                                        + \DP{}{\xi} ( \Dvect{G}^{\Dvect{z}} G^{1/2} \gamma \phi )
! \f]
!
! 右辺について, \f$\Phi=\phi G^{1/2} \gamma^2\f$なる変数変換を行う.
!
! \f[
!   G^{1/2} \nabala_{h0} {\gamma \phi} =   \tilde{\nabla}_{h0} ( \dfrac{\Phi}{\gamma} )
!                                        + \DP{}{\xi} ( \Dvect{G}^{\Dvect{z}} \dfrac{\Phi}{\gamma} )
! \f]
!
! 手続きの引数として与える構造型 Field_IcGrid3D の変数は, 物理場の変数\f$ \phi \f$(速度, 内部エネルギー)そのものではなく,
! それにメトリックと密度を含めた, 上で導入された\f$ \Phi \f$である.
!
!
module TF_HVDerivateTool

  !* gtool ****************
  !*
  
  ! 種類型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: &
    & DP

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

  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

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

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

  ! 正二十面体格子上の物理場データの管理
  ! Manager of the physical field data on icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, Field_IcGrid2D_Final, &
    & paste_field2D_margin

  ! 3 次元正二十面体格子上の物理場データの管理
  ! Manager of the physical field data 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

  use Derivate_Field_IcGrid2D_Manager, only: &
    & Derivate_Field_IcGrid2D, &
    & divergence_op, gradient_op

  use VDerivate_Field_IcGrid3D_Manager, only: &
    & VDerivate_Field_IcGrid3D, d_dVCoord

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

  use Utility_GeoCoord, only: &
    & local_k

  ! OpenMP
  !
  !
  use omp_lib

  implicit none
  private

  public :: TF_HGradient, TF_HDivergence, TF_VDerivate, TF_VDerivateVFlux

contains


subroutine TF_HGradient( &
  & Phi, ret_del_Phi, &
  & tfCoord, &
  & derivate, &
  & vderivate &
  & )

  type(Field_IcGrid3D), intent(in) :: Phi
  type(Field_IcGrid3D), intent(inout) :: ret_del_Phi
  type(TFCoord_IcGrid3D), intent(in) :: tfCoord
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(in) :: vderivate


  type(Field_IcGrid2D) :: tmp_hgradient_Phi_k
  type(Field_IcGrid3D) :: hgradient_Phi
  type(Field_IcGrid3D) :: tmp_vderivate_Phi
  type(Field_IcGrid3D) :: vderivate_Phi

  type(Field_IcGrid2D), pointer :: Phi_k
  type(Field_IcGrid2D), pointer :: ret_del_Phi_k
  type(Field_IcGrid2D), pointer :: hgradient_Phi_k
  type(Field_IcGrid2D), pointer :: tmp_vderivate_Phi_k
  type(Field_IcGrid2D), pointer :: vderivate_Phi_k
  type(Field_IcGrid2D), pointer :: G_Z_k

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

  real(DP) :: Phi_rijk, gamma_
  real(DP) :: radius, xi

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

  icgrid3D => get_icgrid3D(Phi)
  icgrid2D => get_icgrid2D(icgrid3D)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)
  
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)
  KMax = get_vertical_level_num(icgrid3D)

  radius = get_IcRadius(icgrid2D)

  call Field_IcGrid2D_Init(tmp_hgradient_Phi_k, icgrid2D, 'tmp_hgradient_Phi_k', 1)
  call Field_IcGrid3D_Init(hgradient_Phi, icgrid3D, 'hdvi_Phi', 3)
  call Field_IcGrid3D_Init(tmp_vderivate_Phi, icgrid3D, 'tmp_vderivate_Phi', 3)
  call Field_IcGrid3D_Init(vderivate_Phi, icgrid3D, 'vderivate_Phi', 3)

  do k=1, KMax

    Phi_k => get_field_layer(Phi, k)
    hgradient_Phi_k => get_field_layer(hgradient_Phi, k)
    vderivate_Phi_k => get_field_layer(vderivate_Phi, k)
    tmp_vderivate_Phi_k => get_field_layer(tmp_vderivate_Phi, k)
    G_Z_k => get_field_layer(tfCoord%G_Z, k)
    
    xi = icgrid3D%verticalGrid(k)

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

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

          Phi_rijk = Phi_k%val(1, i,j,rcID)

          tmp_hgradient_Phi_k%val(1, i,j,rcID) = Phi_rijk / gamma_
          tmp_vderivate_Phi_k%val(1:3, i,j,rcID) = Phi_rijk * G_Z_k%val(1:3, i,j,rcID) / gamma_

        end do
      end do
    end do

    call paste_field2D_margin(tmp_hgradient_Phi_k)
    call gradient_op(derivate, tmp_hgradient_Phi_k, hgradient_Phi_k)

  end do

  call d_dVCoord(vderivate, tmp_vderivate_Phi, vderivate_Phi)

  !$omp parallel do private(ret_del_Phi_k, hgradient_Phi_k, vderivate_Phi_k)
  do k=1, kMax
    ret_del_Phi_k => get_field_layer(ret_del_Phi, k)
    hgradient_Phi_k => get_field_layer(hgradient_Phi, k)
    vderivate_Phi_k => get_field_layer(vderivate_Phi, k)

    ret_del_Phi_k%val(:,:,:,:) = hgradient_Phi_k%val(:,:,:,:) + vderivate_Phi_k%val(:,:,:,:)

  end do

  call Field_IcGrid2D_Final(tmp_hgradient_Phi_k)
  call Field_IcGrid3D_Final(hgradient_Phi)
  call Field_IcGrid3D_Final(tmp_vderivate_Phi)
  call Field_IcGrid3D_Final(vderivate_Phi)

end subroutine TF_HGradient

subroutine TF_HDivergence( &
  & Phi, ret_del_Phi, &
  & tfCoord, &
  & derivate, &
  & vderivate &
  & )

  type(Field_IcGrid3D), intent(in) :: Phi
  type(Field_IcGrid3D), intent(inout) :: ret_del_Phi
  type(TFCoord_IcGrid3D), intent(in) :: tfCoord
  type(Derivate_Field_IcGrid2D), intent(inout) :: derivate
  type(VDerivate_Field_IcGrid3D), intent(in) :: vderivate

  type(Field_IcGrid3D) :: hdiv_Phi
  type(Field_IcGrid3D) :: tmp_vderivate_Phi
  type(Field_IcGrid3D) :: vderivate_Phi

  type(Field_IcGrid2D), pointer :: Phi_k
  type(Field_IcGrid2D), pointer :: ret_del_Phi_k
  type(Field_IcGrid2D) :: tmp_hdiv_Phi_k
  type(Field_IcGrid2D), pointer :: hdiv_Phi_k
  type(Field_IcGrid2D), pointer :: tmp_vderivate_Phi_k
  type(Field_IcGrid2D), pointer :: vderivate_Phi_k
  type(Field_IcGrid2D), pointer :: G_Z_k

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

  real(DP) :: Phi_rijk(3), gamma_
  real(DP) :: radius
  real(DP) :: s_pos(3), xi

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

  icgrid3D => get_icgrid3D(Phi)
  icgrid2D => get_icgrid2D(icgrid3D)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)
  
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)
  KMax = get_vertical_level_num(icgrid3D)

  radius = get_IcRadius(icgrid2D)

  call Field_IcGrid2D_Init(tmp_hdiv_Phi_k, icgrid2D, 'tmp_hdvi_Phi_k', 3)
  call Field_IcGrid3D_Init(hdiv_Phi, icgrid3D, 'hdvi_Phi', 1)
  call Field_IcGrid3D_Init(tmp_vderivate_Phi, icgrid3D, 'tmp_vderivate_Phi', 1)
  call Field_IcGrid3D_Init(vderivate_Phi, icgrid3D, 'vderivate_Phi', 1)

  do k=1, KMax

    Phi_k => get_field_layer(Phi, k)
    hdiv_Phi_k => get_field_layer(hdiv_Phi, k)
    vderivate_Phi_k => get_field_layer(vderivate_Phi, k)
    tmp_vderivate_Phi_k => get_field_layer(tmp_vderivate_Phi, k)
    G_Z_k => get_field_layer(tfCoord%G_Z, k)

    xi = icgrid3D%verticalGrid(k)
    
   do rcID=rcIDMin, rcIDMax
     !$omp parallel do private(i, gamma_, Phi_rijk)
     do j=jEMin, jEMax
       do i=iEMin, iEMax

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

          Phi_rijk(1:3) = Phi_k%val(1:3, i,j,rcID)

          tmp_hdiv_Phi_k%val(1:3, i,j,rcID) = Phi_rijk(1:3) / gamma_
          tmp_vderivate_Phi_k%val(1, i,j,rcID) = dot_product(Phi_rijk(1:3), G_Z_k%val(1:3, i,j,rcID)) / gamma_

        end do
      end do
    end do

    call paste_field2D_margin(tmp_hdiv_Phi_k)
    call divergence_op(derivate, tmp_hdiv_Phi_k, hdiv_Phi_k)

  end do

  call d_dVCoord(vderivate, tmp_vderivate_Phi, vderivate_Phi)

  !$omp parallel do private(ret_del_Phi_k, hdiv_Phi_k, vderivate_Phi_k)
  do k=1, kMax
    ret_del_Phi_k => get_field_layer(ret_del_Phi, k)
    hdiv_Phi_k => get_field_layer(hdiv_Phi, k)
    vderivate_Phi_k => get_field_layer(vderivate_Phi, k)

    ret_del_Phi_k%val(1,:,:,:) = hdiv_Phi_k%val(1,:,:,:) + vderivate_Phi_k%val(1,:,:,:)

  end do

  call Field_IcGrid2D_Final(tmp_hdiv_Phi_k)
  call Field_IcGrid3D_Final(hdiv_Phi)
  call Field_IcGrid3D_Final(tmp_vderivate_Phi)
  call Field_IcGrid3D_Final(vderivate_Phi)

end subroutine TF_HDivergence

subroutine TF_VDerivate( &
  & Phi, ret_del_Phi, &
  & tfCoord, vderivate &
  & )

  type(Field_IcGrid3D), intent(in) :: Phi
  type(Field_IcGrid3D), intent(out) :: ret_del_Phi
  type(TFCoord_IcGrid3D), intent(in) :: tfCoord
  type(VDerivate_Field_IcGrid3D), intent(in) :: vderivate

  type(Field_IcGrid3D) :: tmp_vderivate_Phi

  type(Field_IcGrid2D), pointer :: Phi_k
  type(Field_IcGrid2D), pointer :: ret_del_Phi_k
  type(Field_IcGrid2D), pointer :: tmp_vderivate_Phi_k

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

  integer :: rcID, i, j, k
  integer :: iEMin, iEMax, jEMin, jEMax, KMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  integer :: rank

  icgrid3D => get_icgrid3D(Phi)
  icgrid2D => get_icgrid2D(icgrid3D)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)
  
  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out)
  KMax = get_vertical_level_num(icgrid3D)
  rank = get_field3D_rank(Phi)

  call Field_IcGrid3D_Init(tmp_vderivate_Phi, icgrid3D, 'tmp_vderivate_Phi', rank)

  do k=1, KMax

    Phi_k => get_field_layer(Phi, k)
    tmp_vderivate_Phi_k => get_field_layer(tmp_vderivate_Phi, k)

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

          tmp_vderivate_Phi_k%val(1:rank, i,j,rcID) = Phi_k%val(1:rank, i,j,rcID) / tfCoord%G_1div2%val(1, i,j,rcID)

        end do
      end do
    end do

  end do

  call d_dVCoord(vderivate, tmp_vderivate_Phi, ret_del_Phi)


  call Field_IcGrid3D_Final(tmp_vderivate_Phi)

end subroutine TF_VDerivate

subroutine TF_VDerivateVFlux(xyr_W, xyz_Phi, ret_xyz_VDelVFlux, tfcoordz)
  
  type(Field_IcGrid3D), intent(in) :: xyr_W
  type(Field_IcGrid3D), intent(in) :: xyz_Phi
  type(Field_IcGrid3D), intent(inout) :: ret_xyz_VDelVFlux
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  
  type(Field_IcGrid3D) :: xyr_Omega
  type(Field_IcGrid2D), pointer :: xyr_Omega_k, xyr_W_kU, xyr_W_kL, xyr_W_k
  type(Field_IcGrid2D), pointer :: xyz_Phi_kU, xyz_Phi_k, xyz_Phi_kL
  type(Field_IcGrid2D), pointer :: ret_xyz_VDelVFlux_k

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

  integer :: rcID, i, j, k, rankID
  integer :: iEMin, iEMax, jEMin, jEMax, KzMax, KrMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  integer :: rank
  
  real(DP), allocatable :: Xir(:), Xiz(:)
  real(DP) :: DelXi_k
  real(DP) :: xyz_Phi_rijkU(3), xyz_Phi_rijkL(3)
  real(DP) :: Phi_rijKz(3), Phi_rijKrU(3), Phi_rijKrL(3)
  
  real(DP) :: DelXiUU, DelXiU
  real(DP) :: DelXiL, DelXiLL

  icgrid3Dz => get_icgrid3D(xyz_Phi)
  icgrid3Dr => get_icgrid3D(xyr_W)
  icgrid2D => get_icgrid2D(icgrid3Dz)
  lcRCInfo => get_LocalProcRectInfo(icgrid2D)
  

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

  KzMax = get_vertical_level_num(icgrid3Dz)
  KrMax = get_vertical_level_num(icgrid3Dr)
  rank = get_field3D_rank(xyz_Phi)
  
  allocate(Xir(KrMax))
  allocate(Xiz(KzMax))

  Xir(:) = icgrid3Dr%verticalGrid(:)
  Xiz(:) = icgrid3Dz%verticalGrid(:)  
  
  !
  !
  
  do k=2, KzMax-1
  
    xyr_W_kU => get_field_layer(xyr_W, k+1)
    xyr_W_kL => get_field_layer(xyr_W, k)
    xyz_Phi_kU => get_field_layer(xyz_Phi, k+1)
    xyz_Phi_k => get_field_layer(xyz_Phi, k)
    xyz_Phi_kL => get_field_layer(xyz_Phi, k-1)
    ret_xyz_VDelVFlux_k => get_field_layer(ret_xyz_VDelVFlux, k)
    
    DelXi_k = Xir(k+1) - Xir(k)
    
    do rcID=rcIDMin, rcIDMax
      !$omp parallel do private( Phi_rijKz, Phi_rijKrU, Phi_rijKrL) 
      do j=jEMin, jEMax
        do i=iEMin, iEMax
        
            Phi_rijKz(1:rank) = xyz_Phi_k%val(1:rank, i,j,rcID)
            Phi_rijKrU(1:rank) = 0.5d0 * ( xyz_Phi_kU%val(1:rank, i,j,rcID) + Phi_rijKz(1:rank) )
            Phi_rijKrL(1:rank) = 0.5d0 * ( xyz_Phi_kL%val(1:rank, i,j,rcID) + Phi_rijKz(1:rank) )
          
            ret_xyz_VDelVFlux_k%val(1:rank, i,j,rcID) = &
              & ( &
              &     xyr_W_kU%val(1, i,j,rcID) * Phi_rijKrU(1:rank) &
            &   - xyr_W_kL%val(1, i,j,rcID) * Phi_rijKrL(1:rank) &
            & ) / ( DelXi_k * tfcoordz%G_1div2%val(1, i,j,rcID) )

        end do
      end do
    end do
    
  end do

  !
  !
  
  xyr_W_kU => get_field_layer(xyr_W, 2)
  xyr_W_kL => get_field_layer(xyr_W, 1)
  xyz_Phi_kU => get_field_layer(xyz_Phi, 2)
  xyz_Phi_k => get_field_layer(xyz_Phi, 1)
  DelXiUU = Xiz(2) - Xiz(1)
  DelXiU = Xiz(1) - Xir(1)

  DelXi_k = Xir(2) - Xir(1)
  
  do rcID=rcIDMin, rcIDMax
    !$omp parallel do private(xyz_Phi_rijkL)
    do j=jEMin, jEMax
      do i=iEMin, iEMax
        
        xyz_Phi_rijkL(1:rank) = - ( xyz_Phi_kU%val(:, i,j,rcID) - xyz_Phi_k%val(:, i,j,rcID) ) * DelXiU / DelXiUU &
          &                     + xyz_Phi_k%val(:, i,j,rcID)  
          
        ret_xyz_VDelVFlux%field2Ds(1)%val(:, i,j,rcID) = &
            & ( &
            &     xyr_W_kU%val(1, i,j,rcID) * ( xyz_Phi_kU%val(:, i,j,rcID) + xyz_Phi_k%val(:, i,j,rcID) ) &
            &   - xyr_W_kL%val(1, i,j,rcID) * ( xyz_Phi_k%val(:, i,j,rcID) + xyz_Phi_rijkL(1:rank) ) &
            & ) / ( 2.0d0 * DelXi_k * tfcoordz%G_1div2%val(1, i,j,rcID) )
          
      end do
    end do
  end do

  
  !
  !
  
  xyr_W_kU => get_field_layer(xyr_W, KrMax)
  xyr_W_kL => get_field_layer(xyr_W, KrMax-1)
  xyz_Phi_k => get_field_layer(xyz_Phi, KzMax)
  xyz_Phi_kL => get_field_layer(xyz_Phi, KzMax-1)
  DelXiL = Xir(KrMax) - Xiz(KzMax)
  DelXiLL = Xiz(KzMax) - Xiz(KzMax-1)
  
  DelXi_k = Xir(KrMax) - Xir(KrMax-1)
  
  do rcID=rcIDMin, rcIDMax
    !$omp parallel do private(xyz_Phi_rijkU)
    do j=jEMin, jEMax
      do i=iEMin, iEMax

        xyz_Phi_rijkU(1:rank) =   ( xyz_Phi_k%val(:, i,j,rcID) - xyz_Phi_kL%val(:, i,j,rcID) ) * DelXiL / DelXiLL &
          &                     + xyz_Phi_k%val(:, i,j,rcID)  
          
        ret_xyz_VDelVFlux%field2Ds(kzMax)%val(:, i,j,rcID) = &
            & ( &
            &     xyr_W_kU%val(1, i,j,rcID) * ( xyz_Phi_rijkU(1:rank) + xyz_Phi_k%val(:, i,j,rcID) ) &
            &   - xyr_W_kL%val(1, i,j,rcID) * ( xyz_Phi_k%val(:, i,j,rcID) + xyz_Phi_kL%val(:, i,j,rcID) ) &
            & ) / ( 2.0d0 * DelXi_k * tfcoordz%G_1div2%val(1, i,j,rcID) )
        
        
      end do
    end do
  end do
    
end subroutine TF_VDerivateVFlux

end module TF_HVDerivateTool
