!> 正二十面体格子上の物理場データを, 地理座標格子上の物理場データへと補間により変換する手続きを提供するモジュール. 
!
!
module geographic_mesh_interpolation

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

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

  ! 種類型パラメータ
  ! Kind type parameter
  !
  use dc_types, only: &
    & DP, & ! 倍精度実数型. Double precision
    & TOKEN

  !
  !
  !
  use gtool_history

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

  ! 数学
  ! Mathematic
  !
  use igmcore_math, only: &
    & PI, ToRadians, ToDegrees

  ! 座標変換
  ! Coordinate conversion
  !
  use igmcore_coordinate_conversion, only: &
    & geo_to_orth_pos, orth_to_geo_pos

  ! 球面三角法
  ! Spherical trigonometry
  !
  use igmcore_spherical_trigonometry, only: geodesic_arc_length
  
  ! 正二十面体格子データを管理するクラス
  ! Manager of icosahedral grid data
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_EffSize_Min, get_EffSize_Max, get_idMin, get_icRadius, get_glevel, &
    & RC_REGIONS_NUM

  ! 正二十面格子上の物理場の管理
  ! Manager of physical field on icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, get_field2D_name

  ! 正二十面体格子データの管理
  ! Manager of Icosahedral grid data
  !
  use IcGrid3D_FVM_Manager, only: &
    & IcGrid3D_FVM, &
    & get_vertical_level_num, get_icgrid2D

  !
  !
  !
  use Field_IcGrid3D_Manager, only: &
    & Field_IcGrid3D, &
    & get_IcGrid3D, get_field_layer

  !* igmconvert ****
  !*

  ! 近傍格子点のキャッシュ情報の管理
  !
  !
  use cache_partition_manager, only: &
    & init_cache_partition_manager, &
    & register_icgridPt, &
    &  get_nearest_pts_candiate, partitions

  ! クイックソート
  ! Quick sort
  !
  use quick_sort, only: &
    & qsort

  !
  !
  !
  use gridpoint_interpolation, only: &
    & init_gridpoint_interpolation, &
    & calc_cache_djk_Wk, interpolation, get_epsilon

  ! OpenMP
  !
  !
  use omp_lib

  ! 宣言文 ; Declartion statement
  !
  implicit none
  private


  ! 公開手続き
  ! Public procedure
  !
  public :: geographic_mesh_interpolation_Init, &
    & interpolate_Field_IcGrid2D, interpolate_Field_IcGrid3D, &
    & get_latitudes, get_longitudes

  ! interface 文 ; interface statements
  !
  !
  interface geographic_mesh_interpolation_Init
    module procedure IcGrid2DGeoMeshInterpolation_Init
    module procedure IcGrid3DGeoMeshInterpolation_Init
  end interface geographic_mesh_interpolation_Init

  ! 非公開変数
  ! Private variables
  !
  
  !> 正二十面体格子の管理を行う構造型 IcGrid2D_FVM の変数へのポインタ.
  type(IcGrid2D_FVM), pointer :: icgrid2D

  !> IcGrid2D_FVM クラスのメンバー変数 rcs_AGrid (正二十面体格子の座標データ)を一次元化した配列.
  !! rcs_AGrid における格子点の余白部分や重なり合う要素については, 一次元かの際に取り除かれる. 
  real(DP), allocatable :: icgrid1D(:,:)

  !> icgrid1D の要素数. 
  integer :: icgrid1D_num

  !> 経度方向の格子点数. 
  integer :: lon_num
  
  !> 緯度方向の格子点数. 
  integer :: lat_num

  !> 経度方向の格子座標データを保持する配列. 
  real(DP), allocatable :: lons(:)

  !> 緯度方向の格子座標データを保持する配列. 
  real(DP), allocatable :: lats(:)

  !> 補間によって得られる地理座標格子上の物理場データを保持する配列. 
  real(DP), pointer :: geo_field2Ds(:,:,:,:)
  integer :: geo_field2D_num

  !> 補間によって得られる 3D 地理座標格子上の物理場データを保持する配列.
  real(DP), pointer :: geo_field3Ds(:,:,:,:,:)
  integer :: geo_field3D_num

  !> 近傍格子点のインデックスを保持する配列. 
  integer, allocatable :: near_ptsID(:,:,:)
  
  !> 補間に用いられる近傍格子点の数. 
  integer :: used_pt_num = 10

  !> 補間時に実際に用いられる格子点数を補間点ごとに格納したした配列. 
  !! 
  integer, allocatable :: actual_used_pt_num(:,:)

  !> １つのパーティション内に登録できる格子点の数. 
  integer :: cp_pt_numMax

  !> 経度方向のパーティション数. 
  integer :: cplonNum

  !> 緯度方向のパーティション数. 
  integer :: cplatNum

  !>
  logical :: Wk_calc_flag

contains

!
!> geographic_mesh_interpolation モジュールの初期化. 
!!
!! @param[in] lonGridNum
!! @param[in] latGridNum
!! @param[in] icgrid_ref 構造型 IcGrid2D_FVM の変数.
!!
subroutine IcGrid2DGeoMeshInterpolation_Init( &
  & lonGridNum, latGridNum, field2D_num, &
  & icgrid2D_ref )

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: lonGridNum
  integer, intent(in) :: latGridNum
  type(IcGrid2D_FVM), target, intent(in) :: icgrid2D_ref
  integer, intent(in) :: field2D_num

  ! 作業変数
  ! Work variables
  !
  integer :: cplonID, cplatID
  integer :: rcID, i, j
  integer :: EMin, EMax, EffSize
  real(DP) :: icgrid_density
  real(DP) :: ic_radius

  ! 実行文 ; Executable statement
  !
  lon_num = lonGridNum
  lat_num = latGridNum
  geo_field2D_num = field2D_num

  icgrid2D => icgrid2D_ref
  EMin = get_EffSize_Min(icgrid2D)
  EMax = get_EffSize_Max(icgrid2D)
  EffSize = EMax - EMin + 1
  ic_radius = get_IcRadius(icgrid2D)

  ! 地理座標格子の座標データを保持する配列のメモリを確保する. 
  allocate(lons(lon_num))
  allocate(lats(lat_num))
  
  ! 補間によって得られる地理座標格子上の物理場データを保持する配列のメモリを確保する. 
  allocate(geo_field2Ds(geo_field2D_num, lon_num, lat_num, 1))
  
  icgrid1D_num = 10 * EffSize**2 - 10 * 2* ( EffSize - 1 ) - 8
  allocate(icgrid1D(icgrid1d_num, 3))

  allocate(near_ptsID(lon_num, lat_num, used_pt_num))
  allocate(actual_used_pt_num(lon_num, lat_num))

  ! 地理座標格子を準備する. 
  call prepair_geo_coordinate()
    
  ! 正二十面体格子の座標データを保持する IcGridField クラスのメンバー変数 rcs_AGrid を, 
  ! 一次元配列化したものを取得し, icgrid1D に格納する. 
  call flatten_icgridOrfield_array(icgrid2D%rcs_AGrid, icgrid1D, get_idMin(icgrid2D))
 
  ! 全領域をいくつかのパーティションに分けることで, 近傍格子点探索の高速化を図る. 
  ! これには, cache_partition_manager モジュールを用いる.  

  ! 補間点に近い格子点を探索する際に用いる, 領域のパーティションをいくつ作成するか決定する. 
  icgrid_density = 10 * ( 2**( get_glevel(icgrid2D)+1 ) + 1 ) / 2.0d0 / ( 2.0d0 * PI )
  cplonNum = int( 2.0d0 * PI * dcos(PI/3.0d0) * icgrid_density / sqrt(real(used_pt_num)) )
  cplatNum = int( sqrt(3.0d0) / 2.0d0 * PI * icgrid_density / sqrt(real(used_pt_num)) )

  !write(*,*) 'icgrid density--', icgrid_density, 'cplonNum:', cplonNum, 'cplatNum:', cplatNum
  
  ! 一つのパーティションが近傍点探索のキャッシュ情報として登録できる格子点数を, 
  ! 1 個の補間点に使用する近傍格子点数の 5 倍に設定する.
  cp_pt_numMax = used_pt_num * 2
  !write(*,*) cplonNum * cplatNum * cp_pt_numMax
  !'cp_pt_numMax=', cp_pt_numMax, cplonNum * cplatNum


  ! chache_partition_manager モジュールの初期化を行う.  
  call init_cache_partition_manager( cplonNum, cplatNum, cp_pt_numMax)

  ! 正二十面体格子の全格子点を, 各パーティションに登録する. 
  call register_icgridPt(icgrid1D, icgrid1D_num)
  !write(*,*) cp_pt_numMax

  ! 全ての補間点に対して近傍点を取得する. 
  call calc_near_icgrid_ptsID()

  ! gridpoint_interpolation モジュールを初期化する.
  call init_gridpoint_interpolation( &
    & lons(:), lats(:), 2.0d0*PI/(lon_num - 1), PI/(lat_num - 1), &
    & used_pt_maxnum = used_pt_num, &
    & interpolation_pts_num = lon_num*(lat_num-2), cache_active = .true. &
    & )

  Wk_calc_flag = .false.

end subroutine IcGrid2DGeoMeshInterpolation_Init

!
!> geographic_mesh_interpolation モジュールの初期化.
!!
!! @param[in] lonGridNum
!! @param[in] latGridNum
!! @param[in] icgrid_ref
!!
subroutine IcGrid3DGeoMeshInterpolation_Init( &
  & lonGridNum, latGridNum, &
  & field2D_num, field3D_num, &
  & icgrid3D_ref )

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: lonGridNum
  integer, intent(in) :: latGridNum
  type(IcGrid3D_FVM), target, intent(in) :: icgrid3D_ref
  integer, intent(in) :: field2D_num
  integer, intent(in) :: field3D_num

  !
  !
  !
  integer :: verticalLv_num

  !
  !

  geo_field3D_num = field3D_num
  verticalLv_num = get_vertical_level_num(icgrid3D_ref)

  ! 補間によって得られる地理座標格子上の物理場データを保持する配列のメモリを確保する.
  allocate(geo_field3Ds(geo_field3D_num, verticalLv_num, lonGridNum, latGridNum, 1))

  call IcGrid2DGeoMeshInterpolation_Init( &
    & lonGridNum, latGridNum, &
    & field2D_num, &
    & get_icgrid2D(icgrid3D_ref) &
    & )

end subroutine IcGrid3DGeoMeshInterpolation_Init

! 
!> 正二十面格子上の 2 次元物理場に対して, 緯度経度座標格子への補間を行う.
!!
!! 実際には, この手続き内で各物理場に対して interpolate_geofield2D を呼び出すことで, 補間は行われる.
!!
!! @param[in] fields
!! @return
!!
function interpolate_Field_IcGrid2D( &
  & field2Ds &
  & ) result(interpolated_fields)

  ! 宣言文 ; Declaration statement
  ! 
  type(Field_IcGrid2D), intent(in) :: field2Ds(:)
  real(DP), pointer :: interpolated_fields(:,:,:,:)

  ! 作業変数
  ! Work variables
  !
  integer :: fieldID

  ! 実行文 ; Executable statement
  !

  do fieldID=1, geo_field2D_num
    write(*,*) 'interpolate:', get_field2D_name(field2Ds(fieldID))
    call interpolate_geofield2D(field2Ds(fieldID), geo_field2Ds(fieldID,:,:,:), icgrid1d_num, used_pt_num)
  end do

  interpolated_fields => geo_field2Ds

end function interpolate_Field_IcGrid2D

!
!>
!!
!! @param[in] fields
!! @param[in] field_num
!!
function interpolate_Field_IcGrid3D(field3Ds) result(interpolated_fields)

  ! 宣言文 ; Declaration statement
  !
  type(Field_IcGrid3D), intent(in) :: field3Ds(:)
  real(DP), pointer :: interpolated_fields(:,:,:,:,:)

  ! 作業変数
  ! Work variables
  !
  type(Field_IcGrid2D), pointer :: field3D_k
  integer :: fieldID
  integer :: verticalLv_num
  integer :: k

  ! 実行文 ; Executable statement
  !

  verticalLv_num = get_vertical_level_num(get_icgrid3D(field3Ds(1)))

  do fieldID=1, geo_field3D_num
    do k=1, verticalLv_num
      field3D_k => get_field_layer(field3Ds(fieldID), k)


      write(*,*) 'interpolate:', get_field2D_name(field3D_k)
      call interpolate_geofield2D(field3D_k, geo_field3Ds(fieldID,k,:,:,:), icgrid1d_num, used_pt_num)
    end do
  end do

  interpolated_fields => geo_field3Ds

end function interpolate_Field_IcGrid3D


!
!>
!!
!! @param[in] field
!! @param[in,out] geo_field
!! @param[in] icgrid1dNum
!! @param[in] usedPtNumMax
!! 
subroutine interpolate_geofield2D( &
  & field, &
  & geo_field, &
  & icgrid1dNum, usedPtNumMax &
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Field_IcGrid2D), intent(in) :: field
  real(DP), intent(inout) :: geo_field(:,:,:)
  integer, intent(in) :: icgrid1dNum
  integer, intent(in) :: usedPtNumMax
  
  ! 作業変数
  ! Work variables
  !
  real(DP) :: field1d(icgrid1dNum,1)
  integer :: lonlat_near_ptsID(usedPtNumMax)
  real(DP) :: s_near_pts(usedPtNumMax,3)

  real(DP) :: s_gp_pt(3)
  real(DP) :: ic_radius
  integer :: usePtNum
  real(DP) :: epsilon

  integer :: latID, lonID
  integer :: used_ptID
  integer :: icgrid1dID

  ! 実行文 ; Executable statement
  !

  ic_radius = get_IcRadius(icgrid2D)
  epsilon = get_epsilon( )

  ! Field_IcGrid2D クラスのオブジェクトが持つ物理場データを保持する配列を一次元化したものを取得する.
  call flatten_icgridOrfield_array(field%val, field1d, get_idMin(icgrid2D))
  
  ! 両極の物理場データを格納する. 
  geo_field(:,1,1) = field1d(1,1)
  geo_field(:,lat_num,1) = field1d(2,1)

  ! 各補間点上において, 近傍格子点の物理場データを用いて, 補間点上の物理場の値を補間する.
 
  !$omp parallel do private(latID, usePtNum, lonlat_near_ptsID, used_ptID, icgrid1dID, s_near_pts, s_gp_pt) 
  do lonID=1, lon_num
    do latID=2, lat_num-1
 
      ! 対象となる補間点が, 実際に利用する近傍格子点の総数を取得する. 
      usePtNum = actual_used_pt_num(lonID, latID)
      ! 正二十面体格子の座標データを保持する配列 icgrid1d において, 対応する近傍格子点の ID を取得する.  
      lonlat_near_ptsID(:) = near_ptsID(lonID, latID, :)
      
      ! 近傍格子点の座標データをデカルト座標から地理座標へ変換する. 
      !
      do used_ptID=1, usePtNum
        icgrid1dID = lonlat_near_ptsID(used_ptID)
        
        if ( icgrid1dID == 1 ) then
          s_near_pts(used_ptID, :) = (/ 90d0, 0d0, ic_radius /)
        else if ( icgrid1dID == 2 ) then
          s_near_pts(used_ptID, :) = (/ -90d0, 0d0, ic_radius /)
        else
          s_near_pts(used_ptID, :) = orth_to_geo_pos( icgrid1d(icgrid1dID,:) )
        end if
      end do
      
      ! 補間点の座標データを地理座標として取得する. 
      s_gp_pt(:) = (/ ToRadians(lons(lonID)), ToRadians(lats(latID)), ic_radius /)
      
      ! interpolation を呼び出す前(初回)に, 補間に必要な係数をあらかじめ計算しておく.
      if ( .not. Wk_calc_flag ) then
           call calc_cache_djk_Wk( (latID-2) * lon_num + lonID, s_gp_pt(:), s_near_pts(1:usePtNum,:), usePtNum)
       end if

      ! 近傍格子点における物理場データを用いて, 補間点上の物理場の値を補間する. 
      geo_field(lonID, latID, 1) = &
        & interpolation( &
        &   field1d(lonlat_near_ptsID(1:usePtNum),:), s_gp_pt(:), s_near_pts(1:usePtNum,:), usePtNum, &
        &   cache_jID=(latID-2) * lon_num + lonID )

      !!!!!!!!!!!!!!!!!!!!!!1
      ! f2003 code isNan
      !====================
      if ( isNan(geo_field(lonID,latID,1)) ) then
         write(*,*) '-<Nan report>-----------------', 'lon=', int(lons(lonID)), ',lat=', int(lats(latID)), usePtNum
         do used_ptID=1,usePtNum
           write(*,*) used_ptID, &
            & geodesic_arc_length(ic_radius, geo_to_orth_pos(s_gp_pt), geo_to_orth_pos(s_near_pts(used_ptID,:))), &
            & field1d(lonlat_near_ptsID(used_ptID),:)
         end do
      end if

      !if ( geo_field(lonID,latID,1) > 10000d00 ) then
      !   write(*,*) '-<Error report>-----------------', geo_field(lonID,latID,1)
      !   write(*,*) 'lon=', int(lons(lonID)), 'lat=', int(lats(latID)), usePtNum
      !   do used_ptID=1,usePtNum
      !     write(*,*) used_ptID, geodesic_arc_length(ic_radius, geo_to_orth_pos(s_gp_pt), geo_to_orth_pos(s_near_pts(used_ptID,:))), field1d(lonlat_near_ptsID(used_ptID),:)
      !   end do
       !
      !end if
  
     end do
  end do
  
  ! 2 回目以降の補間には, 1 回目に計算した補間係数をキャシュとして利用する.
  if ( .not. Wk_calc_flag ) then
    Wk_calc_flag = .true.
  end if

end subroutine interpolate_geofield2D

!
!> 緯度経度座標格子の経度方向の格子座標を取得する.
!!
!! @return 経度格子座標を保持した配列.
!!
function get_longitudes() result(longitudes)

  ! 宣言文 ; Declartion statements
  !

  real(DP) :: longitudes(lon_num)

  ! 実行文 ; Executable statement
  !

  longitudes = lons

end function get_longitudes

!
!> 緯度経度座標格子の緯度方向の格子座標を取得する.
!!
!! @return 緯度格子座標を保持した配列.
!!
function get_latitudes() result(latitudes)

  ! 宣言文 ; Declartion statements
  !
  real(DP) :: latitudes(lat_num)

  ! 実行文 ; Executable statement
  !

  latitudes = lats

end function get_latitudes

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 非公開手続き
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!> 地理座標格子の座標データを保持する配列を初期化する. 
!!
!! 
subroutine prepair_geo_coordinate()
  
  ! 作業変数
  ! Work variables
  !
  integer :: lonID
  integer :: latID
  real(DP) :: del_lat
  real(DP) :: del_lon

  ! 実行文 ; Executable statement
  !

  ! 格子の緯度間隔と経度間隔を求める. 
  del_lat = 180d0 / ( lat_num - 1 )
  del_lon = 360d0 / ( lon_num - 1 )

  ! 各格子点の緯度を求める. 
  lats(1) = 90d0
  do latID=2,lat_num
    lats(latID) = lats(latID-1) - del_lat
  end do
  
  ! 各格子点の経度を求める. 
  lons(1) = -180d0
  do lonID=2,lon_num
    lons(lonID) = lons(lonID-1) + del_lon
  end do
  
end subroutine prepair_geo_coordinate


!
!> 構造型 IcGrid2D_FVM の成分である配列 rcs_AGrid (正二十面体格子の座標データを保持),
!! あるいは構造型 Field_IcGrid2D の成分である配列 val (物理場データを保持)を,
!! 一次元化し結果を引数 flat_array に格納する. 
!!
!! @param[in]    icgridOrfield
!! @param[inout] flat_array
!! @param[in]    idMin
!!
subroutine flatten_icgridOrfield_array( &
  & icgridOrfield, &
  & flat_array,    &
  & idMin &
  & )

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: idMin
  real(DP), intent(in) :: icgridOrfield(:,idMin:,idMin:,:)
  real(DP), intent(inout) :: flat_array(:,:)

  ! 作業変数 
  ! Work variables
  !
  integer :: rcID,i,j
  integer :: EMin, EMax 
  integer :: flatArrayID

  ! 実行文 ; Executable statement
  !

  EMin = get_EffSize_Min(icgrid2D)
  EMax = get_EffSize_Max(icgrid2D)

  ! 一次元配列の一, 二要素目はそれぞれ北極と南極の配列要素を格納する. 
  flat_array(1,:) = icgridOrfield(:,EMin,EMin,1)
  flat_array(2,:) = icgridOrfield(:,EMax,EMax,10)
  
  ! 三要素目以降は, 残りの配列要素を順に格納する. 
  flatArrayID = 3
  do rcID=1, RC_REGIONS_NUM
    do i=EMin+1,EMax
      do j=EMin,EMax-1
        flat_array(flatArrayID, :) = icgridOrfield(:,i,j,rcID)
        flatArrayID = flatArrayID + 1
      end do
    end do
  end do

end subroutine flatten_icgridOrfield_array

!
!>
!
subroutine calc_near_icgrid_ptsID()

  ! 作業変数
  ! Work variables
  !
  integer, allocatable :: nearest_ptsID_candiate(:)
  real(DP), allocatable :: nearest_pts(:,:)
  real(DP), allocatable :: nearest_distances(:)
  integer, allocatable :: sorted_indecies(:)
  
  integer :: candiate_num
  integer :: lonID, latID
  integer :: candiateID
  real(DP) :: ic_radius
  real(DP) :: targetPts(3)
  integer :: cut_id
  real(DP) :: cpDellat

  ! 実行文 ; Executable statement
  !

  allocate(nearest_ptsID_candiate(cp_pt_numMax * 9))
  allocate(nearest_pts(cp_pt_numMax, 3))
  allocate(sorted_indecies(cp_pt_numMax * 9))
  allocate(nearest_distances(cp_pt_numMax * 9))

  ic_radius = get_IcRadius(icgrid2D)

  do lonID=1, lon_num
    do latID=2, lat_num-1
               
      nearest_ptsID_candiate(:) =  -1
      call get_nearest_pts_candiate(lons(lonID), lats(latID), nearest_ptsID_candiate, candiate_num)
    
    
      do candiateID=1, candiate_num
       targetPts(:) = geo_to_orth_pos( (/ ToRadians(lons(lonID)), ToRadians(lats(latID)), ic_radius /) )

       nearest_distances(candiateID) = &
          & geodesic_arc_length(ic_radius, targetPts(:), icgrid1D(nearest_ptsID_candiate(candiateID), :))
       
       targetPts(:) = orth_to_geo_pos(icgrid1D(nearest_ptsID_candiate(candiateID), :))
     end do
     
     call qsort( nearest_distances(1:candiate_num), sorted_indecies(1:candiate_num), candiate_num)
     
       
     if ( candiate_num < used_pt_num ) then
       cut_id = candiate_num
     else
       cut_id = used_pt_num
     end if

     near_ptsID(lonID, latID, 1:cut_id) = nearest_ptsID_candiate(sorted_indecies(1:cut_id))
     actual_used_pt_num(lonID, latID) = cut_id
    
    end do
  end do

end subroutine calc_near_icgrid_ptsID

end module geographic_mesh_interpolation
