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 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

 subroutine SLTTIrrHerIntK13( iexmin, iexmax, jexmin, jexmax, x_ExtLon, y_ExtLat, xyz_DPLon, xyz_DPLat, xyz_ExtQMixB, xyz_ExtQMixB_dlon, xyz_ExtQMixB_dlat, xyz_ExtQMixB_dlonlat, SLTTIntHor, xyz_QMixA )
    ! 水平方向の２次元補間。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
    real(DP) :: InvDeltaLat   ! 1.0/deltalat
    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 latitude of departure point is within an extended array
    call SLTTLagIntChkDPLat( SLTTIntHor, iexmin, iexmax, jexmin, jexmax, y_ExtLat, xyz_DPLat )


 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
          write( 6, * ) "ERROR : GIVE CORRECT KEYWORD FOR <SLTTIntHor> IN NAMELIST"
          stop
                    
      end select



    enddo
  enddo
enddo

end subroutine SLTTIrrHerIntK13

  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 
    if ((f2 - f1)*(f4 - f3)>=0.0_DP) then
    fout = min(max(f2,f3), max(min(f2,f3), fout))
    endif

    end function SLTTIrrHerIntQui1DNonUni 

  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     
    if ((f2 - f1)*(f4 - f3)>=0.0_DP) then
      fout = min(max(f2,f3), max(min(f2,f3), fout))
    endif

  end function SLTTIrrHerIntQui1DUni

  subroutine SLTTLagIntCubCalcFactHor( iexmin, iexmax, jexmin, jexmax, x_ExtLon, y_ExtLat, xyz_MPLon, xyz_MPLat, xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify )
  ! 水平２次元ラグランジュ３次補間の係数計算
  ! 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 ) )

    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, 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

!             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

!#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_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 SLTTLagIntCubCalcFactVer( xyz_DPSigma, xyza_lcifz, xyz_kk )
  ! 鉛直１次元ラグランジュ３次補間のための係数計算
  ! 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 SLTTLagIntCubIntHor( iexmin, iexmax, jexmin, jexmax, xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify, xyz_ExtArr, xyz_MPArr )
    ! 水平２次元ラグランジュ３次補間
    ! 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 SLTTLagIntCubIntVer( xyz_Arr, xyza_lcifz, xyz_kk, xyz_ArrA )
  ! 鉛直１次元ラグランジュ３次補間
  ! 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

  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     
    if ((f2 - f1)*(f4 - f3)>=0.0_DP) then
      fout = min(max(f2,f3), max(min(f2,f3), fout))
    endif

  end function SLTTIrrHerIntQui1DUniLon

  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

  subroutine SLTTLagIntChkDPLat( SLTTIntHor, iexmin, iexmax, jexmin, jexmax, y_ExtLat, xyz_DPLat )

    !
    ! 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 : 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 : 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