module BoundaryCondition_Configurator

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

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

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

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

  use igmcore_linear_algebra, only: &
    & vec_length, vec_normalize, cross

  use igmcore_coordinate_conversion, only: &
    & orth_to_geo_pos
    
  ! 2 次元正二十面体格子の管理
  ! Manager of 2D icosahedral grid
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_LocalProcRectInfo, get_IcRadius

  ! 正二十面格子における広域矩形領域のサブドメインである局所矩形領域(各 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
  !
  use IcGrid3D_FVM_Manager, only: &
    & IcGrid3D_FVM, &
    & get_vertical_level_num, get_icgrid2D

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

  ! 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


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

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

  use Utility_GeoCoord, only: &
    & local_k
    
  implicit none
  private

  public :: apply_UpperLowerBoundary_SlipCondition, apply_UpperLowerBoundary_PressCondition

contains

subroutine apply_UpperLowerBoundary_SlipCondition( &
  & xyz_Vh, xyr_W, xyz_E, xyz_R, xyz_P, &
  & xyz_RhoBasic, &
  & tfcoordz, Rd, Cv &
  & )
  
  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: xyz_Vh
  type(Field_IcGrid3D), intent(inout) :: xyr_W
  type(Field_IcGrid3D), intent(inout) :: xyz_E
  type(Field_IcGrid3D), intent(inout) :: xyz_R
  type(Field_IcGrid3D), intent(inout) :: xyz_P
  type(Field_IcGrid3D), intent(inout) :: xyz_RhoBasic
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  real(DP), intent(in) :: Rd
  real(DP), intent(in) :: Cv

  ! 作業変数
  ! Local variables
  !
  integer :: rcID, i, j
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo

  real(DP), pointer :: Xiz(:), Xir(:)
  type(IcGrid3D_FVM), pointer :: icgrid3Dz, icgrid3Dr
  type(Field_icGrid3D) :: xyz_Rho

  real(DP) :: Grad_xyk_hs_rij(3), xyz_W_rijk(3), xyr_VH_rijk(3), ps_rijk
  
  integer :: k, KzMax, KrMax
  
  ! 格子に関する情報を取得する.
  ! Get the information about grid.

  icgrid3Dz => get_icgrid3D(xyz_Vh)
  icgrid3Dr => get_icgrid3D(xyr_W)

  lcRCInfo => get_LocalProcRectInfo(get_icgrid2D(icgrid3Dr))

  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)

  call Field_IcGrid3D_Init(xyz_Rho, icgrid3Dz, 'xyz_Rho', 1)
  
  Xiz => icgrid3Dz%verticalGrid(:)
  Xir => icgrid3Dr%verticalGrid(:)
  
  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jEMax
      do i=iEMin, iEMax
         
         !
         !
         
         ! 
        do k=1, 3
          Grad_xyk_hs_rij(:) = - tfcoordz%G_1div2%val(1, i,j,rcID) * tfcoordz%G_Z%field2Ds(1)%val(1:3, i,j,rcID)          
          xyz_W_rijk(k) = dot_product(xyz_Vh%field2Ds(k)%val(:,i,j,rcID), Grad_xyk_hs_rij(:))
        end do 
        
        xyr_W%field2Ds(1)%val(1, i,j,rcID) = ( 3.0d0 * xyz_W_rijk(1) - xyz_W_rijk(2) ) / 2d0
        !xyr_W%field2Ds(KrMax)%val(1, i,j,rcID) = xyr_W%field2Ds(KrMax-1)%val(1, i,j,rcID) 
        
        !xyz_P%field2Ds(1)%val(1, i,j,rcID) = xyz_P%field2Ds(2)%val(1, i,j,rcID)
        !xyz_E%field2Ds(1)%val(1, i,j,rcID) = xyz_E%field2Ds(2)%val(1, i,j,rcID)
        
      end do
    end do
  end do

  call Field_IcGrid3D_Final(xyz_Rho)
  
end subroutine apply_UpperLowerBoundary_SlipCondition

subroutine apply_UpperLowerBoundary_PressCondition( &
  & xyz_P, xyz_E, xyr_W, xyr_GW, xyz_R, &
  & tfcoordz, Grav, Rd, Cv &
  & )
  
  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid3D), intent(inout) :: xyz_P
  type(Field_IcGrid3D), intent(inout) :: xyz_E
  type(Field_IcGrid3D), intent(in) :: xyr_W
  type(Field_IcGrid3D), intent(in) :: xyr_GW
  type(Field_IcGrid3D), intent(in) :: xyz_R
  type(TFCoord_IcGrid3D), intent(in) :: tfcoordz
  real(DP), intent(in) :: Grav
  real(DP), intent(in) :: Rd
  real(DP), intent(in) :: Cv

  ! 作業変数
  ! Local variables
  !
  integer :: rcID, i, j
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  type(LocalProcRectInfo), pointer :: lcRCInfo

  real(DP), pointer :: Xiz(:), Xir(:)
  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(IcGrid3D_FVM), pointer :: icgrid3Dz, icgrid3Dr
  type(Field_icGrid3D) :: xyz_Rho

  real(DP) :: Grad_xyk_hs_rij(3), xyz_W_rijk(3), xyr_VH_rijk(3), xyz_ps
  real(DP) :: DelzP, verticalUniVec(3), s_pos(3), xyz_GW
  
  integer :: k, KzMax, KrMax
  
  ! 格子に関する情報を取得する.
  ! Get the information about grid.

  icgrid3Dz => get_icgrid3D(xyz_P)
  icgrid3Dr => get_icgrid3D(xyr_W)
  icgrid2D => get_icgrid2D(icgrid3Dz)

  lcRCInfo => get_LocalProcRectInfo(get_icgrid2D(icgrid3Dr))

  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)
  Xiz => icgrid3Dz%verticalGrid(:)
  
  do rcID=rcIDMin, rcIDMax
    do j=jEMin, jEMax
      do i=iEMin, iEMax
        s_pos(:) = orth_to_geo_pos(icgrid2D%rcs_AGrid(:,i,j,rcID))
        verticalUnivec(:) = local_k(s_pos(1), s_pos(2))
        
        xyz_GW = 0.5d0 * ( xyr_GW%field2Ds(1)%val(1,i,j,rcID) + xyr_GW%field2Ds(2)%val(1,i,j,rcID) )       
        DelzP = - tfcoordz%G_1div2%val(1,i,j,rcID) * ( xyz_R%field2Ds(1)%val(1, i,j,rcID) * Grav + xyz_GW )
         
        xyz_ps = Rd / Cv * xyz_E%field2Ds(1)%val(1, i,j,rcID) - xyz_P%field2Ds(1)%val(1, i,j,rcID)
      !  xyz_P%field2Ds(1)%val(1, i,j,rcID) = ( 4d0 * xyz_P%field2Ds(2)%val(1, i,j,rcID) & 
      !    & -  xyz_P%field2Ds(3)%val(1, i,j,rcID) - DelzP * 2d0 * ( Xiz(2) - Xiz(1) ) ) / 3d0  
      !  xyz_E%field2Ds(1)%val(1, i,j,rcID) = Cv / Rd * ( xyz_ps + xyz_P%field2Ds(1)%val(1, i,j,rcID) ) 


        xyz_GW = 0.5d0 * ( xyr_GW%field2Ds(krMax-1)%val(1,i,j,rcID) + xyr_GW%field2Ds(krMax)%val(1,i,j,rcID) )       
        DelzP = - tfcoordz%G_1div2%val(1,i,j,rcID) * ( xyz_R%field2Ds(KzMax)%val(1, i,j,rcID) * Grav + xyz_GW )
         
        xyz_ps = Rd / Cv * xyz_E%field2Ds(KzMax)%val(1, i,j,rcID) - xyz_P%field2Ds(kzMax)%val(1, i,j,rcID)
!        xyz_P%field2Ds(KzMax)%val(1, i,j,rcID) = xyz_P%field2Ds(KzMax-1)%val(1, i,j,rcID) + DelzP * ( Xiz(KzMax) - Xiz(KzMax-1) ) 
!        xyz_E%field2Ds(KzMax)%val(1, i,j,rcID) = Cv / Rd * ( xyz_ps + xyz_P%field2Ds(KzMax)%val(1, i,j,rcID) ) 


      end do
    end do
  end do

end subroutine apply_UpperLowerBoundary_PressCondition


end module BoundaryCondition_Configurator
