!>
!
!
module gridpoint_interpolation

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

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

  ! 数学
  ! Mathematic
  !
  use igmcore_math, only: &
    & PI, toRadians

  ! 宣言文 ; Declaration statement
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: init_gridpoint_interpolation, calc_cache_djk_Wk, interpolation, get_epsilon

  ! 非公開変数
  ! Private variables
  !
  
  !>
  logical :: cache_flag = .false.
  
  !>
  logical :: cache_calc_flag = .false.

  !>
  real(DP) :: epsilon

  !>
  real(DP), allocatable :: cache_djk(:,:)

  !>
  real(DP), allocatable :: cache_Wk(:,:)

  !>
  real(DP), allocatable :: cache_Wtotal(:)

contains

!
!>
!!
subroutine init_gridpoint_interpolation( &
  & lons, lats, del_lon, del_lat, &
  & used_pt_maxnum, interpolation_pts_num, &  ! (in)
  & cache_active )

  ! 宣言文 ; Declaration statement
  !
  real(DP), intent(in) :: lons(:)
  real(DP), intent(in) :: lats(:)
  real(DP), intent(in) :: del_lon
  real(DP), intent(in) :: del_lat
  integer, intent(in), optional :: used_pt_maxnum
  integer, intent(in), optional :: interpolation_pts_num
  logical, intent(in), optional :: cache_active

  ! 実行文 ; Executable statement
  !
  
  epsilon = calc_epsilon(lons(:), lats(:), del_lon, del_lat)
  
  if ( present(cache_active) ) then
    cache_flag = cache_active
    allocate(cache_djk(interpolation_pts_num, used_pt_maxnum))
    allocate(cache_Wk(interpolation_pts_num, used_pt_maxnum))
    allocate(cache_Wtotal(interpolation_pts_num))

  end if

end subroutine init_gridpoint_interpolation

!
!>
!!
subroutine calc_cache_djk_Wk( &
  & k, gp_pt, near_pts, used_pt_num &
  & )

  ! 宣言文 ; Decllaration statement
  !
  integer, intent(in) :: k
  integer, intent(in) :: used_pt_num
  real(DP), intent(in)  :: gp_pt(3)
  real(DP), intent(in) ::  near_pts(used_pt_num, 3)

  ! 作業変数
  ! Work variables
  !
  real(DP) :: S_k(used_pt_num)
  real(DP) :: T_k(used_pt_num)
  integer l

  ! 実行文 ; Executable statement
  !
  
  !
  call calc_d_jk(cache_djk(k,1:used_pt_num), gp_pt(1), gp_pt(2), near_pts(:,1), near_pts(:,2), used_pt_num)
 
  !
  call calc_S_k(S_k(:), cache_djk(k, 1:used_pt_num), used_pt_num)
  call calc_T_k(T_k(:), S_k(:), cache_djk(k, 1:used_pt_num), near_pts(:,:), used_pt_num)
  call calc_W_k(cache_Wk(k, 1:used_pt_num), S_k(:), T_k(:), used_pt_num)
  
  !
  cache_Wtotal(k) = sum(cache_Wk(k,:))

end subroutine calc_cache_djk_Wk

!
!>
!!
function interpolation( &
  & z_k, gp_pt, near_pts, used_pt_num, cache_jID &
  & ) result(z_j)

  ! 宣言文 ; Decllaration statement
  !
  integer, intent(in) :: used_pt_num
  real(DP), intent(in)  :: gp_pt(3)
  real(DP), intent(in) ::  near_pts(used_pt_num, 3)
  real(DP), intent(in) :: z_k(used_pt_num)
  real(DP) z_j
  integer, intent(in), optional :: cache_jID 

  ! 作業変数
  ! Work variables
  !
  integer :: m, l
  real(DP) :: d_jk(used_pt_num)
  real(DP) :: S_k(used_pt_num)
  real(DP) :: T_k(used_pt_num)
  real(DP) :: W_k(used_pt_num), W_total
  real(DP) :: Delz_k(used_pt_num)

  ! 実行文 ; Executable statement
  !
  
  if ( cache_flag ) then
    d_jk(:) = cache_djk(cache_jID, 1:used_pt_num)
  else
    !
    call calc_d_jk(d_jk(:), gp_pt(1), gp_pt(2), near_pts(:,1), near_pts(:,2), used_pt_num)
  end if

  !
  if ( d_jk(1) < epsilon ) then    
    m=1
    do l=2,used_pt_num
      if( d_jk(l) < epsilon) then
        m = m + 1
      else
        exit
      end if
    end do
    
    z_j = sum(z_k(1:m)) / m
    return
  end if

  if( cache_flag  ) then
    !
    W_k(:) = cache_Wk(cache_jID, 1:used_pt_num)
    W_total = cache_Wtotal(cache_jID)

  else
    !
    call calc_S_k(S_k(:), d_jk(:), used_pt_num)
    call calc_T_k(T_k(:), S_k(:), d_jk(:), near_pts(:,:), used_pt_num)
    call calc_W_k(W_k(:), S_k(:), T_k(:), used_pt_num)
    W_total = sum(W_k(:))
  end if
  
  !
  call calc_delz_k(delz_k(:), z_k, W_k(:), d_jk(:), gp_pt, near_pts(:,:), used_pt_num)
  !
  z_j = sum(W_k(:) * ( z_k(:) + Delz_k(:) )) / W_total
  
end function interpolation

function get_epsilon() result(eps)
  real(DP) eps

  eps = epsilon

end function get_epsilon

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

!
!>
!!
function calc_epsilon( &
  & lons, lats, del_lon, del_lat &
  & ) result(eps)
  
  ! 宣言文 ; Declaration statement
  !
  real(DP), intent(in) :: lons(:)
  real(DP), intent(in) :: lats(:)
  real(DP), intent(in) :: del_lon
  real(DP), intent(in) :: del_lat
  real(DP) eps

  ! 作業変数
  ! Work variables
  !
  real(DP) :: max_lat, min_lat

  ! 実行文 ; Executable statement
  !

  max_lat = maxval(lats(:))
  min_lat = minval(lats(:))

  if ( max_lat * min_lat > 0.0d0 ) then
    eps = 0.01d0 * max( 0.5d0 * del_lon * ( dcos(toRadians(max_lat)) + dcos(toRadians(min_lat)) ), del_lat)
  else
    eps = 0.01d0 * max( 0.5d0 * del_lon * ( dcos(dabs(toRadians(max_lat))) + 1.0d0 ), del_lat)    
  end if
  
  !write(*,*) epsilon, max_lat, min_lat, del_lat
end function calc_epsilon

!
!>
!!
subroutine calc_d_jk( &
  & d_jk,     &
  & lambda_j, &
  & theta_j, lambda_k, theta_k, used_pt_num &
  & )

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: used_pt_num
  real(DP), intent(inout) :: d_jk(used_pt_num)
  real(DP), intent(in) :: lambda_j
  real(DP), intent(in) :: theta_j
  real(DP), intent(in) :: lambda_k(used_pt_num)
  real(DP), intent(in) :: theta_k(used_pt_num)

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

  ! 実行文 ; Executable statement
  !

  do k=1, used_pt_num
    d_jk(k) = ds_p1_p2( lambda_j, theta_j, lambda_k(k), theta_k(k) )
  end do

end subroutine calc_d_jk

!
!>
!!
function ds_p1_p2( &
  & lambda1, theta1, lambda2, theta2 &
  & ) result(ds)

  ! 宣言文 ; Declaration statement
  !
  real(DP), intent(in) :: lambda1
  real(DP), intent(in) :: theta1
  real(DP), intent(in) :: lambda2
  real(DP), intent(in) :: theta2
  real(DP) ds

  ! 実行文 ; Executable statement
  !

  ds = dacos( &
       &        dsin(theta1) * dsin(theta2) &
       &      + dcos(theta1) * dcos(theta2) * dcos(lambda1 - lambda2) &
       &   )

end function ds_p1_p2

!
!>
!!
subroutine calc_S_k( &
  & S_k, &
  & d_jk, used_pt_num )

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: used_pt_num
  real(DP), intent(inout) :: S_k(used_pt_num)
  real(DP), intent(in) :: d_jk(used_pt_num)

  ! 作業変数
  ! Work variables
  !
  real(DP) :: r_j
  integer :: k

  ! 実行文 ; Executable statement
  !
  
  r_j = d_jk(used_pt_num)

  do k=1, used_pt_num
    if( d_jk(k) < r_j / 3.0d0 ) then
      S_k(k) = 1.0d0 / d_jk(k)

    else if( r_j < d_jk(k) ) then 
      S_k(k) = 0.0d0

    else 
      S_k(k) = 27.0d0 * ( d_jk(k)/r_j - 1 )**2 / ( 4.0d0 * r_j )

    end if
  end do

end subroutine calc_S_k

!
!>
!!
subroutine calc_T_k( &
  & T_k, &
  & S_k, d_jk, near_pts, used_pt_num &
  & )

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: used_pt_num
  real(DP), intent(inout) :: T_k(used_pt_num)
  real(DP), intent(in) :: S_k(used_pt_num)
  real(DP), intent(in) ::  near_pts(used_pt_num, 3)
  real(DP), intent(in) :: d_jk(used_pt_num)

  ! 作業変数
  ! Work variables
  !
  real(DP) :: tk_sum
  real(DP) ::  d_kl
  real(DP) :: cos_theta_jkl
  integer :: k, l

  ! 実行文 ; Executable statement
  !
  
  do k=1, used_pt_num
    
    tk_sum = 0.0d0

    do l=1, used_pt_num
      if ( k /= l ) then
        d_kl = ds_p1_p2(near_pts(k,1), near_pts(k,2), near_pts(l,1), near_pts(l,2))
        cos_theta_jkl = ( dcos(d_kl) - dcos(d_jk(k)) * dcos(d_jk(l)) ) / ( dsin(d_jk(k)) * dsin(d_jk(l)) )
        tk_sum = tk_sum + S_k(l) * ( 1 - cos_theta_jkl ) 
      end if
    end do

    T_k(k) = tk_sum
  end do
end subroutine calc_T_k

!
!>
!!
subroutine calc_W_k( &
  & W_k, S_k, T_k, used_pt_num &
  & )

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: used_pt_num
  real(DP), intent(inout) :: W_k(used_pt_num)
  real(DP), intent(in) :: S_k(used_pt_num)
  real(DP), intent(in) :: T_k(used_pt_num)

  ! 作業変数
  ! Work variables
  !
  integer :: k, l
  real(DP) :: S_total
  real(DP) :: d_delz_k_dlambda(used_pt_num)

  ! 実行文 ; Executable statement
  !
  
  do k=1, used_pt_num
    
    S_total = 0.0d0
    do l = 1, used_pt_num
      if ( k /= l ) then
        S_total = S_total + S_k(l)
      end if
    end do
    
    W_k(k) = S_k(k)**2 * ( 1.0d0 + T_k(k) / S_total )
  end do

end subroutine calc_W_k

!
!>
!!
subroutine calc_delz_k( &
  & delz_k,               &
  & z_k, W_k, d_jk, gp_pt, near_pts, used_pt_num &
  & )

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: used_pt_num
  real(DP), intent(inout) :: delz_k(used_pt_num)
  real(DP), intent(in) :: z_k(used_pt_num)
  real(DP), intent(in) :: W_k(used_pt_num)
  real(DP), intent(in) :: d_jk(used_pt_num)
  real(DP), intent(in) :: gp_pt(3)
  real(DP), intent(in) :: near_pts(used_pt_num, 3)

  ! 作業変数
  ! Work variables
  !
  real(DP) :: d_delz_k_dlambda(used_pt_num)
  real(DP) :: d_delz_k_dtheta(used_pt_num)
  real(DP) :: W_total
  integer :: k, l
  
  real(DP) :: tmp
  real(DP) :: d_delz_k_dlambda_temp
  real(DP) :: d_delz_k_dtheta_temp
  real(DP) :: d_delz_d_norm(used_pt_num)
  real(DP) :: v

  ! 実行文 ; Executable statement
  !

  do k=1, used_pt_num

    W_total = 0.0d0
    d_delz_k_dlambda_temp = 0.0d0
    d_delz_k_dtheta_temp = 0.0d0

    do l=1, used_pt_num
      if ( l /= k ) then
        W_total = W_total + W_k(l)

        tmp = W_k(l) * ( z_k(l) - z_k(k) ) / ds_p1_p2(near_pts(k,1), near_pts(k,2), near_pts(l,1), near_pts(l,2))**2

        d_delz_k_dlambda_temp =  d_delz_k_dlambda_temp &
          &                     +  tmp * ds_lambda(near_pts(l,1), near_pts(l,2), near_pts(k,1), near_pts(k,2)) 
 
        d_delz_k_dtheta_temp = d_delz_k_dtheta_temp &
          &                    +  tmp * ds_theta(near_pts(l,1), near_pts(l,2), near_pts(k,1), near_pts(k,2))

      end if  
    end do
    
    d_delz_k_dlambda(k) = d_delz_k_dlambda_temp / W_total
    d_delz_k_dtheta(k) = d_delz_k_dtheta_temp / W_total
    d_delz_d_norm(k) = dsqrt(d_delz_k_dlambda(k)**2 + d_delz_k_dtheta(k)**2)

  end do

  if ( maxval(d_delz_d_norm(:)) < 1.0e-10 ) then
    !
    delz_k(:) = 0.0d0

  else
    v = 0.1d0 * ( maxval(z_k(:)) - minval(z_k(:)) ) / maxval(d_delz_d_norm(:))

    do k=1, used_pt_num

      delz_k(k) = ( d_delz_k_dlambda(k) * ds_lambda(gp_pt(1), gp_pt(2), near_pts(k,1), near_pts(k,2)) &
        &           + d_delz_k_dtheta(k) * ds_theta(gp_pt(1), gp_pt(2), near_pts(k,1), near_pts(k,2)) &
        &         ) * v / ( v + d_jk(k) )
    end do

  end if

end subroutine calc_delz_k

!
!>
!!
function ds_lambda( &
  & p1_lambda, p1_theta, p2_lambda, p2_theta &
  & ) result(val)

  ! 宣言文 ; Declaration statement
  !
  real(DP), intent(in) :: p1_lambda
  real(DP), intent(in) :: p1_theta
  real(DP), intent(in) :: p2_lambda
  real(DP), intent(in) :: p2_theta
  real(DP) val

  ! 実行文 ; Executable statement
  !

  if ( dabs(p1_lambda - p2_lambda) > PI ) then
    val = - ( p1_lambda + p2_lambda ) * dcos(p2_theta)
  else
    val = ( p1_lambda - p2_lambda ) * dcos(p2_theta)
  end if
  
end function ds_lambda

!
!>
!!
function ds_theta( &
  & p1_lambda, p1_theta, p2_lambda, p2_theta &
  & ) result(val)

  ! 宣言文 ; Declaration statement
  !
  real(DP), intent(in) :: p1_lambda
  real(DP), intent(in) :: p1_theta
  real(DP), intent(in) :: p2_lambda
  real(DP), intent(in) :: p2_theta
  real(DP) val

  ! 実行文 ; Executable statement
  !
  
  val = p1_theta - p2_theta

end function ds_theta

end module gridpoint_interpolation
