module Utility_GeoCoord

  use dc_types, only: &
    & DP

  use igmcore_math, only: &
    & PI

  use igmcore_coordinate_conversion, only: &
    & orth_to_geo_pos, &
    & geo_to_orth_vec, orth_to_geo_vec

  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_EffSize_Min, get_EffSize_Max, &
    & 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

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

  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D

  use Field_IcGrid3D_Manager, only: &
    & Field_IcGrid3D, &
    & Field_IcGrid3D_Init, Field_IcGrid3D_Final, &
    & get_IcGrid3D, get_field_layer

  implicit none
  private

public :: local_k
public :: VecField_GeoCoordConversion, HVelocity_GeoCoordConversion

contains
subroutine VecField_GeoCoordConversion( &
  & xyz_Vec, xyz_SLon, xyz_SLat, xyz_SV &
  & )

  type(Field_IcGrid3D), intent(in) :: xyz_Vec
  type(Field_IcGrid3D), intent(out) :: xyz_SLon
  type(Field_IcGrid3D), intent(out) :: xyz_SLat
  type(Field_IcGrid3D), intent(out) :: xyz_SV

  type(IcGrid2D_FVM), pointer :: icgrid2D
  type(IcGrid3D_FVM), pointer :: icgrid3Dz
  type(Field_IcGrid2D), pointer :: xyz_Vec_k
  type(Field_IcGrid2D), pointer :: xyz_SLon_k, xyz_SLat_k, xyz_SV_k
  real(DP) :: s_pos(3)
  real(DP) :: s_vec(3)

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

  icgrid3Dz => get_icgrid3D(xyz_Vec)
  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)

  do k=1, kzMax
    xyz_Vec_k => get_field_layer(xyz_Vec, k)
    xyz_SLon_k => get_field_layer(xyz_SLon, k)
    xyz_SLat_k => get_field_layer(xyz_SLat, k)
    xyz_SV_k => get_field_layer(xyz_SV, 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))
          s_vec(:) = orth_to_geo_vec(xyz_Vec_k%val(:, i,j,rcID), s_pos(:))

          xyz_SLon_k%val(1, i,j,rcID) = s_vec(1)
          xyz_SLat_k%val(1, i,j,rcID) = s_vec(2)
          xyz_SV_k%val(1, i,j,rcID) = s_vec(3)

        end do
      end do
    end do
  end do


end subroutine VecField_GeoCoordConversion

subroutine HVelocity_GeoCoordConversion(xyz_Vh, U, V)

  type(Field_IcGrid3D), intent(in) :: xyz_Vh
  type(Field_IcGrid3D), intent(out) :: U
  type(Field_IcGrid3D), intent(out) :: V

  type(Field_IcGrid3D) :: dummyW

  call Field_IcGrid3D_Init(dummyW, get_icgrid3D(xyz_Vh), 'dummyW', 1)

  call VecField_GeoCoordConversion( &
    & xyz_Vh, U, V, dummyW &
    & )

  call Field_IcGrid3D_Final(dummyW)

end subroutine HVelocity_GeoCoordConversion

function local_k(lon, lat) result(k)

  implicit none

  real(DP), intent(in) :: lon
  real(DP), intent(in) :: lat
  real(DP) :: k(3)

  if ( lat == PI / 2.0d0 ) then
    k(:) = (/ 0.0d0, 0.0d0, 1.0d0 /)
  else if  ( lat == - PI / 2.0d0 ) then
    k(:) = (/ 0.0d0, 0.0d0, -1.0d0 /)
  else
    k(:) = geo_to_orth_vec( (/ 0.0d0, 0.0d0, 1.0d0 /), (/ lon, lat, 1.0d0 /) )
  end if

end function local_k

end module Utility_GeoCoord
