!= 上流点探索
!
!= Finding departure point
!
! Authors::   Hiroki KASHIMURA, Yoshiyuki O. TAKAHASHI, Shin-ichi Takehiro
! Version::   
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2013-2026. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module sltt_dp
  !
  != 上流点探索 (for ISPACK3/spml2)
  !
  != Finding Departure Point (for ISPACK3/spml2)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  !== Procedures List
  !
  ! SLTTDPHor     :: 水平方向の上流点探索
  ! SLTTDPVer     :: 鉛直方向の上流点探索
  ! ---------------------     :: ------------
  ! SLTTDPHor     :: Finding DP in Horizontal
  ! SLTTDPVer     :: Finding DP in Vertical
  !
  !== NAMELIST
  !
  ! NAMELIST#
  !
  !== References
  ! * Williamson, D. L., and Rasch, P. J., 1989:
  !   Two-dimensional semi-Lagrangian transport with shape-preserving interpolation.
  !   <i> Mon. Wea. Rev.</i>, <b>117</b>, 102-129.
  

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

  ! 格子点設定
  ! Grid points settings
  !
#ifdef LIB_MPI
  use gridset, only:       &
    &                imax, & ! 経度格子点数.
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数.
                             ! Number of grid points in latitude
    &                kmax, & ! 鉛直層数.
                             ! Number of vertical level
    &          imax_global   ! 経度格子点数 (全球).
                             ! Number of grid points in longitude on whole globe
  use ua_mpi_module_base, only : kc
#else
  use gridset, only:       &
    &                imax, & ! 経度格子点数.
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数.
                             ! Number of grid points in latitude
    &                kmax, & ! 鉛直層数.
                             ! Number of vertical level
    &           kc => kmax,& ! 並列鉛直層数.
                             ! Number of vertical level for parallel compuation
    &          imax_global   ! 経度格子点数 (全球).
                             ! Number of grid points in longitude on whole globe
#endif

  implicit none

  private

  public :: SLTTDPHor
  public :: SLTTDPVer

contains

  !--------------------------------------------------------------------------------------

  subroutine SLTTDPHor(                           &
    & DelTime, x_Lon, y_Lat, y_SinLat, y_CosLat,  & ! (in)
    & jexmin, jexmax,                             & ! (in)
    & x_ExtLon, y_ExtLat, xyz_ExtU, xyz_ExtV,     & ! (in)
    & xyz_DPLon, xyz_DPLat                        & ! (out)
    & )
    ! 水平方向の上流点探索
    ! Finding DP in Horizontal

    use sltt_const, only : PIx2, dtjw, nloop_dp_h


    real(DP), intent(in ) :: DelTime
                              ! $\Delta t$         
    real(DP), intent(in ) :: x_Lon(0:imax-1)
                              ! $\lambda$ longitude 
    real(DP), intent(in ) :: y_Lat(1:jmax)
                              ! $\varphi$ latitude
    real(DP), intent(in ) :: y_SinLat(1:jmax)
                              ! $\sin\varphi$
    real(DP), intent(in ) :: y_CosLat(1:jmax)
                              ! $\cos\varphi$
    integer , intent(in ) :: jexmin
    integer , intent(in ) :: jexmax
    real(DP), intent(in ) :: x_ExtLon(0:imax_global-1)
                              ! 経度の拡張配列
                              ! Extended array of Lon
    real(DP), intent(in ) :: y_ExtLat(jexmin:jexmax)
                              ! 緯度の拡張配列
                              ! Extended array of Lat    
    real(DP), intent(in ) :: xyz_ExtU(0:imax_global-1, jexmin:jexmax, 1:kc)    
                              ! 東西風速の拡張配列
                              ! Extended array of zonal velocity
    real(DP), intent(in ) :: xyz_ExtV(0:imax_global-1, jexmin:jexmax, 1:kc)
                              ! 南北風速の拡張配列
                              ! Extended array of meridional velocity
    real(DP), intent(out) :: xyz_DPLon(0:imax_global-1, 1:jmax, 1:kc)
                              ! 上流点の経度
                              ! Lon of departure point
    real(DP), intent(out) :: xyz_DPLat(0:imax_global-1, 1:jmax, 1:kc)
                              ! 上流点の緯度
                              ! Lat of departure point

    !
    ! local variables
    !
    real(DP) :: xyz_MPLon (0:imax_global-1, 1:jmax, 1:kc)
                              ! 中間点の経度
                              ! Lon of mid-point
    real(DP) :: xyz_MPLat (0:imax_global-1, 1:jmax, 1:kc)
                              ! 中間点の緯度
                              ! Lat of mid-point
    real(DP) :: xyz_MPLonN(0:imax_global-1, 1:jmax, 1:kc)
                              ! 中間点の経度（現在推定値）
                              ! Lon of mid-point (temporal estimation)
    real(DP) :: xyz_MPLatN(0:imax_global-1, 1:jmax, 1:kc)
                              ! 中間点の緯度（現在推定値）
                              ! Lat of mid-point (now)
    logical  :: FlagWindInt   ! 初期化フラッグ
                              ! Flag for initialization
    logical  :: xyz_FlagEstDP(0:imax_global-1, 1:jmax, 1:kc)

    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer :: t              ! 推定回数の DO ループ用作業変数
                              ! Work variables for DO loop for DP estimation


    ! 初期化処理
    ! initialization
    !
    xyz_FlagEstDP = .true.
    do k = 1, kc
      do j = 1, jmax
        do i = 0, imax_global-1
          xyz_MPLonN(i,j,k) = x_ExtLon(i)
          xyz_MPLatN(i,j,k) = y_Lat(j)
        end do
      end do
    end do


    FlagWindInt = .false.
    mp_search_loop : do t = 1, nloop_dp_h
      if( mod( t, 5 ) .eq. 0 ) &
        & write( 6, * ) "### Loop for searching departure points : ", t


      xyz_MPLon = xyz_MPLonN
      xyz_MPLat = xyz_MPLatN


      call SLTTDPHorCore(                                                   &
        & DelTime, x_Lon, y_Lat, y_SinLat, y_CosLat,                        & ! (in)
        & jexmin, jexmax,                                                   & ! (in)
        & x_ExtLon, y_ExtLat, xyz_ExtU, xyz_ExtV, xyz_MPLon, xyz_MPLat,     & ! (in)
        & xyz_MPLonN, xyz_MPLatN,                                           & ! (inout)
        & xyz_FlagEstDP, FlagWindInt                                        & ! (in)
        & )
      FlagWindInt = .true.

    end do mp_search_loop

    xyz_FlagEstDP = .true.
    FlagWindInt   = .true.

    call SLTTDPHorCore(                                                   &
      & 2.0_DP * DelTime, x_Lon, y_Lat, y_SinLat, y_CosLat,               & ! (in)
      & jexmin, jexmax,                                                   & ! (in)
      & x_ExtLon, y_ExtLat, xyz_ExtU, xyz_ExtV, xyz_MPLon, xyz_MPLat,     & ! (in)
      & xyz_DPLon, xyz_DPLat,                                             & ! (inout)
      & xyz_FlagEstDP, FlagWindInt                                        & ! (in)
      & )


  end subroutine SLTTDPHor

  !--------------------------------------------------------------------------------------

  subroutine SLTTDPHorCore(                                             &
    & DelTime, x_APLon, y_APLat, y_APSinLat, y_APCosLat,                & ! (in)
    & jexmin, jexmax,                                                   & ! (in)
    & x_ExtLon, y_ExtLat, xyz_ExtU, xyz_ExtV, xyz_MPLon, xyz_MPLat,     & ! (in)
    & xyz_DPLon, xyz_DPLat,                                             & ! (inout)
    & xyz_FlagEstDP, FlagWindInt                                        & ! (in)
    & )
    ! 水平方向の上流点探索のコア部分
    ! Finding DP in Horizontal (Core)

    use constants0 , only : PI
    use constants  , only : RPlanet
    use mpi_wrapper, only : myrank
    use sltt_const , only : PIx2, PIH
    use sltt_lagint, only : &
      & SLTTLagIntCubCalcFactHor , SLTTLagIntCubIntHor


    real(DP), intent(in   ) :: DelTime
                              ! $\Delta t$    
    real(DP), intent(in   ) :: x_APLon(0:imax-1)
                              ! 到着点（グリッド上）の経度
                              ! Lon of arrival point (which is on grid)
    real(DP), intent(in   ) :: y_APLat(1:jmax )
                              ! 到着点（グリッド上）の緯度
                              ! Lat of arrival point (which is on grid)
    real(DP), intent(in   ) :: y_APSinLat(1:jmax)
                              ! 到着点の sin(緯度)
                              ! sin(lat) of arrival point
    real(DP), intent(in   ) :: y_APCosLat(1:jmax)
                              ! 到着点の cos(緯度)
                              ! cos(lat) of arrival point    
    integer , intent(in   ) :: jexmin
    integer , intent(in   ) :: jexmax
    real(DP), intent(in   ) :: x_ExtLon(0:imax_global-1)
                              ! 経度の拡張配列
                              ! Extended array of Lon
    real(DP), intent(in   ) :: y_ExtLat(jexmin:jexmax)
                              ! 緯度の拡張配列
                              ! Extended array of Lat
    real(DP), intent(in   ) :: xyz_ExtU(0:imax_global-1, jexmin:jexmax, 1:kc)
                              ! 東西風速の拡張配列
                              ! Extended array of zonal wind
    real(DP), intent(in   ) :: xyz_ExtV(0:imax_global-1, jexmin:jexmax, 1:kc)
                              ! 南北風速の拡張配列
                              ! Extended array of meridional wind
    real(DP), intent(in   ) :: xyz_MPLon(0:imax_global-1, 1:jmax, 1:kc)
                              ! 中間点の経度
                              ! Lon of mid-point
    real(DP), intent(in   ) :: xyz_MPLat(0:imax_global-1, 1:jmax, 1:kc)
                              ! 中間点の緯度
                              ! Lat of mid-point
    real(DP), intent(inout) :: xyz_DPLon(0:imax_global-1, 1:jmax, 1:kc)
                              ! 上流点の経度
                              ! Lon of departure point    
    real(DP), intent(inout) :: xyz_DPLat(0:imax_global-1, 1:jmax, 1:kc)
                              ! 上流点の緯度
                              ! Lat of departure point    
    logical , intent(in   ) :: xyz_FlagEstDP(0:imax_global-1, 1:jmax, 1:kc)
                              ! 上流点探索のフラグ
                              ! Flag for finding departure point
    logical , intent(in   ) :: FlagWindInt
                              ! 初回フラグ
                              ! Initial flag

    !
    ! local variables
    !
    real(DP) :: xyz_MPU(0:imax_global-1, 1:jmax, 1:kc)
                              ! 中間点の東西風速
                              ! Zonal wind at mid-point
    real(DP) :: xyz_MPV(0:imax_global-1, 1:jmax, 1:kc)
                              ! 中間点の南北風速
                              ! Meridional wind at mid-point
    real(DP) :: MPCosLat      ! 中間点のcos(lat)
                              ! cos(lat) at mid-point

    real(DP) :: xyza_lcifx(0:imax_global-1, 1:jmax, 1:kc, -1:2)
    real(DP) :: xyza_lcify(0:imax_global-1, 1:jmax, 1:kc, -1:2)
    integer :: xyz_ii(0:imax_global-1, 1:jmax, 1:kc)
    integer :: xyz_jj(0:imax_global-1, 1:jmax, 1:kc)
                              ! ラグランジュ補間のための作業変数
                              ! Working variables for Lagrange interpolation


    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数

    real(DP) :: InvRPlanet 

    InvRPlanet = 1.0_DP/RPlanet


    if( FlagWindInt ) then

      ! ラグランジュ３次補間
      ! Cubic interpolation
      call SLTTLagIntCubCalcFactHor(                    &
        & jexmin, jexmax,                               & ! (in)
        & x_ExtLon, y_ExtLat, xyz_MPLon, xyz_MPLat,     & ! (in)
        & xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify        & ! (out)
        & )

#ifdef AXISYMMETRY
      xyza_lcifx = 0.25D0    ! 1 次元ラグランジュ補間計算へ
#endif
      
      ! 中間点の風速推定
      ! estimating wind velocities at mid-point
      !
      call SLTTLagIntCubIntHor(                             &
        & jexmin, jexmax,                                   & ! (in)
        & xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify, xyz_ExtU, & ! (in)
        & xyz_MPU                                           & ! (out)
        & )
      call SLTTLagIntCubIntHor(                             &
        & jexmin, jexmax,                                   & ! (in)
        & xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify, xyz_ExtV, & ! (in)
        & xyz_MPV                                           & ! (out)
        & )
    else

      !
      ! 推定風速の更新
      !
      xyz_MPU(0:imax_global-1,1:jmax,:) = xyz_ExtU(0:imax_global-1,1:jmax,:)
      xyz_MPV(0:imax_global-1,1:jmax,:) = xyz_ExtV(0:imax_global-1,1:jmax,:)

    end if
    
    do k = 1, kc
      do j = 1, jmax
        do i = 0, imax_global-1

          if( xyz_FlagEstDP(i,j,k) ) then

!            MPCosLat = cos( xyz_MPLat(i,j,k) )
!!$            xyz_DPLon(i,j,k) = x_APLon(i) &
            xyz_DPLon(i,j,k) = x_ExtLon(i) &
              & - DelTime * xyz_MPU(i,j,k) / cos( xyz_MPLat(i,j,k) ) * InvRPlanet !/ RPlanet
            xyz_DPLat(i,j,k) = y_APLat(j) &
              & - DelTime * xyz_MPV(i,j,k)* InvRPlanet ! / RPlanet

          end if

        end do
      end do
   end do
   
   ! 普通の緯度経度の範囲に直す
   xyz_DPLon = modulo(xyz_DPLon,PIx2)
    
  end subroutine SLTTDPHorCore

  !--------------------------------------------------------------------------------------

  subroutine SLTTDPVer(      &
    & DelTime, xyr_SigmaDot, &
    & xyz_DPSigma            &
    & )
    ! 鉛直方向の上流点探索
    ! Finding DP in Vertical

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only : r_Sigma, z_Sigma

    use sltt_const, only : nloop_dp_v


    real(DP), intent(in ) :: DelTime
                              ! $\Delta t$
    real(DP), intent(in ) :: xyr_SigmaDot(0:imax-1, 1:jmax, 0:kmax)
                              ! 鉛直流速（SigmaDot）    
    real(DP), intent(out) :: xyz_DPSigma (0:imax-1, 1:jmax, 1:kmax)
                              ! 上流点高度
                              ! Sigma of the departure point        


    !
    ! local variables
    !
    real(DP) :: xyz_MPSigma   (0:imax-1, 1:jmax, 1:kmax)
                              ! 中間点の高度
                              ! Sigma of mid-point
    real(DP) :: xyz_MPSigmaN  (0:imax-1, 1:jmax, 1:kmax)
                              ! 中間点の高度（現在推定値）
                              ! Lat of mid-point (temporal estimation)
    real(DP) :: xyz_MPSigmaDot(0:imax-1, 1:jmax, 1:kmax)
                              ! 中間点のSigmaDot
                              ! SigmaDot of mid-point
    real(DP) :: xyza_lcifz(0:imax-1, 1:jmax, 1:kmax, -1:2)
                              ! ラグランジュ補間のための作業変数
                              ! Work variable for Lagrange interpolation
    integer :: xyz_kk(0:imax-1, 1:jmax, 1:kmax)
                              ! ラグランジュ補間のための作業変数
                              ! Work variable for Lagrange interpolation
    

    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k, kk, k2       ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer :: t              ! 推定回数の DO ループ用作業変数
                              ! Work variables for DO loop for DP estimation


    ! 初期化処理
    ! initialization
    !
    do k = 1, kmax
      xyz_MPSigmaN(:,:,k) = z_Sigma(k)
    end do

    ! 上流点推定のループ
    ! Loop for finding departure point
    mp_search_loop : do t = 1, nloop_dp_v


      xyz_MPSigma = xyz_MPSigmaN


      ! 中間点の鉛直流速をラグランジュ３次補間で求める 
      ! vertical wind velocity at mid-point is estimated by using Lagrange cubic 
      ! interpolation
      !
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1

            !
            ! Routine for dcpam
            !
            ! Departure points, xyz_MPSigma(:,:,k), must be located between 
            ! r_Sigma(kk) > xyz_MPSigma(k) > r_Sigma(kk+1).
            ! Further, 1 <= kk <= kmax-2.
            !

            !
            ! economical method
            !
            if( xyz_MPSigma(i,j,k) > r_Sigma(k) ) then
              k_search_1 : do k2 = k, 1, -1
                if( r_Sigma(k2) > xyz_MPSigma(i,j,k) ) exit k_search_1
              end do k_search_1
              xyz_kk(i,j,k) = k2
            else
              k_search_2 : do k2 = min( k+1, kmax ), kmax
                if( r_Sigma(k2) < xyz_MPSigma(i,j,k) ) exit k_search_2
              end do k_search_2
              xyz_kk(i,j,k) = k2 - 1
            end if
            xyz_kk(i,j,k) = min( max( xyz_kk(i,j,k), 1 ), kmax-2 )

          end do
        end do
      end do
      !
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            kk = xyz_kk(i,j,k)
            xyza_lcifz(i,j,k,-1) =                         &
              &     ( xyz_MPSigma(i,j,k) - r_Sigma(kk  ) ) &
              &   * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+1) ) &
              &   * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+2) ) &
              & / ( ( r_Sigma(kk-1) - r_Sigma(kk  ) )      &
              &   * ( r_Sigma(kk-1) - r_Sigma(kk+1) )      &
              &   * ( r_Sigma(kk-1) - r_Sigma(kk+2) ) )
            xyza_lcifz(i,j,k, 0) =                         &
              &     ( xyz_MPSigma(i,j,k) - r_Sigma(kk-1) ) &
              &   * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+1) ) &
              &   * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+2) ) &
              & / ( ( r_Sigma(kk  ) - r_Sigma(kk-1) )      &
              &   * ( r_Sigma(kk  ) - r_Sigma(kk+1) )      &
              &   * ( r_Sigma(kk  ) - r_Sigma(kk+2) ) )
            xyza_lcifz(i,j,k, 1) =                         &
              &     ( xyz_MPSigma(i,j,k) - r_Sigma(kk-1) ) &
              &   * ( xyz_MPSigma(i,j,k) - r_Sigma(kk  ) ) &
              &   * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+2) ) &
              & / ( ( r_Sigma(kk+1) - r_Sigma(kk-1) )      &
              &   * ( r_Sigma(kk+1) - r_Sigma(kk  ) )      &
              &   * ( r_Sigma(kk+1) - r_Sigma(kk+2) ) )
            xyza_lcifz(i,j,k, 2) =                         &
              &     ( xyz_MPSigma(i,j,k) - r_Sigma(kk-1) ) &
              &   * ( xyz_MPSigma(i,j,k) - r_Sigma(kk  ) ) &
              &   * ( xyz_MPSigma(i,j,k) - r_Sigma(kk+1) ) &
              & / ( ( r_Sigma(kk+2) - r_Sigma(kk-1) )      &
              &   * ( r_Sigma(kk+2) - r_Sigma(kk  ) )      &
              &   * ( r_Sigma(kk+2) - r_Sigma(kk+1) ) )
          end do
        end do
      end do
      !
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            kk = xyz_kk(i,j,k)
            xyz_MPSigmaDot(i,j,k) =                             &
              &   xyza_lcifz(i,j,k,-1) * xyr_SigmaDot(i,j,kk-1) &
              & + xyza_lcifz(i,j,k, 0) * xyr_SigmaDot(i,j,kk  ) &
              & + xyza_lcifz(i,j,k, 1) * xyr_SigmaDot(i,j,kk+1) &
              & + xyza_lcifz(i,j,k, 2) * xyr_SigmaDot(i,j,kk+2)
          end do
        end do
      end do



      do k = 1, kmax
        xyz_MPSigmaN(:,:,k) = z_Sigma(k) - xyz_MPSigmaDot(:,:,k) * DelTime
      end do
      xyz_MPSigmaN = min( max( xyz_MPSigmaN, r_Sigma(kmax) ), r_Sigma(0) )


    end do mp_search_loop

    ! 上流点高度の推定
    ! estimating departure point
    !
    do k = 1, kmax
      xyz_DPSigma(:,:,k) = z_Sigma(k) - xyz_MPSigmaDot(:,:,k) * 2.0_DP * DelTime
    end do


  end subroutine SLTTDPVer


end module sltt_dp