!= 物質移流計算 (セミラグランジュ法）で用いる補間法
!
!= Interpolation methods for Semi-Lagrangian method
!
! Authors::   Hiroki KASHIMURA, Yoshiyuki O. TAKAHASHI
! Version::   
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2013. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module sltt_lagint
  !
  != セミラグランジュ法 で用いる補間法
  !
  != Interpolation methods for Semi-Lagrangian method
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! セミラグランジュ法で用いる補間を演算するモジュールです. 
  ! スペクトル変換・高精度補間に由来する人工的な短波を除去するために Sun et al. (1996) の
  ! 単調フィルタを応用したものを部分的に用いている。 
  !
  ! This is a Interpolation module. Semi-Lagrangian method (Enomoto 2008 modified)
  ! Monotonicity filter (Sun et al 1996) is partly used.
  !
  !== Procedures List
  !
  ! SLTTLagIntCubCalcFactHor :: 水平２次元ラグランジュ３次補間用の係数計算
  ! SLTTLagIntCubIntHor      :: 水平２次元ラグランジュ３次補間（上流点探索で用いる）
  ! SLTTLagIntCubCalcFactVer :: 鉛直１次元ラグランジュ３次補間用の係数計算
  ! SLTTLagIntCubIntVer      :: 鉛直１次元ラグランジュ３次補間（上流点探索で用いる）
  ! SLTTIrrHerIntK13         :: 水平２次元変則エルミート５次補間
  ! SLTTIrrHerIntQui2DHor    :: 水平２次元変則エルミート５次補間（コア部分）
  ! SLTTIrrHerIntQui1DUni    :: １次元変則エルミート５次補間（等間隔格子）
  ! SLTTIrrHerIntQui1DNonUni :: １次元変則エルミート５次補間（不等間隔格子）  
  ! SLTTIrrHerIntQui1DUniLon :: １次元変則エルミート５次補間（等間隔：経度方向専用）  
  ! SLTTHerIntCub1D          :: １次元エルミート３次補間
  ! SLTTHerIntCub2D          :: ２次元エルミート３次補間
  ! ---------------------    :: ------------
  ! SLTTLagIntCubCalcFactHor :: Calculation of factors for 2D Lagrangian Cubic Interpolation
  ! SLTTLagIntCubIntHor      :: 2D Lagrangian Cubic Interpolation used in finding DP in horizontal
  ! SLTTLagIntCubCalcFactVer :: Calculation of factors for 1D Lagrangian Cubic Interpolation
  ! SLTTLagIntCubIntVer      :: 1D Lagrangian Cubic Interpolation used in finding DP in vertical
  ! SLTTHerIntK13            :: Horizontal 2D Irregular Hermite Quintic Interpolation 
  ! SLTTIrrHerIntQui2DHor    :: Horizontal 2D Irregular Hermite Quintic Interpolation (Core)
  ! SLTTIrrHerIntQui1DUni    :: 1D Irregular Hermite Quintic Interpolation for uniform grids
  ! SLTTIrrHerIntQui1DNonUni :: 1D Irregular Hermite Quintic Interpolation for non-uniform grids
  ! SLTTIrrHerIntQui1DUniLon :: 1D Irregular Hermite Quintic Interpolation for uniform longitude grids
  ! SLTTHerIntCub1D          :: 1D Hermite Cubic Interpolation
  ! SLTTHerIntCub2D          :: 2D Hermite Cubic Interpolation
  !
  !== NAMELIST
  !
  ! NAMELIST#
  !
  !== References
  ! * Enomoto, T., 2008: 
  !   Bicubic Interpolation with Spectral Derivatives. 
  !   <i>SOLA</i>, <b>4</b>, 5-8. doi:10.2151/sola.2008-002
  !
  ! * Sun, W.-Y., Yeh, K.-S., and Sun, R.-Y., 1996: 
  !   A simple semi-Lagrangian scheme for advection equations. 
  !   <i>Quarterly Journal of the Royal Meteorological Society</i>, 
  !   <b>122(533)</b>, 1211-1226. doi:10.1002/qj.49712253310
  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP,  & ! 倍精度実数型. Double precision.
    &                 TOKEN  ! キーワード.   Keywords. 

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only:       &
    &                imax, & ! 経度格子点数.
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数.
                             ! Number of grid points in latitude
    &                kmax    ! 鉛直層数.
                             ! Number of vertical level
  use composition, only:                              &
    &                    ncmax,                       &
                             ! 成分の数
                             ! Number of composition
    &                    CompositionInqFlagAdv

  ! 座標データ設定
  ! Axes data settings
  !
  use axesset, only: &
      & DeltaLon, InvDeltaLon ! 経度格子点間隔とその逆数
                              ! Interval of Longitude grids and its inverse

  implicit none

  private

!  public :: SLTTLagIntQuadCalcFactHor
!  public :: SLTTLagIntQuadIntHor
!  public :: SLTTLagIntQuadCalcFactVer
!  public :: SLTTLagIntQuadIntVer
  public :: SLTTLagIntCubCalcFactHor
  public :: SLTTLagIntHorMaxMin
  public :: SLTTLagIntCubIntHor
  public :: SLTTLagIntCubCalcFactVer
  public :: SLTTLagIntCubIntVer
  public :: SLTTIrrHerIntK13
  public :: SLTTIrrLinInt
  public :: SLTTHerIntCub1D
  public :: SLTTHerIntCub2D  
  public :: SLTTIrrHerIntQui1DUni
  public :: SLTTIrrHerIntQui1DNonUni
!  public :: SLTTHerIntQui1D  
!  public :: judgeSun1996


  character(*), parameter:: module_name = 'sltt_lagint'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: sltt_lagint.F90,v 1.4 2013/09/21 14:42:08 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

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

!  subroutine SLTTLagIntQuadCalcFactHor(             &
!    & x_ExtLon, y_ExtLat, xyz_MPLon, xyz_MPLat,     & ! (in)
!    & xyz_ii, xyz_jj, xyza_lqifx, xyza_lqify        & ! (out)
!    & )
!
!    use sltt_const , only : PIx2, jew
!    use mpi_wrapper, only : myrank
!
!    real(DP), intent(in ) :: x_ExtLon(-2+0:imax-1+3)
!    real(DP), intent(in ) :: y_ExtLat(-jew+1:jmax/2+jew)
!    real(DP), intent(in ) :: xyz_MPLon(0:imax-1, 1:jmax/2, 1:kmax)
!    real(DP), intent(in ) :: xyz_MPLat(0:imax-1, 1:jmax/2, 1:kmax)
!    integer , intent(out) :: xyz_ii(0:imax-1, 1:jmax/2, 1:kmax)
!    integer , intent(out) :: xyz_jj(0:imax-1, 1:jmax/2, 1:kmax)
!    real(DP), intent(out) :: xyza_lqifx(0:imax-1, 1:jmax/2, 1:kmax, 0:2)
!    real(DP), intent(out) :: xyza_lqify(0:imax-1, 1:jmax/2, 1:kmax, 0:2)
!
!    !
!    ! local variables
!    !
!    integer :: ii
!    integer :: jj
!
!    integer :: i, j, k, j2
!
!
!    xyz_ii = int( xyz_MPLon / ( PIx2 / imax ) )
!
!    do k = 1, kmax
!      do j = 1, jmax/2
!        do i = 0, imax-1
!
!          if( xyz_MPLat(i,j,k) > y_ExtLat(j) ) then
!            j_search_1 : do j2 = j+1, jmax/2+jew
!              if( y_ExtLat(j2) > xyz_MPLat(i,j,k) ) exit j_search_1
!            end do j_search_1
!            xyz_jj(i,j,k) = j2 - 1
!            if ( xyz_jj(i,j,k) > jmax/2+jew-2 ) xyz_jj(i,j,k) = jmax/2+jew-2
!          else
!            j_search_2 : do j2 = j-1, -jew+1, -1
!              if( y_ExtLat(j2) < xyz_MPLat(i,j,k) ) exit j_search_2
!            end do j_search_2
!            xyz_jj(i,j,k) = j2
!            if ( xyz_jj(i,j,k) < -jew+1 ) xyz_jj(i,j,k) = -jew+1
!          end if
!
!        end do
!      end do
!    end do
!
!#ifdef SLTT_CHECK
!    do k = 1, kmax
!      do j = 1, jmax/2
!        do i = 0, imax-1
!          if( ( xyz_jj(i,j,k) < 0 ) .or. ( xyz_jj(i,j,k) > (jmax+1) ) ) then
!            write( 6, * ) 'Error: in sltt_dp_h0 : ', &
!              'Latitudinal array size is not enough.'
!            write( 6, * ) ' myrank = ', myrank, ' i = ', i, ' j = ', j, ' k = ', k
!            write( 6, * ) 'jj = ', xyz_jj(i,j,k), ' jmax = ', jmax
!            stop
!          end if
!        end do
!      end do
!    end do
!#endif
!
!    !
!    ! calculation of Lagrange cubic interpolation factor 
!    ! for longitudinal direction
!    !
!    do k = 1, kmax
!      do j = 1, jmax/2
!        do i = 0, imax-1
!          ii = xyz_ii(i,j,k)
!          jj = xyz_jj(i,j,k)
!          xyza_lqifx(i,j,k,0) =                             &
!            &     ( xyz_MPLon(i,j,k) - x_ExtLon(ii+1) )     &
!            &   * ( xyz_MPLon(i,j,k) - x_ExtLon(ii+2) )     &
!            & / ( ( x_ExtLon(ii  )   - x_ExtLon(ii+1) )     &
!            &   * ( x_ExtLon(ii  )   - x_ExtLon(ii+2) ) )
!          xyza_lqifx(i,j,k,1) =                             &
!            &     ( xyz_MPLon(i,j,k) - x_ExtLon(ii  ) )     &
!            &   * ( xyz_MPLon(i,j,k) - x_ExtLon(ii+2) )     &
!            & / ( ( x_ExtLon(ii+1)   - x_ExtLon(ii  ) )     &
!            &   * ( x_ExtLon(ii+1)   - x_ExtLon(ii+2) ) )
!          xyza_lqifx(i,j,k,2) =                             &
!            &     ( xyz_MPLon(i,j,k) - x_ExtLon(ii  ) )     &
!            &   * ( xyz_MPLon(i,j,k) - x_ExtLon(ii+1) )     &
!            & / ( ( x_ExtLon(ii+2)   - x_ExtLon(ii  ) )     &
!            &   * ( x_ExtLon(ii+2)   - x_ExtLon(ii+1) ) )
!
!          xyza_lqify(i,j,k,0) =                             &
!            &     ( xyz_MPLat(i,j,k) - y_ExtLat(jj+1) )     &
!            &   * ( xyz_MPLat(i,j,k) - y_ExtLat(jj+2) )     &
!            & / ( ( y_ExtLat(jj  )   - y_ExtLat(jj+1) )     &
!            &   * ( y_ExtLat(jj  )   - y_ExtLat(jj+2) ) )
!          xyza_lqify(i,j,k,1) =                             &
!            &     ( xyz_MPLat(i,j,k) - y_ExtLat(jj  ) )     &
!            &   * ( xyz_MPLat(i,j,k) - y_ExtLat(jj+2) )     &
!            & / ( ( y_ExtLat(jj+1)   - y_ExtLat(jj  ) )     &
!            &   * ( y_ExtLat(jj+1)   - y_ExtLat(jj+2) ) )
!          xyza_lqify(i,j,k,2) =                             &
!            &     ( xyz_MPLat(i,j,k) - y_ExtLat(jj  ) )     &
!            &   * ( xyz_MPLat(i,j,k) - y_ExtLat(jj+1) )     &
!            & / ( ( y_ExtLat(jj+2)   - y_ExtLat(jj) )       &
!            &   * ( y_ExtLat(jj+2)   - y_ExtLat(jj+1) ) )
!        end do
!      end do
!    end do
!
!
!  end subroutine SLTTLagIntQuadCalcFactHor

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

!  subroutine SLTTLagIntQuadIntHor(                        &
!    & xyz_ii, xyz_jj, xyza_lqifx, xyza_lqify, xyz_ExtArr, & ! (in)
!    & xyz_MPArr                                           & ! (out)
!    & )
!
!    use sltt_const , only : jew
!
!    integer , intent(in ) :: xyz_ii(0:imax-1, 1:jmax/2, 1:kmax)
!    integer , intent(in ) :: xyz_jj(0:imax-1, 1:jmax/2, 1:kmax)
!    real(DP), intent(in ) :: xyza_lqifx(   0:imax-1  ,      1:jmax/2    , 1:kmax, 0:2)
!    real(DP), intent(in ) :: xyza_lqify(   0:imax-1  ,      1:jmax/2    , 1:kmax, 0:2)
!    real(DP), intent(in ) :: xyz_ExtArr(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax)
!    real(DP), intent(out) :: xyz_MPArr (   0:imax-1  ,      1:jmax/2    , 1:kmax)
!
!
!    !
!    ! local variables
!    !
!    integer :: ii
!    integer :: jj
!    integer :: kk
!
!    integer :: i, j, k
!
!
!    do k = 1, kmax
!      do j = 1, jmax/2
!        do i = 0, imax-1
!          ii = xyz_ii(i,j,k)
!          jj = xyz_jj(i,j,k)
!          kk = k
!          xyz_MPArr(i,j,k) =                                      &
!            &   xyza_lqify(i,j,k,0)                               &
!            & * ( xyza_lqifx(i,j,k,0) * xyz_ExtArr(ii  ,jj  ,kk)  &
!            &   + xyza_lqifx(i,j,k,1) * xyz_ExtArr(ii+1,jj  ,kk)  &
!            &   + xyza_lqifx(i,j,k,2) * xyz_ExtArr(ii+2,jj  ,kk) )&
!            & + xyza_lqify(i,j,k,1)                               &
!            & * ( xyza_lqifx(i,j,k,0) * xyz_ExtArr(ii  ,jj+1,kk)  &
!            &   + xyza_lqifx(i,j,k,1) * xyz_ExtArr(ii+1,jj+1,kk)  &
!            &   + xyza_lqifx(i,j,k,2) * xyz_ExtArr(ii+2,jj+1,kk) )&
!            & + xyza_lqify(i,j,k,2)                               &
!            & * ( xyza_lqifx(i,j,k,0) * xyz_ExtArr(ii  ,jj+2,kk)  &
!            &   + xyza_lqifx(i,j,k,1) * xyz_ExtArr(ii+1,jj+2,kk)  &
!            &   + xyza_lqifx(i,j,k,2) * xyz_ExtArr(ii+2,jj+2,kk) )
!        end do
!      end do
!    end do
!
!
!  end subroutine SLTTLagIntQuadIntHor
!
!  !--------------------------------------------------------------------------------------
!
!  subroutine SLTTLagIntQuadCalcFactVer(  &
!    & xyz_DPSigma,                & ! (in)
!    & xyza_lqifz, xyz_kk          & ! (out)
!    & )
!
!    ! 座標データ設定
!    ! Axes data settings
!    !
!    use axesset, only : z_Sigma
!
!
!    real(DP), intent(in ) :: xyz_DPSigma (0:imax-1, 1:jmax, 1:kmax)
!    real(DP), intent(out) :: xyza_lqifa_z(0:imax-1, 1:jmax, 1:kmax, 0:2)
!    integer , intent(out) :: xyz_kk    (0:imax-1, 1:jmax, 1:kmax)
!
!
!    !
!    ! local variables
!    !
!
!    integer :: i
!    integer :: j
!    integer :: k
!    integer :: kk
!    integer :: k2
!
!
!    do k = 1, kmax
!      do j = 1, jmax
!        do i = 0, imax-1
!
!          !
!          ! Routine for dcpam
!          !
!          ! Departure points, xyz_DPSigma(:,:,k), must be located between 
!          ! z_Sigma(kk) > xyz_DPSigma(k) > z_Sigma(kk+1).
!          ! Further, 1 <= kk <= kmax-2.
!          !
!
!          !
!          ! economical method
!          !
!          if( xyz_DPSigma(i,j,k) > z_Sigma(k) ) then
!            k_search_1 : do k2 = k, 1, -1
!              if( z_Sigma(k2) > xyz_DPSigma(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( z_Sigma(k2) < xyz_DPSigma(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_lqifa_z(i,j,k,0) =                          &
!            &     ( xyz_DPSigma(i,j,k) - z_Sigma(kk+1) ) &
!            &   * ( xyz_DPSigma(i,j,k) - z_Sigma(kk+2) ) &
!            & / ( ( z_Sigma(kk  )      - z_Sigma(kk+1) ) &
!            &   * ( z_Sigma(kk  )      - z_Sigma(kk+2) ) )
!          xyza_lqifa_z(i,j,k,1) =                          &
!            &     ( xyz_DPSigma(i,j,k) - z_Sigma(kk  ) ) &
!            &   * ( xyz_DPSigma(i,j,k) - z_Sigma(kk+2) ) &
!            & / ( ( z_Sigma(kk+1)      - z_Sigma(kk  ) ) &
!            &   * ( z_Sigma(kk+1)      - z_Sigma(kk+2) ) )
!          xyza_lqifa_z(i,j,k,2) =                          &
!            &     ( xyz_DPSigma(i,j,k) - z_Sigma(kk  ) ) &
!            &   * ( xyz_DPSigma(i,j,k) - z_Sigma(kk+1) ) &
!            & / ( ( z_Sigma(kk+2)      - z_Sigma(kk  ) ) &
!            &   * ( z_Sigma(kk+2)      - z_Sigma(kk+1) ) )
!        end do
!      end do
!    end do
!
!
!  end subroutine SLTTLagIntQuadCalcFactVer
!
!  !--------------------------------------------------------------------------------------
!
!  subroutine SLTTLagIntQuadIntVer( &
!    & xyz_Arr, xyza_lqifz, xyz_kk, & ! (in)
!    & xyz_ArrA                     & ! (out)
!    & )
!
!
!    real(DP), intent(in ) :: xyz_Arr   (0:imax-1, 1:jmax, 1:kmax)
!    real(DP), intent(in ) :: xyza_lqifa_z(0:imax-1, 1:jmax, 1:kmax, 0:2)
!    integer , intent(in ) :: xyz_kk    (0:imax-1, 1:jmax, 1:kmax)
!    real(DP), intent(out) :: xyz_ArrA  (0:imax-1, 1:jmax, 1:kmax)
!
!
!    !
!    ! local variables
!    !
!
!    integer :: i
!    integer :: j
!    integer :: k
!    integer :: kk
!
!
!    do k = 1, kmax
!      do j = 1, jmax
!        do i = 0, imax-1
!          kk = xyz_kk(i,j,k)
!          xyz_ArrA(i,j,k) =                             &
!            &   xyza_lqifa_z(i,j,k,0) * xyz_Arr(i,j,kk  ) &
!            & + xyza_lqifa_z(i,j,k,1) * xyz_Arr(i,j,kk+1) &
!            & + xyza_lqifa_z(i,j,k,2) * xyz_Arr(i,j,kk+2)
!        end do
!      end do
!    end do
!
!
!  end subroutine SLTTLagIntQuadIntVer

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

  subroutine SLTTLagIntCubCalcFactHor(              &
    & iexmin, iexmax, jexmin, jexmax,               & ! (in)
    & x_ExtLon, y_ExtLat, xyz_MPLon, xyz_MPLat,     & ! (in)
    & xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify        & ! (out)
    & )
  ! 水平２次元ラグランジュ３次補間の係数計算
  ! Calculation of factors for 2D Lagrangian cubic interpolation

    use sltt_const , only : PIx2
    use mpi_wrapper, only : myrank

    integer , intent(in ) :: iexmin
    integer , intent(in ) :: iexmax
    integer , intent(in ) :: jexmin
    integer , intent(in ) :: jexmax
    real(DP), intent(in ) :: x_ExtLon(iexmin:iexmax)
    real(DP), intent(in ) :: y_ExtLat(jexmin:jexmax)
    real(DP), intent(in ) :: xyz_MPLon(0:imax-1, 1:jmax/2, 1:kmax)
    real(DP), intent(in ) :: xyz_MPLat(0:imax-1, 1:jmax/2, 1:kmax)
    integer , intent(out) :: xyz_ii(0:imax-1, 1:jmax/2, 1:kmax)
    integer , intent(out) :: xyz_jj(0:imax-1, 1:jmax/2, 1:kmax)
    real(DP), intent(out) :: xyza_lcifx(0:imax-1, 1:jmax/2, 1:kmax, -1:2)
    real(DP), intent(out) :: xyza_lcify(0:imax-1, 1:jmax/2, 1:kmax, -1:2)

    !
    ! local variables
    !
    integer :: ii
    integer :: jj

    integer :: i, j, k, j2

!    integer  :: ns            ! 南北半球の違いに対応するための変数
!    real(DP) :: in_deltalat
!    
!if (y_ExtLat(jmax/4) > 0) then
!    ns = 0
!else
!    ns = jmax/2 - 1
!endif
!in_deltalat = 1.0_DP/(y_ExtLat(jmax/4) - y_ExtLat(jmax/4-1))



    xyz_ii = int( xyz_MPLon / ( PIx2 / imax ) )

    ! Treatment below should not be necessary.
    ! xyz_MPLon is less than 2 PI, so, xyz_ii should be less than imax.
    ! But, in rare cases, xyz_ii is imax in some points.
    ! That means xyz_ii should be zero there. So, it is forced below. 
    do k = 1, kmax
      do j = 1, jmax/2
        do i = 0, imax-1
          if ( xyz_ii(i,j,k) == imax ) then
            xyz_ii(i,j,k) = 0
          end if
        end do
      end do
    end do


    do k = 1, kmax
      do j = 1, jmax/2
        do i = 0, imax-1

          ! comment out 2015/04/10
!!$          if( xyz_MPLat(i,j,k) > y_ExtLat(j) ) then
!!$            j_search_1 : do j2 = j+1, jexmax
!!$              if( y_ExtLat(j2) > xyz_MPLat(i,j,k) ) exit j_search_1
!!$            end do j_search_1
!!$            xyz_jj(i,j,k) = j2 - 1
!!$            if ( xyz_jj(i,j,k) > jexmax-2 ) xyz_jj(i,j,k) = jexmax-2
!!$          else
!!$            j_search_2 : do j2 = j-1, jexmin, -1
!!$              if( y_ExtLat(j2) < xyz_MPLat(i,j,k) ) exit j_search_2
!!$            end do j_search_2
!!$            xyz_jj(i,j,k) = j2
!!$            if ( xyz_jj(i,j,k) < jexmin+1 ) xyz_jj(i,j,k) = jexmin+1
!!$          end if
          ! trial
          if( xyz_MPLat(i,j,k) > y_ExtLat(j) ) then
            j_search_1 : do j2 = j+1, jexmax
              if( y_ExtLat(j2) > xyz_MPLat(i,j,k) ) exit j_search_1
            end do j_search_1
            xyz_jj(i,j,k) = j2 - 1
          else
            j_search_2 : do j2 = j-1, jexmin, -1
              if( y_ExtLat(j2) < xyz_MPLat(i,j,k) ) exit j_search_2
            end do j_search_2
            xyz_jj(i,j,k) = j2
          end if


!             xyz_jj(i,j,k) = int(xyz_MPLat(i,j,k)*in_deltalat) + ns + 1
!             if (y_ExtLat(xyz_jj(i,j,k)) > xyz_MPLat(i, j, k)) then 
!                xyz_jj(i,j,k) = xyz_jj(i,j,k) - 1
!             endif


        end do
      end do
    end do

    do k = 1, kmax
      do j = 1, jmax/2
        do i = 0, imax-1
          ii = xyz_ii(i,j,k)
          if ( ii-1 < iexmin ) then
            call MessageNotify( 'E', module_name, &
              & 'Longitudinal point for interporation factor calculation ' &
              & // 'is out of range of an extended array : Rank %d, (i,j,k) = (%d,%d,%d).', &
              & i = (/ myrank, i, j, k /) )
          end if
          if ( ii+2 > iexmax ) then
            call MessageNotify( 'E', module_name, &
              & 'Longitudinal point for interporation factor calculation ' &
              & // 'is out of range of an extended array : Rank %d, (i,j,k) = (%d,%d,%d).', &
              & i = (/ myrank, i, j, k /) )
          end if
          jj = xyz_jj(i,j,k)
          if ( jj-1 < jexmin ) then
            call MessageNotify( 'E', module_name, &
              & 'Latitudinal point for interporation factor calculation ' &
              & // 'is out of range of an extended array : Rank %d, (i,j,k) = (%d,%d,%d).', &
              & i = (/ myrank, i, j, k /) )
          end if
          if ( jj+2 > jexmax ) then
            call MessageNotify( 'E', module_name, &
              & 'Latitudinal point for interporation factor calculation ' &
              & // 'is out of range of an extended array : Rank %d, (i,j,k) = (%d,%d,%d).', &
              & i = (/ myrank, i, j, k /) )
          end if
        end do
      end do
    end do

    !
    ! calculation of Lagrange cubic interpolation factor 
    ! for longitudinal direction
    !
    do k = 1, kmax
      do j = 1, jmax/2
        do i = 0, imax-1
          ii = xyz_ii(i,j,k)
          jj = xyz_jj(i,j,k)
          xyza_lcifx(i,j,k,-1) =                            &
            &     ( xyz_MPLon(i,j,k) - x_ExtLon(ii  ) )     &
            &   * ( xyz_MPLon(i,j,k) - x_ExtLon(ii+1) )     &
            &   * ( xyz_MPLon(i,j,k) - x_ExtLon(ii+2) )     &
            &   * (-InvDeltaLon**3)/6.0_DP                  ! economical 
!            & / ( ( x_ExtLon(ii-1)   - x_ExtLon(ii  ) )     &
!            &   * ( x_ExtLon(ii-1)   - x_ExtLon(ii+1) )     &
!            &   * ( x_ExtLon(ii-1)   - x_ExtLon(ii+2) ) )
          xyza_lcifx(i,j,k, 0) =                            &
            &     ( xyz_MPLon(i,j,k) - x_ExtLon(ii-1) )     &
            &   * ( xyz_MPLon(i,j,k) - x_ExtLon(ii+1) )     &
            &   * ( xyz_MPLon(i,j,k) - x_ExtLon(ii+2) )     &
            &   * (0.5_DP*InvDeltaLon**3)                   ! economical 
!            & / ( ( x_ExtLon(ii  )   - x_ExtLon(ii-1) )     &
!            &   * ( x_ExtLon(ii  )   - x_ExtLon(ii+1) )     &
!            &   * ( x_ExtLon(ii  )   - x_ExtLon(ii+2) ) )
          xyza_lcifx(i,j,k, 1) =                            &
            &     ( xyz_MPLon(i,j,k) - x_ExtLon(ii-1) )     &
            &   * ( xyz_MPLon(i,j,k) - x_ExtLon(ii  ) )     &
            &   * ( xyz_MPLon(i,j,k) - x_ExtLon(ii+2) )     &
            &   * (-0.5_DP*InvDeltaLon**3)                  ! economical 
!            & / ( ( x_ExtLon(ii+1)   - x_ExtLon(ii-1) )     &
!            &   * ( x_ExtLon(ii+1)   - x_ExtLon(ii  ) )     &
!            &   * ( x_ExtLon(ii+1)   - x_ExtLon(ii+2) ) )
          xyza_lcifx(i,j,k, 2) =                            &
            &     ( xyz_MPLon(i,j,k) - x_ExtLon(ii-1) )     &
            &   * ( xyz_MPLon(i,j,k) - x_ExtLon(ii  ) )     &
            &   * ( xyz_MPLon(i,j,k) - x_ExtLon(ii+1) )     &
            &   * (InvDeltaLon**3)/6.0_DP                   ! economical 
!            & / ( ( x_ExtLon(ii+2)   - x_ExtLon(ii-1) )     &
!            &   * ( x_ExtLon(ii+2)   - x_ExtLon(ii  ) )     &
!            &   * ( x_ExtLon(ii+2)   - x_ExtLon(ii+1) ) )

          xyza_lcify(i,j,k,-1) =                            &
            &     ( xyz_MPLat(i,j,k) - y_ExtLat(jj  ) )     &
            &   * ( xyz_MPLat(i,j,k) - y_ExtLat(jj+1) )     &
            &   * ( xyz_MPLat(i,j,k) - y_ExtLat(jj+2) )     &
            & / ( ( y_ExtLat(jj-1)   - y_ExtLat(jj  ) )     &
            &   * ( y_ExtLat(jj-1)   - y_ExtLat(jj+1) )     &
            &   * ( y_ExtLat(jj-1)   - y_ExtLat(jj+2) ) )
          xyza_lcify(i,j,k, 0) =                            &
            &     ( xyz_MPLat(i,j,k) - y_ExtLat(jj-1) )     &
            &   * ( xyz_MPLat(i,j,k) - y_ExtLat(jj+1) )     &
            &   * ( xyz_MPLat(i,j,k) - y_ExtLat(jj+2) )     &
            & / ( ( y_ExtLat(jj  )   - y_ExtLat(jj-1) )     &
            &   * ( y_ExtLat(jj  )   - y_ExtLat(jj+1) )     &
            &   * ( y_ExtLat(jj  )   - y_ExtLat(jj+2) ) )
          xyza_lcify(i,j,k, 1) =                            &
            &     ( xyz_MPLat(i,j,k) - y_ExtLat(jj-1) )     &
            &   * ( xyz_MPLat(i,j,k) - y_ExtLat(jj  ) )     &
            &   * ( xyz_MPLat(i,j,k) - y_ExtLat(jj+2) )     &
            & / ( ( y_ExtLat(jj+1)   - y_ExtLat(jj-1) )     &
            &   * ( y_ExtLat(jj+1)   - y_ExtLat(jj  ) )     &
            &   * ( y_ExtLat(jj+1)   - y_ExtLat(jj+2) ) )
          xyza_lcify(i,j,k, 2) =                            &
            &     ( xyz_MPLat(i,j,k) - y_ExtLat(jj-1) )     &
            &   * ( xyz_MPLat(i,j,k) - y_ExtLat(jj  ) )     &
            &   * ( xyz_MPLat(i,j,k) - y_ExtLat(jj+1) )     &
            & / ( ( y_ExtLat(jj+2)   - y_ExtLat(jj-1) )     &
            &   * ( y_ExtLat(jj+2)   - y_ExtLat(jj  ) )     &
            &   * ( y_ExtLat(jj+2)   - y_ExtLat(jj+1) ) )
        end do
      end do
    end do


  end subroutine SLTTLagIntCubCalcFactHor

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

  subroutine SLTTLagIntCubIntHor(                         &
    & iexmin, iexmax, jexmin, jexmax,                     & ! (in)
    & xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify, xyz_ExtArr, & ! (in)
    & xyz_MPArr                                           & ! (out)
    & )
    ! 水平２次元ラグランジュ３次補間
    ! 2D Lagrangian cubic interpolation


    integer , intent(in ) :: iexmin
    integer , intent(in ) :: iexmax
    integer , intent(in ) :: jexmin
    integer , intent(in ) :: jexmax
    integer , intent(in ) :: xyz_ii(0:imax-1, 1:jmax/2, 1:kmax)
    integer , intent(in ) :: xyz_jj(0:imax-1, 1:jmax/2, 1:kmax)
    real(DP), intent(in ) :: xyza_lcifx(     0:imax-1,      1:jmax/2, 1:kmax, -1:2)
    real(DP), intent(in ) :: xyza_lcify(     0:imax-1,      1:jmax/2, 1:kmax, -1:2)
    real(DP), intent(in ) :: xyz_ExtArr(iexmin:iexmax, jexmin:jexmax, 1:kmax)
    real(DP), intent(out) :: xyz_MPArr (     0:imax-1,      1:jmax/2, 1:kmax)


    !
    ! local variables
    !
    integer :: ii
    integer :: jj
    integer :: kk

    integer :: i, j, k


    do k = 1, kmax
      do j = 1, jmax/2
        do i = 0, imax-1
          ii = xyz_ii(i,j,k)
          jj = xyz_jj(i,j,k)
          kk = k
          xyz_MPArr(i,j,k) =                                       &
            &   xyza_lcify(i,j,k,-1)                               &
            & * ( xyza_lcifx(i,j,k,-1) * xyz_ExtArr(ii-1,jj-1,kk)  &
            &   + xyza_lcifx(i,j,k, 0) * xyz_ExtArr(ii  ,jj-1,kk)  &
            &   + xyza_lcifx(i,j,k, 1) * xyz_ExtArr(ii+1,jj-1,kk)  &
            &   + xyza_lcifx(i,j,k, 2) * xyz_ExtArr(ii+2,jj-1,kk) )&
            & + xyza_lcify(i,j,k, 0)                               &
            & * ( xyza_lcifx(i,j,k,-1) * xyz_ExtArr(ii-1,jj  ,kk)  &
            &   + xyza_lcifx(i,j,k, 0) * xyz_ExtArr(ii  ,jj  ,kk)  &
            &   + xyza_lcifx(i,j,k, 1) * xyz_ExtArr(ii+1,jj  ,kk)  &
            &   + xyza_lcifx(i,j,k, 2) * xyz_ExtArr(ii+2,jj  ,kk) )&
            & + xyza_lcify(i,j,k, 1)                               &
            & * ( xyza_lcifx(i,j,k,-1) * xyz_ExtArr(ii-1,jj+1,kk)  &
            &   + xyza_lcifx(i,j,k, 0) * xyz_ExtArr(ii  ,jj+1,kk)  &
            &   + xyza_lcifx(i,j,k, 1) * xyz_ExtArr(ii+1,jj+1,kk)  &
            &   + xyza_lcifx(i,j,k, 2) * xyz_ExtArr(ii+2,jj+1,kk) )&
            & + xyza_lcify(i,j,k, 2)                               &
            & * ( xyza_lcifx(i,j,k,-1) * xyz_ExtArr(ii-1,jj+2,kk)  &
            &   + xyza_lcifx(i,j,k, 0) * xyz_ExtArr(ii  ,jj+2,kk)  &
            &   + xyza_lcifx(i,j,k, 1) * xyz_ExtArr(ii+1,jj+2,kk)  &
            &   + xyza_lcifx(i,j,k, 2) * xyz_ExtArr(ii+2,jj+2,kk) )
        end do
      end do
    end do


  end subroutine SLTTLagIntCubIntHor

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

  subroutine SLTTLagIntCubCalcFactVer(  &
    & xyz_DPSigma,                & ! (in)
    & xyza_lcifz, xyz_kk          & ! (out)
    & )
  ! 鉛直１次元ラグランジュ３次補間のための係数計算
  ! Calculation of factors for 1D Lagrangian cubic interpolation
  
    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only : z_Sigma


    real(DP), intent(in ) :: xyz_DPSigma (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyza_lcifz(0:imax-1, 1:jmax, 1:kmax, -1:2)
    integer , intent(out) :: xyz_kk    (0:imax-1, 1:jmax, 1:kmax)


    !
    ! local variables
    !

    integer :: i
    integer :: j
    integer :: k
    integer :: kk
    integer :: k2


    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1

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

          !
          ! economical method
          !
          if( xyz_DPSigma(i,j,k) > z_Sigma(k) ) then
            k_search_1 : do k2 = k, 2, -1
              if( z_Sigma(k2) > xyz_DPSigma(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( z_Sigma(k2) < xyz_DPSigma(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), 2 ), kmax-2 )
!          xyz_kk(i,j,k) = min( max( xyz_kk(i,j,k), 1 ), kmax-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)
          xyza_lcifz(i,j,k,-1) =                         &
            &     ( xyz_DPSigma(i,j,k) - z_Sigma(kk  ) ) &
            &   * ( xyz_DPSigma(i,j,k) - z_Sigma(kk+1) ) &
            &   * ( xyz_DPSigma(i,j,k) - z_Sigma(kk+2) ) &
            & / ( ( z_Sigma(kk-1)      - z_Sigma(kk  ) ) &
            &   * ( z_Sigma(kk-1)      - z_Sigma(kk+1) ) &
            &   * ( z_Sigma(kk-1)      - z_Sigma(kk+2) ) )
          xyza_lcifz(i,j,k, 0) =                         &
            &     ( xyz_DPSigma(i,j,k) - z_Sigma(kk-1) ) &
            &   * ( xyz_DPSigma(i,j,k) - z_Sigma(kk+1) ) &
            &   * ( xyz_DPSigma(i,j,k) - z_Sigma(kk+2) ) &
            & / ( ( z_Sigma(kk  )      - z_Sigma(kk-1) ) &
            &   * ( z_Sigma(kk  )      - z_Sigma(kk+1) ) &
            &   * ( z_Sigma(kk  )      - z_Sigma(kk+2) ) )
          xyza_lcifz(i,j,k, 1) =                         &
            &     ( xyz_DPSigma(i,j,k) - z_Sigma(kk-1) ) &
            &   * ( xyz_DPSigma(i,j,k) - z_Sigma(kk  ) ) &
            &   * ( xyz_DPSigma(i,j,k) - z_Sigma(kk+2) ) &
            & / ( ( z_Sigma(kk+1)      - z_Sigma(kk-1) ) &
            &   * ( z_Sigma(kk+1)      - z_Sigma(kk  ) ) &
            &   * ( z_Sigma(kk+1)      - z_Sigma(kk+2) ) )
          xyza_lcifz(i,j,k, 2) =                         &
            &     ( xyz_DPSigma(i,j,k) - z_Sigma(kk-1) ) &
            &   * ( xyz_DPSigma(i,j,k) - z_Sigma(kk  ) ) &
            &   * ( xyz_DPSigma(i,j,k) - z_Sigma(kk+1) ) &
            & / ( ( z_Sigma(kk+2)      - z_Sigma(kk-1) ) &
            &   * ( z_Sigma(kk+2)      - z_Sigma(kk  ) ) &
            &   * ( z_Sigma(kk+2)      - z_Sigma(kk+1) ) )
        end do
      end do
    end do


  end subroutine SLTTLagIntCubCalcFactVer

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

  subroutine SLTTLagIntCubIntVer(  &
    & xyz_Arr, xyza_lcifz, xyz_kk, & ! (in)
    & xyz_ArrA                     & ! (out)
    & )
  ! 鉛直１次元ラグランジュ３次補間
  ! 1D Lagrangian cubic interpolation

    real(DP), intent(in ) :: xyz_Arr   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyza_lcifz(0:imax-1, 1:jmax, 1:kmax, -1:2)
    integer , intent(in ) :: xyz_kk    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyz_ArrA  (0:imax-1, 1:jmax, 1:kmax)


    !
    ! local variables
    !

    integer :: i
    integer :: j
    integer :: k
    integer :: kk


    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          kk = xyz_kk(i,j,k)
          xyz_ArrA(i,j,k) =                              &
            &   xyza_lcifz(i,j,k,-1) * xyz_Arr(i,j,kk-1) &
            & + xyza_lcifz(i,j,k, 0) * xyz_Arr(i,j,kk  ) &
            & + xyza_lcifz(i,j,k, 1) * xyz_Arr(i,j,kk+1) &
            & + xyza_lcifz(i,j,k, 2) * xyz_Arr(i,j,kk+2)
        end do
      end do
    end do


  end subroutine SLTTLagIntCubIntVer

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

  subroutine SLTTLagIntHorMaxMin(                                &
    & iexmin, iexmax, jexmin, jexmax,                            & ! (in)
    & x_ExtLon, y_ExtLat, xyz_DPLon, xyz_DPLat, xyzf_ExtQMixB,   & ! (in)
    & xyzf_QMixMinA, xyzf_QMixMaxA                               & ! (out)
    & )
    ! 水平方向の２次元補間。
    ! 2D linear interpolation

    use mpi_wrapper, only : myrank

    integer , intent(in ) :: iexmin
    integer , intent(in ) :: iexmax
    integer , intent(in ) :: jexmin
    integer , intent(in ) :: jexmax
    real(DP), intent(in ) :: x_ExtLon(iexmin:iexmax)
                               ! 経度の拡張配列
                               ! Extended array of Lon
    real(DP), intent(in ) :: y_ExtLat(jexmin:jexmax)
                               ! 緯度の拡張配列
                               ! Extended array of Lat
    real(DP), intent(in ) :: xyz_DPLon(0:imax-1, 1:jmax/2, 1:kmax)
                               ! 上流点の経度
                               ! Lon at departure point
    real(DP), intent(in ) :: xyz_DPLat(0:imax-1, 1:jmax/2, 1:kmax)
                               ! 上流点の緯度
                               ! Lat at departure point
    real(DP), intent(in ) :: xyzf_ExtQMixB(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax)
                               ! 物質混合比の拡張配列
                               ! Extended array of mix ration of tracers
    real(DP), intent(out) :: xyzf_QMixMinA(     0:imax-1,      1:jmax/2, 1:kmax, 1:ncmax)
                               ! 上流点を囲む 4 格子点の最小値
                               ! Minimum mixing ratio of tracers at 4 grid points 
                               ! surrounding a departure point
    real(DP), intent(out) :: xyzf_QMixMaxA(     0:imax-1,      1:jmax/2, 1:kmax, 1:ncmax)
                               ! 上流点を囲む 4 格子点の最大値
                               ! Maximum mixing ratio of tracers at 4 grid points 
                               ! surrounding a departure point


    real(DP) :: DeltaLat      ! 緯度グリッド幅; grid width in meridional direction
    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:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in number of composition

    integer:: ii
    integer:: jj

    integer :: isten          ! 上流点を囲む4格子点の南西の点の座標番号（東西方向）
	                          ! Youngest index of grid points around the departure point (i-direction)
    integer :: jsten          ! 上流点を囲む4格子点の南西の点の座標番号（南北方向）
	                          ! Youngest index of grid points around the departure point (j-direction)
    real(DP) :: xy_Z(0:1,0:1)     ! 上流点を囲む 4 格子点上の混合比を格納する作業変数
                                  ! Work variable containing mixing ratio at 4 grid points around DP



    do k = 1, kmax
      do j = 1, jmax/2
        do i = 0, imax-1
          ! 上流点を囲む4点を探す
          ! Determine four grid points around the departure point
          isten = int(xyz_DPLon(i,j,k)*InvDeltaLon)
          do jj = jexmin, jexmax
            if (y_ExtLat(jj) > xyz_DPLat(i, j, k)) then
              jsten = jj-1
              exit
            endif
          enddo
          !MPIに対応できない↓
!      jsten = int(xyz_DPLat(i,j,k)*in_deltalat) + ns + 1
!      if (y_ExtLat(jsten) > xyz_DPLat(i, j, k)) then 
!        jsten = jsten - 1
!      endif

          ! Check indices
          if ( ( isten < iexmin ) .or. ( isten > iexmax ) ) then
            call MessageNotify( 'E', module_name, &
              & 'Departure point is out of range of an extended array for linear interporation ' &
              & // 'in longitudinal direction : Rank %d, i = %d.', &
              & i = (/ myrank, isten /) )
          end if
          if ( ( jsten < jexmin ) .or. ( jsten > jexmax ) ) then
            call MessageNotify( 'E', module_name, &
              & 'Departure point is out of range of an extended array for linear interporation ' &
              & // 'in latitudinal direction : Rank %d, j = %d.', &
              & i = (/ myrank, jsten /) )
          end if

          do n = 1, ncmax

            do jj = 0, 1
              do ii = 0, 1
                xy_Z(ii,jj) = xyzf_ExtQMixB(isten+ii,jsten+jj,k,n)
              end do
            end do

            xyzf_QMixMinA(i,j,k,n) = min( xy_Z(0,0), xy_Z(1,0), xy_Z(0,1), xy_Z(1,1) )
            xyzf_QMixMaxA(i,j,k,n) = max( xy_Z(0,0), xy_Z(1,0), xy_Z(0,1), xy_Z(1,1) )

          end do

        end do
      end do
    end do

  end subroutine SLTTLagIntHorMaxMin

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

  subroutine SLTTIrrLinInt(                                      &
    & iexmin, iexmax, jexmin, jexmax,                            & ! (in)
    & x_ExtLon, y_ExtLat, xyz_DPLon, xyz_DPLat, xyzf_ExtQMixB,   & ! (in)
    & xyzf_QMixA                                                 & ! (out)
    & )
    ! 水平方向の２次元補間。
    ! 2D linear interpolation

    use mpi_wrapper, only : myrank

    integer , intent(in ) :: iexmin
    integer , intent(in ) :: iexmax
    integer , intent(in ) :: jexmin
    integer , intent(in ) :: jexmax
    real(DP), intent(in ) :: x_ExtLon(iexmin:iexmax)
                               ! 経度の拡張配列
                               ! Extended array of Lon
    real(DP), intent(in ) :: y_ExtLat(jexmin:jexmax)
                               ! 緯度の拡張配列
                               ! Extended array of Lat
    real(DP), intent(in ) :: xyz_DPLon(0:imax-1, 1:jmax/2, 1:kmax)
                               ! 上流点の経度
                               ! Lon at departure point
    real(DP), intent(in ) :: xyz_DPLat(0:imax-1, 1:jmax/2, 1:kmax)
                               ! 上流点の緯度
                               ! Lat at departure point
    real(DP), intent(in ) :: xyzf_ExtQMixB(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax)
                               ! 物質混合比の拡張配列
                               ! Extended array of mix ration of tracers
    real(DP), intent(out) :: xyzf_QMixA (   0:imax-1  ,      1:jmax/2    , 1:kmax, 1:ncmax)
                               ! 次ステップの物質混合比
                               ! Mix ration of tracers at next time-step


    real(DP) :: DeltaLat      ! 緯度グリッド幅; grid width in meridional direction
    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:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in number of composition

    integer:: ii
    integer:: jj

    integer :: isten          ! 上流点を囲む4格子点の南西の点の座標番号（東西方向）
	                          ! Youngest index of grid points around the departure point (i-direction)
    integer :: jsten          ! 上流点を囲む4格子点の南西の点の座標番号（南北方向）
	                          ! Youngest index of grid points around the departure point (j-direction)
    real(DP) :: xy_Z(0:1,0:1)     ! 上流点を囲む 4 格子点上の混合比を格納する作業変数
                                  ! Work variable containing mixing ratio at 4 grid points around DP
    real(DP) :: y_Z(0:1)          ! 上流点を囲む 2 格子点上の混合比を格納する作業変数
                                  ! Work variable containing mixing ratio at 2 grid points in latitudinal direction


    do k = 1, kmax
      do j = 1, jmax/2
        do i = 0, imax-1
          ! 上流点を囲む4点を探す
          ! Determine four grid points around the departure point
          isten = int(xyz_DPLon(i,j,k)*InvDeltaLon)
          do jj = jexmin, jexmax
            if (y_ExtLat(jj) > xyz_DPLat(i, j, k)) then
              jsten = jj-1
              exit
            endif
          enddo
          !MPIに対応できない↓
!      jsten = int(xyz_DPLat(i,j,k)*in_deltalat) + ns + 1
!      if (y_ExtLat(jsten) > xyz_DPLat(i, j, k)) then 
!        jsten = jsten - 1
!      endif

          ! Check indices
          if ( ( isten < iexmin ) .or. ( isten > iexmax ) ) then
            call MessageNotify( 'E', module_name, &
              & 'Departure point is out of range of an extended array for linear interporation ' &
              & // 'in longitudinal direction : Rank %d, i = %d.', &
              & i = (/ myrank, isten /) )
          end if
          if ( ( jsten < jexmin ) .or. ( jsten > jexmax ) ) then
            call MessageNotify( 'E', module_name, &
              & 'Departure point is out of range of an extended array for linear interporation ' &
              & // 'in latitudinal direction : Rank %d, j = %d.', &
              & i = (/ myrank, jsten /) )
          end if

          do n = 1, ncmax

            do jj = 0, 1
              do ii = 0, 1
                xy_Z(ii,jj) = xyzf_ExtQMixB(isten+ii,jsten+jj,k,n)
              end do
            end do
            DeltaLat = y_ExtLat(jsten+1) - y_ExtLat(jsten)

            do jj = 0, 1
              y_Z(jj) =   ( xy_Z(1,jj) - xy_Z(0,jj) ) / DeltaLon   &
                &         * ( xyz_DPLon(i,j,k) - x_ExtLon(isten) ) + xy_Z(0,jj)
            end do
            xyzf_QMixA(i,j,k,n) =   ( y_Z(1) - y_Z(0) ) / DeltaLat           &
              &                     * ( xyz_DPLat(i,j,k) - y_ExtLat(jsten) ) &
              &                   + y_Z(0)

          end do

        end do
      end do
    end do

  end subroutine SLTTIrrLinInt

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

  subroutine SLTTIrrHerIntK13(              &
   & iexmin, iexmax, jexmin, jexmax,                            & ! (in)
   & x_ExtLon, y_ExtLat, xyz_DPLon, xyz_DPLat, xyz_ExtQMixB,    & ! (in)
   & xyz_ExtQMixB_dlon, xyz_ExtQMixB_dlat, xyz_ExtQMixB_dlonlat,& ! (in) 
!    & xyz_ExtQMixB_dlon2, xyz_ExtQMixB_dlat2,        & ! (in) fxx, fyy
!    & xyz_ExtQMixB_dlon2lat, xyz_ExtQMixB_dlonlat2,  & ! (in) fxxy, fxyy
!    & xyz_ExtQMixB_dlon2lat2,                                  & ! (in) fxxyy       	
   & SLTTIntHor,                                                        & ! (in) 
   & xyz_QMixA        & ! (out)
   & )
    ! 水平方向の２次元補間。Enomoto (2008, SOLA)で提案された「スペクトルで計算した微分値を用いた双３次補間」を
    ! 発展させた方法、スペクトル微分を用いた変則エルミート５次補間を方向分離して行う。５次補間のたびに Sun et al. (1996, MWR)
    ! の単調フィルタを修正したものを適用する。
    ! 2D Interpolation. Spectral transformation is used for calculation of derivatives, which are used
    ! for Irregular Hermite quintic interpolation. The original idea of using Spectral transformation for derivatives
    ! is presented by Enomoto (2008, SOLA).
    ! Monotonicity filter presented by Sun et al. (1996, MWR) is partly modified and used after each interpolation.

    use sltt_const , only : PIx2
    use mpi_wrapper, only : myrank
    use gridset, only: lmax    ! スペクトルデータの配列サイズ
                               ! Size of array for spectral data

    integer , intent(in ) :: iexmin
    integer , intent(in ) :: iexmax
    integer , intent(in ) :: jexmin
    integer , intent(in ) :: jexmax
    real(DP), intent(in ) :: x_ExtLon(iexmin:iexmax)
                               ! 経度の拡張配列
                               ! Extended array of Lon
    real(DP), intent(in ) :: y_ExtLat(jexmin:jexmax)
                               ! 緯度の拡張配列
                               ! Extended array of Lat
    real(DP), intent(in ) :: xyz_DPLon(0:imax-1, 1:jmax/2, 1:kmax)
                               ! 上流点の経度
                               ! Lon at departure point
    real(DP), intent(in ) :: xyz_DPLat(0:imax-1, 1:jmax/2, 1:kmax)
                               ! 上流点の緯度
                               ! Lat at departure point
    real(DP), intent(in ) :: xyz_ExtQMixB(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax)
                               ! 物質混合比の拡張配列
                               ! Extended array of mix ration of tracers
    real(DP), intent(in) :: xyz_ExtQMixB_dlon(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax)
                               ! 物質混合比の経度微分の拡張配列
                               ! Extended array of zonal derivative of the mix ration
    real(DP), intent(in) :: xyz_ExtQMixB_dlat(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax)
                               ! 物質混合比の緯度微分の拡張配列
                               ! Extended array of meridional derivative of the mix ration
    real(DP), intent(in) :: xyz_ExtQMixB_dlonlat(iexmin:iexmax, jexmin:jexmax, 1:kmax, 1:ncmax)            
                               ! 物質混合比の緯度経度微分の拡張配列
                               ! Extended array of zonal and meridional derivative of the mix ration
    character(TOKEN), intent(in):: SLTTIntHor
                               ! 水平方向の補間方法を指定するキーワード
                               ! Keyword for Interpolation Method for Horizontal direction
    real(DP), intent(out) :: xyz_QMixA (   0:imax-1  ,      1:jmax/2    , 1:kmax, 1:ncmax)
                               ! 次ステップの物質混合比
                               ! Mix ration of tracers at next time-step

!---fxx, fyy, fxxy, fxyy, fxxyy
!    real(DP), intent(inout) :: xyz_ExtQMixB_dlon2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax)
!    real(DP), intent(inout) :: xyz_ExtQMixB_dlat2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax)
!    real(DP), intent(inout) :: xyz_ExtQMixB_dlon2lat(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax)            
!    real(DP), intent(inout) :: xyz_ExtQMixB_dlonlat2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax)                
!    real(DP), intent(inout) :: xyz_ExtQMixB_dlon2lat2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax)                



    real(DP) :: DeltaLat      ! 緯度グリッド幅; grid width in meridional direction
    integer:: i, ii           ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j, jj           ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in number of composition

    integer :: isten          ! 上流点を囲む4格子点の南西の点の座標番号（東西方向）
	                          ! Youngest index of grid points around the departure point (i-direction)
    integer :: jsten          ! 上流点を囲む4格子点の南西の点の座標番号（南北方向）
	                          ! Youngest index of grid points around the departure point (j-direction)
    integer  :: num           ! 配列配置のための作業変数
	                          ! Work variable for array packing
    real(DP) :: a_z(16)       ! 上流点を囲む16格子点上の混合比を格納する作業変数
                              ! Work variable containing mixing ratio at 16 grid points around DP
    real(DP) :: a_zx(16)      ! 上流点を囲む16格子点上の混合比の経度微分を格納する作業変数
                              ! Work variable containing zonal derivative of mixing ratio at 16 grid points around DP
    real(DP) :: a_zy(16)      ! 上流点を囲む16格子点上の混合比の緯度微分を格納する作業変数
                              ! Work variable containing meridional derivative of mixing ratio at 16 grid points around DP
    real(DP) :: a_zxy(16)     ! 上流点を囲む16格子点上の混合比の緯度経度微分を格納する作業変数
                              ! Work variable containing zonal and meridional derivative of mixing ratio at 16 grid points around DP
!    real(DP):: a_zxx(16), a_zyy(16), a_zxxy(16), a_zxyy(16), a_zxxyy(16)


    ! Check whether a longitude of departure point is within an extended array
    call SLTTLagIntChkDPLon(       &
      & SLTTIntHor,                      & ! (in)
      & iexmin, iexmax, jexmin, jexmax,  & ! (in)
      & x_ExtLon, xyz_DPLon              & ! (in)
      & )
    ! Check whether a latitude of departure point is within an extended array
    call SLTTLagIntChkDPLat(             &
      & SLTTIntHor,                      & ! (in)
      & iexmin, iexmax, jexmin, jexmax,  & ! (in)
      & y_ExtLat, xyz_DPLat              & ! (in)
      & )


 do k = 1, kmax
   do j = 1, jmax/2
     do i = 0, imax-1
       ! 上流点を囲む4点を探す
       ! Determine four grid points around the departure point
       isten = int(xyz_DPLon(i,j,k)*InvDeltaLon)
       do jj = jexmin, jexmax
         if (y_ExtLat(jj) > xyz_DPLat(i, j, k)) then
           jsten = jj-1
           exit
         endif
       enddo
      !MPIに対応できない↓
!      jsten = int(xyz_DPLat(i,j,k)*in_deltalat) + ns + 1
!      if (y_ExtLat(jsten) > xyz_DPLat(i, j, k)) then 
!        jsten = jsten - 1
!      endif


        ! 水平方向の補間方法の選択
        select case (SLTTIntHor)

        case ("HQ") ! 方向分離型2次元変則エルミート5次補間;  2D Hermite quintic interpolation
          do n = 1, ncmax
            ! ２次元補間のための配列配置（ワイド）
            ! Array packing for 2D interpolation
            do jj = -1, 2
              num = (jj+1)*4 + 2
              do ii = -1, 2
                a_z(ii+num) = xyz_ExtQMixB(isten+ii, jsten+jj, k, n)
                a_zx(ii+num) = xyz_ExtQMixB_dlon(isten+ii, jsten+jj, k, n)
                a_zy(ii+num) = xyz_ExtQMixB_dlat(isten+ii, jsten+jj, k, n)
                a_zxy(ii+num) = xyz_ExtQMixB_dlonlat(isten+ii, jsten+jj, k, n)
              enddo
            enddo

            ! 方向分離型2次元変則エルミート5次補間 
            ! 2D Hermite quintic interpolation
            xyz_QMixA(i,j,k,n) = SLTTIrrHerIntQui2DHor(a_z, a_zx, a_zy, a_zxy, &
              & y_ExtLat(jsten-1)-y_ExtLat(jsten), y_ExtLat(jsten+1)-y_ExtLat(jsten), y_ExtLat(jsten+2)-y_ExtLat(jsten),&
              &   xyz_DPLon(i, j, k)-x_ExtLon(isten),  xyz_DPLat(i, j, k)-y_ExtLat(jsten) )
          enddo


        case("HC") !方向分離型2次元エルミート３次補間; 2D Hermite Cubic Interpolation
          do n = 1, ncmax
          ! ２次元補間のための配列配置
          !    4           3
          !         x
          !    1           2           
            a_z(1) = xyz_ExtQMixB(isten, jsten, k, n)
            a_zx(1) = xyz_ExtQMixB_dlon(isten, jsten, k, n)
            a_zy(1) = xyz_ExtQMixB_dlat(isten, jsten, k, n)
            a_zxy(1) = xyz_ExtQMixB_dlonlat(isten, jsten, k, n)
!               a_zxx(1) = xyz_ExtQMixB_dlon2(isten, jsten, k, n)
!               a_zyy(1) = xyz_ExtQMixB_dlat2(isten, jsten, k, n)
!               a_zxxy(1) = xyz_ExtQMixB_dlon2lat(isten, jsten, k, n)
!               a_zxyy(1) = xyz_ExtQMixB_dlonlat2(isten, jsten, k, n)
!               a_zxxyy(1) = xyz_ExtQMixB_dlon2lat2(isten, jsten, k, n)

            a_z(2) = xyz_ExtQMixB(isten+1, jsten, k, n)
            a_zx(2) = xyz_ExtQMixB_dlon(isten+1, jsten, k, n)
            a_zy(2) = xyz_ExtQMixB_dlat(isten+1, jsten, k, n)
            a_zxy(2) = xyz_ExtQMixB_dlonlat(isten+1, jsten, k, n)
!               a_zxx(2) = xyz_ExtQMixB_dlon2(isten+1, jsten, k, n)
!               a_zyy(2) = xyz_ExtQMixB_dlat2(isten+1, jsten, k, n)
!               a_zxxy(2) = xyz_ExtQMixB_dlon2lat(isten+1, jsten, k, n)
!               a_zxyy(2) = xyz_ExtQMixB_dlonlat2(isten+1, jsten, k, n)
!               a_zxxyy(2) = xyz_ExtQMixB_dlon2lat2(isten+1, jsten, k, n)

            a_z(3) = xyz_ExtQMixB(isten+1, jsten+1, k, n)
            a_zx(3) = xyz_ExtQMixB_dlon(isten+1, jsten+1, k, n)
            a_zy(3) = xyz_ExtQMixB_dlat(isten+1, jsten+1, k, n)
            a_zxy(3) = xyz_ExtQMixB_dlonlat(isten+1, jsten+1, k, n)
!               a_zxx(3) = xyz_ExtQMixB_dlon2(isten+1, jsten+1, k, n)
!               a_zyy(3) = xyz_ExtQMixB_dlat2(isten+1, jsten+1, k, n)
!               a_zxxy(3) = xyz_ExtQMixB_dlon2lat(isten+1, jsten+1, k, n)
!               a_zxyy(3) = xyz_ExtQMixB_dlonlat2(isten+1, jsten+1, k, n)
!               a_zxxyy(3) = xyz_ExtQMixB_dlon2lat2(isten+1, jsten+1, k, n)

            a_z(4) = xyz_ExtQMixB(isten, jsten+1, k, n)
            a_zx(4) = xyz_ExtQMixB_dlon(isten, jsten+1, k, n)
            a_zy(4) = xyz_ExtQMixB_dlat(isten, jsten+1, k, n)
            a_zxy(4) = xyz_ExtQMixB_dlonlat(isten, jsten+1, k, n)
!               a_zxx(4) = xyz_ExtQMixB_dlon2(isten, jsten+1, k, n)
!               a_zyy(4) = xyz_ExtQMixB_dlat2(isten, jsten+1, k, n)
!               a_zxxy(4) = xyz_ExtQMixB_dlon2lat(isten, jsten+1, k, n)
!               a_zxyy(4) = xyz_ExtQMixB_dlonlat2(isten, jsten+1, k, n)
!               a_zxxyy(4) = xyz_ExtQMixB_dlon2lat2(isten, jsten+1, k, n)

            DeltaLat = y_ExtLat(jsten+1) - y_ExtLat(jsten)

            !方向分離型2次元エルミート３次補間
            xyz_QMixA(i,j,k,n) = SLTTHerIntCub2D(a_z(1:4), a_zx(1:4), a_zy(1:5), a_zxy(1:4), DeltaLon, DeltaLat,&
                                &  xyz_DPLon(i, j, k)-x_ExtLon(isten), xyz_DPLat(i, j, k)-y_ExtLat(jsten) )

!方向分離型2次元エルミート5次補間
!    	xyz_QMixA(i,j,k) = hqint2D(a_z, a_zx, a_zy, a_zxy, a_zxx, a_zyy, a_zxxy, a_zxyy, a_zxxyy, deltalon, deltalat, &
!		&   xyz_DPLon(i, j, k)-x_ExtLon(isten),  xyz_DPLat(i, j, k)-y_ExtLat(jsten) )
         enddo

        case DEFAULT
          call MessageNotify( 'E', module_name, &
            & 'GIVE CORRECT KEYWORD FOR <SLTTIntHor> IN NAMELIST.' )
      end select



    enddo
  enddo
enddo

end subroutine SLTTIrrHerIntK13


  function SLTTHerIntCub2D(f, fx, fy, fxy, dx, dy, Xix, Xiy) result (fout)
    !２次元エルミート３次補間  
    ! 2D Hermite Cubic Interpolation  
    !    4----b------3
    !         |
    !         X
    !         |
    !    1----a------2
    ! Xix:点1と点aとの間隔、Xiy:点aと点Xとの間隔
    ! Xix:distance between 1 and a, Xiy:distance between a and X

    implicit none
    real(DP), dimension(4), intent(in) :: f, fx, fy, fxy
    real(DP), intent(in) :: dx, dy, Xix, Xiy
    real(DP) :: fout

    !------Internal variables-------
    real(DP) :: fa, fb, fya, fyb


    ! 点1と点2から点aでの値を求める
    ! interpolate a from 1 and 2 
    fa = SLTTHerIntCub1D(f(1), f(2), fx(1), fx(2), dx, Xix)
    fya = SLTTHerIntCub1D(fy(1), fy(2), fxy(1), fxy(2), dx, Xix)

    ! 点4と点3から点bでの値を求める
    ! interpolate b from 4 and 3 
    fb = SLTTHerIntCub1D(f(4), f(3), fx(4), fx(3), dx, Xix)
    fyb = SLTTHerIntCub1D(fy(4), fy(3), fxy(4), fxy(3), dx, Xix)

    ! 点aと点bから点Xをでの値を求める
    ! interpolate X from a and b
    fout = SLTTHerIntCub1D(fa, fb, fya, fyb, dy, Xiy)

  end function SLTTHerIntCub2D


  function SLTTHerIntCub1D(f1, f2, g1, g2, dx, Xi) result (fout)
    !エルミート３次補間
    ! 1D Hermite Cubic Interpolation
    !    1-----x------2
    ! f:関数値、g:微分値、dx:点１と点２の間隔、Xi:点１と補間する点xとの間隔
    ! f:function value, g:derivative, dx:distance between 1 and 2, Xi:distance between 1 and x
    implicit none
    real(DP), intent(in) :: f1, f2, g1, g2
    real(DP), intent(in) :: dx, Xi
    real(DP) :: fout
    !------Internal variables-------
    real(DP):: a, b
    real(DP):: indx

    indx = 1.0_DP/dx

    a = (g1 + g2)*indx*indx + 2.0_DP*(f1 - f2)*indx*indx*indx
    b = 3.0_DP*(f2 - f1)*indx*indx - (2.0_DP*g1 + g2)*indx
    fout = a*Xi*Xi*Xi + b*Xi*Xi + g1*Xi + f1

  end function SLTTHerIntCub1D



!function centdif4(f,x) result (g3out)
!! 不等間隔格子での4次精度の中心差分。
!! 関数値 f1, f2, f3, f4, f5 から 点3の x微分 g3 をもとめる。
!real(DP), intent(in) :: f(5), x(5)
!real(DP) :: g3out
!
!!---内部変数---
!real(DP) :: s1, s2, t1, t2, gtmp1, gtmp2
!
!s1 = x(3) - x(2)
!t1 = x(4) - x(3)
!s2 = x(3) - x(1)
!t2 = x(5) - x(3)
!!準備｜不等間隔格子での2次精度中心差分(http://ruby.gfd-dennou.org/products/numru-derivative/derivative/doc/document.pdf)
!gtmp1 = (s1*s1*f(4) + (t1*t1 - s1*s1)*f(3) - t1*t1*f(2))/(s1*t1*(s1 + t1))
!gtmp2 = (s2*s2*f(5) + (t2*t2 - s2*s2)*f(3) - t2*t2*f(1))/(s2*t2*(s2 + t2))
!
!if (gtmp1 == 0.0) then !2次精度中心差分で傾きゼロの点は、そのまま出力する
!g3out = 0.0_DP
!else
!g3out = (gtmp1*s2*t2 - gtmp2*s1*t1)/(s2*t2 - s1*t1) 
!endif
!
!end function centdif4


!  function hqint2D(f, fx, fy, fxy, fxx, fyy, fxxy, fxyy, fxxyy, dx, dy, Xix, Xiy) result (fout)
!!２次元エルミート５次補間（２階微分使用）
!! 2D Hermite Quintic Interpolation (using 2nd derivatives)
!!    4----b------3
!!         |
!!         X
!!         |
!!    1----a------2
!! Xix:点1と点aとの間隔、Xiy:点aと点Xとの間隔
!! Xix:distance between 1 and a, Xiy:distance between a and X
!
!    implicit none
!    real(DP), dimension(4), intent(in) :: f, fx, fy, fxy, fxx, fyy, fxxy, fxyy, fxxyy
!    real(DP), intent(in) :: dx, dy, Xix, Xiy
!    real(DP) :: fout
!!------interlan variables-------
!    real(DP) :: fa, fb, fya, fyb, fyya, fyyb
!
!
!! 点1と点2から点aでの値を求める
!! interpolate a from 1 and 2 
!fa = SLTTHerIntQui1D(f(1), f(2), fx(1), fx(2), fxx(1), fxx(2), dx, Xix)
!fya = SLTTHerIntQui1D(fy(1), fy(2), fxy(1), fxy(2), fxxy(1), fxxy(2), dx, Xix)
!fyya = SLTTHerIntQui1D(fyy(1), fyy(2), fxyy(1), fxyy(2), fxxyy(1), fxxyy(2), dx, Xix)
!
!
!! 点4と点3から点bでの値を求める
!! interpolate b from 4 and 3 
!fb = SLTTHerIntQui1D(f(4), f(3), fx(4), fx(3), fxx(4), fxx(3), dx, Xix)
!fyb = SLTTHerIntQui1D(fy(4), fy(3), fxy(4), fxy(3), fxxy(4), fxxy(3), dx, Xix)
!fyyb = SLTTHerIntQui1D(fyy(4), fyy(3), fxyy(4), fxyy(3), fxxyy(4), fxxyy(3), dx, Xix)
!
!! 点aと点bから点Xをでの値を求める
!! interpolate X from a and b 
!fout = SLTTHerIntQui1D(fa, fb, fya, fyb, fyya, fyyb, dy, Xiy)
!
!    end function hqint2D



  function SLTTIrrHerIntQui2DHor(f, fx, fy, fxy, dy21, dy23, dy24, Xix, Xiy) result (fout)
    ! ２次元変則エルミート５次補間（２階微分不使用）
    ! 2D Hermite Quintic Interpolation (without 2nd derivatives)
    !
    ! 13-14-d-15--16
    ! |   | | |   |
    ! 9--10-c-11--12
    ! |   | x |   |
    ! 5---6-b-7---8
    ! |   | | |   |
    ! 1---2-a-3---4
    !
    ! Xix:点6と点aとの間隔、Xiy:点aと点Xとの間隔
    ! dy21=lat(5)-lat(1), dy23=lat(5)-lat(9), dy24=lat(5)-lat(13) 
    ! Xix:distance between 6 and a, Xiy:distance between b and x
    ! dy21=lat(5)-lat(1), dy23=lat(5)-lat(9), dy24=lat(5)-lat(13) 

    implicit none
    real(DP), dimension(16), intent(in) :: f, fx, fy, fxy
    real(DP), intent(in) :: dy21, dy23, dy24, Xix, Xiy
    real(DP) :: fout
    !------internal variables-------
    real(DP) :: fa, fb, fc, fd, fyb, fyc, fya, fyd, lmin, lmax

    ! 点1-4から点aでの値を求める
    ! interpolate a from points 1-4
    fa = SLTTIrrHerIntQui1DUniLon(f(1), f(2), f(3), f(4), fx(2), fx(3),  Xix)

    ! 点5-8から点bでの値を求める
    ! interpolate b from points 5-8
    fb = SLTTIrrHerIntQui1DUniLon(f(5), f(6), f(7), f(8), fx(6), fx(7),  Xix)
    fyb = SLTTIrrHerIntQui1DUniLon(fy(5), fy(6), fy(7), fy(8), fxy(6), fxy(7),  Xix)

    ! 点9-12から点cでの値を求める
    ! interpolate c from points 9-12
    fc = SLTTIrrHerIntQui1DUniLon(f(9), f(10), f(11), f(12), fx(10), fx(11),  Xix)
    fyc = SLTTIrrHerIntQui1DUniLon(fy(9), fy(10), fy(11), fy(12), fxy(10), fxy(11),  Xix)

    ! 点13-16から点dでの値を求める
    ! interpolate d from points 13-16
    fd = SLTTIrrHerIntQui1DUniLon(f(13), f(14), f(15), f(16), fx(14), fx(15),  Xix)

    ! 点a-dから点Xをでの値を求める
    ! interpolate X from points a-d
    fout = SLTTIrrHerIntQui1DNonUni(fa, fb, fc, fd, fyb, fyc, dy21, dy23, dy24, Xiy)
    !等間隔格子に近似しても、精度はほとんど落ちない。計算は軽くなる。
    !fout = SLTTIrrHerIntQui1DUni(fa, fb, fc, fd, fyb, fyc, dy23, Xiy)

  end function SLTTIrrHerIntQui2DHor



!  function SLTTHerIntQui1D(f1, f2, g1, g2, h1, h2, dx, Xi) result (fout)
!! エルミート５次補間（２階微分使用）
!! 1D Hermite Quintic Interpolation (using 2nd derivatives)
!!    1-----x------2
!! f:関数値、g:微分値、h:2階微分、dx:点１と点２の間隔、Xi:点１と補間する点xとの間隔
!! f(x) = a_5*x^5 + a_4*x^4 + a_3*x^3 + a_2*x^2 + a_1*x + a_0 で補間する
!! f:function value, g:derivative, h:2nd derivative, 
!! dx:distance between 1 and 2, Xi:distance between 1 and x
!! f(x) = a_5*x^5 + a_4*x^4 + a_3*x^3 + a_2*x^2 + a_1*x + a_0
!
!    implicit none
!    real(DP), intent(in) :: f1, f2, g1, g2, h1, h2
!    real(DP), intent(in) :: dx, Xi
!    real(DP) :: fout
!!------内部変数-------
!    real(DP):: a(0:5)
!    real(DP):: indx
!
!indx = 1.0_DP/dx
!
!a(0) = f1
!a(1) = g1
!a(2) = 0.5_DP*h1
!a(3) = 10.0_DP*(-f1+f2)*indx*indx*indx + (-6.0_DP*g1-4.0_DP*g2)*indx*indx + (-1.5_DP*h1+0.5_DP*h2)*indx
!a(4) = 15.0_DP*(f1-f2)*indx*indx*indx*indx + (8.0_DP*g1+7.0_DP*g2)*indx*indx*indx + (1.5_DP*h1-h2)*indx*indx
!a(5) = 6.0_DP*(-f1+f2)*indx*indx*indx*indx*indx + 3.0_DP*(-g1-g2)*indx*indx*indx*indx + 0.5_DP*(-h1+h2)*indx*indx*indx
!
!fout = a(5)*Xi*Xi*Xi*Xi*Xi + a(4)*Xi*Xi*Xi*Xi + a(3)*Xi*Xi*Xi + a(2)*Xi*Xi + a(1)*Xi + a(0)
!
!    end function SLTTHerIntQui1D



  function SLTTIrrHerIntQui1DUni(f1, f2, f3, f4, g2, g3, dx, Xi) result (fout)
    ! 変則エルミート５次補間（２階微分不使用）
    ! 1D Hermite Quintic Interpolation (without 2nd derivatives)
    ! 等間隔格子の場合
    ! equal separation
    !    1---2--x-3---4
    ! f:関数値、g:微分値、、dx:点の間隔、Xi:点２と補間する点xとの間隔
    ! f:function value, g:derivative,
    ! dx:equal separation of each point, Xi:distance between 2 and x
    ! f(x) = a_5*x^5 + a_4*x^4 + a_3*x^3 + a_2*x^2 + a_1*x + a_0 
    ! f1 = f(-X), f2 = f(0), f3 = f(X), f4 = f(2X)

    implicit none
    real(DP), intent(in) :: f1, f2, f3, f4, g2, g3
    real(DP), intent(in) :: dx, Xi
    real(DP) :: fout
    !------internal variables-------
    real(DP):: a(0:5)
    real(DP):: indx

    indx = 1.0_DP/dx

    a(0) = f2
    a(1) = g2
    a(5) = (0.75_DP*f3 - (f1 - f4)/12.0_DP -0.5_DP*( g3 + a(1))*dx -0.75_DP*a(0))*indx**5
    a(3) = -a(5)*dx*dx - a(1)*indx*indx + (f3 - f1)*0.5_DP*indx**3
    a(4) = ( (g3 + a(1))*0.5_DP*dx - f3 +a(0) )*indx**4 - 1.5_DP*a(5)*dx - 0.5_DP*a(3)*indx
    a(2) = ( (f3 + f1)*0.5_DP -a(0))*indx*indx -a(4)*dx*dx 

    !fout = a(5)*Xi*Xi*Xi*Xi*Xi + a(4)*Xi*Xi*Xi*Xi &
    !&    + a(3)*Xi*Xi*Xi + a(2)*Xi*Xi + a(1)*Xi + a(0)
    fout =   (((((a(5)*Xi + a(4))*Xi+ a(3))*Xi)+ a(2))*Xi + a(1))*Xi + a(0)


    ! Monotonic filter
#ifdef SLTTFULLMONOTONIC
    fout = min(max(f2,f3), max(min(f2,f3), fout))
#elif SLTT2D1DMONOTONIC
    ! Do nothing
#else
    if ((f2 - f1)*(f4 - f3)>=0.0_DP) then
      fout = min(max(f2,f3), max(min(f2,f3), fout))
    endif
#endif

  end function SLTTIrrHerIntQui1DUni

  function SLTTIrrHerIntQui1DUniLon(f1, f2, f3, f4, g2, g3, Xi) result (fout)
    ! 変則エルミート５次補間（２階微分不使用）
    ! 1D Hermite Quintic Interpolation (without 2nd derivatives)
    ! 経度方向（等間隔）専用
    ! equal separation
    !    1---2--x-3---4
    ! f:関数値、g:微分値、、dx:点の間隔、Xi:点２と補間する点xとの間隔
    ! f:function value, g:derivative,
    ! dx:equal separation of each point, Xi:distance between 2 and x
    ! f(x) = a_5*x^5 + a_4*x^4 + a_3*x^3 + a_2*x^2 + a_1*x + a_0 
    ! f1 = f(-X), f2 = f(0), f3 = f(X), f4 = f(2X)

    implicit none
    real(DP), intent(in) :: f1, f2, f3, f4, g2, g3
    real(DP), intent(in) :: Xi
    real(DP) :: fout
    !------internal variables-------
    real(DP):: a(0:5)


    a(0) = f2
    a(1) = g2
    a(5) = (0.75_DP*f3 - (f1 - f4)/12.0_DP -0.5_DP*( g3 + a(1))*DeltaLon -0.75_DP*a(0))*InvDeltaLon**5
    a(3) = -a(5)*DeltaLon*DeltaLon - a(1)*InvDeltaLon*InvDeltaLon + (f3 - f1)*0.5_DP*InvDeltaLon**3
    a(4) = ( (g3 + a(1))*0.5_DP*DeltaLon - f3 +a(0) )*InvDeltaLon**4 - 1.5_DP*a(5)*DeltaLon - 0.5_DP*a(3)*InvDeltaLon
    a(2) = ( (f3 + f1)*0.5_DP -a(0))*InvDeltaLon*InvDeltaLon -a(4)*DeltaLon*DeltaLon 

    !fout = a(5)*Xi*Xi*Xi*Xi*Xi + a(4)*Xi*Xi*Xi*Xi &
    !&    + a(3)*Xi*Xi*Xi + a(2)*Xi*Xi + a(1)*Xi + a(0)
    fout =   (((((a(5)*Xi + a(4))*Xi+ a(3))*Xi)+ a(2))*Xi + a(1))*Xi + a(0)


    ! Monotonic filter
#ifdef SLTTFULLMONOTONIC
    fout = min(max(f2,f3), max(min(f2,f3), fout))
#elif SLTT2D1DMONOTONIC
    ! Do nothing
#else
    if ((f2 - f1)*(f4 - f3)>=0.0_DP) then
      fout = min(max(f2,f3), max(min(f2,f3), fout))
    endif
#endif

  end function SLTTIrrHerIntQui1DUniLon



  function SLTTIrrHerIntQui1DNonUni(f1, f2, f3, f4, g2, g3, dx21, dx23, dx24, Xi) result (fout)
    ! 変則エルミート５次補間（２階微分不使用）
    ! 1D Hermite Quintic Interpolation
    ! 不等間隔格子の場合
    ! non-equal separation
    !    1---2--x-3---4
    ! f:関数値、g:微分値、、dx:点の間隔、Xi:点２と補間する点xとの間隔
    ! f:function value, g:derivative
    ! dx21: lon(1)-lon(2), dx23: lon(3)-lon(2), dx24: lon(4)-lon(2), 
    ! f(x) = a_5*x^5 + a_4*x^4 + a_3*x^3 + a_2*x^2 + a_1*x + a_0
    ! f1 = f(dx21), f2 = f(0), f3 = f(dx23), f4 = f(dx24)

    implicit none
    real(8), intent(in) :: f1, f2, f3, f4, g2, g3
    real(8), intent(in) :: dx21, dx23, dx24, Xi
    real(8) :: fout, r, t
    !------Internal variables-------
    real(8):: a(0:5)
    real(8):: indx
    real(8):: Y1, Y3, Y4, Z3, Xi2
    integer :: n

    ! 計算効率化のため、dx21, dx23, dx24 を dx23 で正規化する。
    ! このとき、1階微分値は × dx23 する必要がある。
    r = dx21/dx23
    t = dx24/dx23

    a(0) = f2
    a(1) = g2*dx23

    Y1 = f1 - a(0) -a(1)*r
    Y3 = f3 - a(0) -a(1)
    Y4 = f4 - a(0) -a(1)*t
    Z3 = g3*dx23 - a(1)

    ! 連立方程式
    ! a(5) + a(4) + a(3) + a(2) = Y3
    ! a(5)r^5 + a(4)r^4 + a(3)r^3 + a(2)r^2 = Y1
    ! a(5)t^5 + a(4)t^4 + a(3)t^3 + a(2)t^2 = Y4
    ! 5a(5) + 4a(4) + 3a(3) + 2a(2) = Z3
    ! の解は

    a(5) = Y1/( (r-1.0_DP)*(r-1.0_DP)*r*r*(r-t) ) - Y4 / ( (t-1.0_DP)*(t-1.0_DP)*t*t*(r-t) ) &
    &     - (4.0_DP + 2.0_DP*r*t -3.0_DP*(r+t))*Y3 / ( (r-1.0_DP)*(r-1.0_DP)*(t-1.0_DP)*(t-1.0_DP) )&
    &     + Z3 / ( (r-1.0_DP)*(t-1.0_DP) )

    a(4) = -(t+2.0_DP)*Y1 / ((r-1.0_DP)*(r-1.0_DP)*r*r*(r-t) ) &
    &      +(r+2.0_DP)*Y4 / ((t-1.0_DP)*(t-1.0_DP)*t*t*(r-t) ) &
    &      +(5.0_DP - 3.0_DP*(r*r + r*t + t*t) +2.0_DP*r*t*(r+t))*Y3 / ( (r-1.0_DP)*(r-1.0_DP)*(t-1.0_DP)*(t-1.0_DP) ) &
    &      -(r+t+1.0_DP)*Z3/((r-1.0_DP)*(t-1.0_DP))

    a(3) = -2.0_DP*Y3 + Z3 -3.0_DP*a(5) - 2.0_DP*a(4)

    a(2) = Y3 - a(5) - a(4) - a(3)

    Xi2 = Xi/dx23

    !fout = a(5)*Xi2*Xi2*Xi2*Xi2*Xi2 + a(4)*Xi2*Xi2*Xi2*Xi2 &
    !&    + a(3)*Xi2*Xi2*Xi2 + a(2)*Xi2*Xi2 + a(1)*Xi2 + a(0)
    fout =   (((((a(5)*Xi2 + a(4))*Xi2+ a(3))*Xi2)+ a(2))*Xi2 + a(1))*Xi2 + a(0)

    ! Monotonic filter 
#ifdef SLTTFULLMONOTONIC
    fout = min(max(f2,f3), max(min(f2,f3), fout))
#elif SLTT2D1DMONOTONIC
    ! Do nothing
#else
    if ((f2 - f1)*(f4 - f3)>=0.0_DP) then
      fout = min(max(f2,f3), max(min(f2,f3), fout))
    endif
#endif

    end function SLTTIrrHerIntQui1DNonUni


    
!	function judgeSun1996(f1, f2, f3, f4) result(tf)
!	! Sun et al. (1996, MWR)の単調フィルタ条件の判定
!	! Judge the condition for Sun et al. (1996, MWR) monotonic filter
!	
!	real(8), intent(in) :: f1, f2, f3, f4
!	logical :: tf
!	
!	real(8) :: ineq1
!	
!	ineq1 = (f2 - f1)*(f4 - f3)
!	
!	if((ineq1>=0.0_8)) then
!	tf = .true.
!	else
!	tf = .false.
!	endif
!	
!	end function judgeSun1996


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

  !
  ! Checking latitude of departure point
  !
  subroutine SLTTLagIntChkDPLon(       &
    & SLTTIntHor,                      & ! (in)
    & iexmin, iexmax, jexmin, jexmax,  & ! (in)
    & x_ExtLon, xyz_DPLon              & ! (in)
    & )

    !
    ! MPI
    !
    use mpi_wrapper, only : myrank
                              ! Number of MPI rank

    character(*), intent(in ) :: SLTTIntHor
    integer     , intent(in ) :: iexmin
    integer     , intent(in ) :: iexmax
    integer     , intent(in ) :: jexmin
    integer     , intent(in ) :: jexmax
    real(DP)    , intent(in ) :: x_ExtLon(iexmin:iexmax)
                              ! 緯度の拡張配列
                              ! Extended array of Lat
    real(DP)    , intent(in ) :: xyz_DPLon(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 上流点の緯度
                              ! Lat of departure point

    !
    ! local variables
    !
    integer :: iedge

    integer :: i
    integer :: j
    integer :: k


    select case (SLTTIntHor)
    case ("HQ")
      ! 方向分離型2次元変則エルミート5次補間;  2D Hermite quintic interpolation

      do k = 1, kmax
        do j = 1, jmax/2
          do i = 0, imax-1
            iedge = iexmin + 1
            if ( xyz_DPLon(i,j,k) < x_ExtLon(iedge) ) then
              call MessageNotify( 'E', module_name, &
                & 'Departure point is out of range of an extended array for longitudinal direction : Rank %d, %f < %f.', &
                & i = (/ myrank /), &
                & d = (/ xyz_DPLon(i,j,k), x_ExtLon(iedge) /) )
            end if
            iedge = iexmax - 1
            if ( xyz_DPLon(i,j,k) > x_ExtLon(iedge) ) then
              call MessageNotify( 'E', module_name, &
                & 'Departure point is out of range of an extended array for longitudinal direction : Rank %d, %f > %f.', &
                & i = (/ myrank /), &
                & d = (/ xyz_DPLon(i,j,k), x_ExtLon(iedge) /) )
            end if
          end do
        end do
      end do

    end select


  end subroutine SLTTLagIntChkDPLon

  !----------------------------------------------------------------------------
  !
  ! Checking latitude of departure point
  !
  subroutine SLTTLagIntChkDPLat(       &
    & SLTTIntHor,                      & ! (in)
    & iexmin, iexmax, jexmin, jexmax,  & ! (in)
    & y_ExtLat, xyz_DPLat              & ! (in)
    & )

    !
    ! MPI
    !
    use mpi_wrapper, only : myrank
                              ! Number of MPI rank

    character(*), intent(in ) :: SLTTIntHor
    integer     , intent(in ) :: iexmin
    integer     , intent(in ) :: iexmax
    integer     , intent(in ) :: jexmin
    integer     , intent(in ) :: jexmax
    real(DP)    , intent(in ) :: y_ExtLat(jexmin:jexmax)
                              ! 緯度の拡張配列
                              ! Extended array of Lat
    real(DP)    , intent(in ) :: xyz_DPLat(0:imax-1, 1:jmax/2, 1:kmax)
                              ! 上流点の緯度
                              ! Lat of departure point

    !
    ! local variables
    !
    integer :: jedge

    integer :: i
    integer :: j
    integer :: k


    select case (SLTTIntHor)
    case ("HQ")
      ! 方向分離型2次元変則エルミート5次補間;  2D Hermite quintic interpolation

      do k = 1, kmax
        do j = 1, jmax/2
          do i = 0, imax-1
            jedge = jexmin + 1
            if ( xyz_DPLat(i,j,k) < y_ExtLat(jedge) ) then
              call MessageNotify( 'E', module_name, &
                & 'Departure point is out of range of an extended array for latitudinal direction : Rank %d, %f < %f.', &
                & i = (/ myrank /), &
                & d = (/ xyz_DPLat(i,j,k), y_ExtLat(jedge) /) )
            end if
            jedge = jexmax - 1
            if ( xyz_DPLat(i,j,k) > y_ExtLat(jedge) ) then
              call MessageNotify( 'E', module_name, &
                & 'Departure point is out of range of an extended array for latitudinal direction : Rank %d, %f > %f.', &
                & i = (/ myrank /), &
                & d = (/ xyz_DPLat(i,j,k), y_ExtLat(jedge) /) )
            end if
          end do
        end do
      end do

    end select


  end subroutine SLTTLagIntChkDPLat

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

end module sltt_lagint
