module Vertical_Interpolation

  use dc_types, only: &
    & DP

  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_EffSize_Min, get_EffSize_Max, &
    & RC_REGIONS_NUM

  use IcGrid3D_FVM_Manager, only: &
    & IcGrid3D_FVM, &
    & get_vertical_level_num, get_icgrid2D

  use Field_Grid1D_Manager, only: &
    & Field_Grid1D

  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D

  use Field_IcGrid3D_Manager, only: &
    & Field_IcGrid3D, &
    & get_field3D_rank, &
    & get_icgrid3D, get_field_layer
	
  implicit none
  private

  public :: interpolate_fieldr_to_z, interpolate_fieldz_to_r

  interface interpolate_fieldr_to_z
    module procedure interpolate_field3Dr_to_z
  end interface interpolate_fieldr_to_z

  interface interpolate_fieldz_to_r
    module procedure interpolate_field1Dz_to_r
    module procedure interpolate_field3Dz_to_r
  end interface interpolate_fieldz_to_r

contains


subroutine interpolate_field1Dz_to_r(field1Dz, field1Dr)
  type(Field_Grid1D), intent(in) :: field1Dz
  type(Field_Grid1D), intent(out) :: field1Dr

  real(DP) :: delXiUU
  real(DP) :: delXiU
  real(DP) :: delXiLL
  real(DP) :: delXiL
  real(DP) :: delXi
  real(DP), pointer :: gridr(:), gridz(:)
  integer :: k, krMax

  gridr => field1Dr%grid1D(:)
  gridz => field1Dz%grid1D(:)
  krMax = size(gridz)

  do k=2, krMax-1
    delXi = gridz(k) - gridz(k-1)
    delXiU = gridz(k) - gridr(k)
    delXiL = gridr(k) - gridz(k-1)

    field1Dr%val(k,:) = ( &
            &   delXiL * field1Dz%val(k,:) + delXiU * field1Dz%val(k-1,:) &
            & ) / delXi

  end do

  delXiU = gridz(1) - gridr(1)
  delXiUU = gridr(2) - gridz(1)

  field1Dr%val(1,:) =   field1Dz%val(1,:) &
          &                          - ( field1Dr%val(2,:) - field1Dz%val(1,:) ) * delXiU / delXiUU

  ! kr=krMax
  delXiL = gridr(krMax) - gridz(krMax-1)
  delXiLL = gridz(krMax-1) - gridr(krMax-1)

  field1Dr%val(krMax,:) = field1Dz%val(krMax-1,:) &
          &                          + ( field1Dz%val(krMax-1,:) - field1Dr%val(krMax-1,:) ) * delXiL / delXiLL

end subroutine interpolate_field1Dz_to_r

subroutine interpolate_field3Dr_to_z(field3Dr, field3Dz)
  type(Field_IcGrid3D), intent(in) :: field3Dr
  type(Field_IcGrid3D), intent(inout) :: field3Dz

  integer :: rcID, i, j, k, rankID
  integer :: EMin, EMax, kzMax, rank
  type(IcGrid3D_FVM), pointer :: icgrid3Dr
  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  type(Field_IcGrid2D), pointer :: field3D_krU
  type(Field_IcGrid2D), pointer :: field3D_krL
  type(Field_IcGrid2D), pointer :: field3D_kz
  real(DP) :: delXiU
  real(DP) :: delXiL
  real(DP) :: delXi

  icgrid3Dz => get_icgrid3D(field3Dz)
  icgrid3Dr => get_icgrid3D(field3Dr)
  kzMax = get_vertical_level_num(icgrid3Dz)
!  EMin = get_EffSize_Min(get_icgrid2D(icgrid3Dz))
!  EMax = get_EffSize_Max(get_icgrid2D(icgrid3Dz))
!  rank = get_field3D_rank(field3Dr)

  !$omp parallel do private(field3D_Kz, field3D_KrL, field3D_KrU, delXi, delXiU, delXiL)
  do k=1, KzMax
    field3D_kz => get_field_layer(field3Dz, k)
    field3D_krL => get_field_layer(field3Dr, k)
    field3D_krU => get_field_layer(field3Dr, k+1)

    delXi = icgrid3Dr%verticalGrid(k+1) - icgrid3Dr%verticalGrid(k)
    delXiU = icgrid3Dr%verticalGrid(k+1) - icgrid3Dz%verticalGrid(k)
    delXiL = icgrid3Dz%verticalGrid(k) - icgrid3Dr%verticalGrid(k)

    field3D_kz%val(:,:,:,:) = ( &
              &   delXiL * field3D_krU%val(:,:,:,:) + delXiU * field3D_krL%val(:,:,:,:) &
              & ) / delXi

  end do

end subroutine interpolate_field3Dr_to_z

subroutine interpolate_field3Dz_to_r(field3Dz, field3Dr)
  type(Field_IcGrid3D), intent(in) :: field3Dz
  type(Field_IcGrid3D), intent(inout) :: field3Dr

  integer :: rcID, i, j, k
  integer :: krMax
  type(IcGrid3D_FVM), pointer :: icgrid3Dr
  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  type(Field_IcGrid2D), pointer :: field3D_kzU
  type(Field_IcGrid2D), pointer :: field3D_krU
  type(Field_IcGrid2D), pointer :: field3D_kzL
  type(Field_IcGrid2D), pointer :: field3D_krL
  type(Field_IcGrid2D), pointer :: field3D_kr
  real(DP) :: delXiUU
  real(DP) :: delXiU
  real(DP) :: delXiLL
  real(DP) :: delXiL
  real(DP) :: delXi

  icgrid3Dz => get_icgrid3D(field3Dz)
  icgrid3Dr => get_icgrid3D(field3Dr)
  krMax = get_vertical_level_num(icgrid3Dr)
!  EMin = get_EffSize_Min(get_icgrid2D(icgrid3Dz))
!  EMax = get_EffSize_Max(get_icgrid2D(icgrid3Dz))

  !$omp parallel do private(field3D_Kr, field3D_KzL, field3D_KzU, delXi, delXiU, delXiL)
  do k=2, KrMax-1
    field3D_kr => get_field_layer(field3Dr, k)
    field3D_kzL => get_field_layer(field3Dz, k-1)
    field3D_kzU => get_field_layer(field3Dz, k)

    delXi = icgrid3Dz%verticalGrid(k) - icgrid3Dz%verticalGrid(k-1)
    delXiU = icgrid3Dz%verticalGrid(k) - icgrid3Dr%verticalGrid(k)
    delXiL = icgrid3Dr%verticalGrid(k) - icgrid3Dz%verticalGrid(k-1)

    field3D_kr%val(:,:,:,:) = ( &
        &   delXiL * field3D_kzU%val(:,:,:,:) + delXiU * field3D_kzL%val(:,:,:,:) &
        & ) / delXi

  end do

  ! kr=1
  field3D_kr => get_field_layer(field3Dr, 1)
  field3D_krU => get_field_layer(field3Dr, 2)
  field3D_kzU => get_field_layer(field3Dz, 1)

  delXiU = icgrid3Dz%verticalGrid(1) - icgrid3Dr%verticalGrid(1)
  delXiUU = icgrid3Dr%verticalGrid(2) - icgrid3Dz%verticalGrid(1)

  field3D_kr%val(:,:,:,:) =   ( 1d0 + delXiU/delXiUU ) * field3D_kzU%val(:,:,:,:) &
    &                       - delXiU/delXiUU * field3D_krU%val(:,:,:,:)

  ! kr=krMax
  field3D_kr => get_field_layer(field3Dr, krMax)
  field3D_kzL => get_field_layer(field3Dz, krMax-1)
  field3D_krL => get_field_layer(field3Dr, krMax-1)

  delXiL = icgrid3Dr%verticalGrid(krMax) - icgrid3Dz%verticalGrid(krMax-1)
  delXiLL = icgrid3Dz%verticalGrid(krMax-1) - icgrid3Dr%verticalGrid(krMax-1)

  field3D_kr%val(:,:,:,:) =   field3D_kzL%val(:,:,:,:) &
    &                          + ( field3D_kzL%val(:,:,:,:) - field3D_krL%val(:,:,:,:) ) * delXiL / delXiLL
  
end subroutine interpolate_field3Dz_to_r


end module Vertical_Interpolation
