Class sltt_lagint
In: sltt/sltt_lagint.F90

セミラグランジュ法 で用いる補間法

Interpolation methods for Semi-Lagrangian method

Note that Japanese and English are described in parallel.

セミラグランジュ法で用いる補間を演算するモジュールです. スペクトル変換・高精度補間に由来する人工的な短波を除去するために 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 :水平2次元ラグランジュ3次補間用の係数計算
SLTTLagIntCubIntHor :水平2次元ラグランジュ3次補間(上流点探索で用いる)
SLTTLagIntCubCalcFactVer :鉛直1次元ラグランジュ3次補間用の係数計算
SLTTLagIntCubIntVer :鉛直1次元ラグランジュ3次補間(上流点探索で用いる)
SLTTIrrHerIntK13 :水平2次元変則エルミート5次補間
SLTTIrrHerIntQui2DHor :水平2次元変則エルミート5次補間(コア部分)
SLTTIrrHerIntQui1DUni :1次元変則エルミート5次補間(等間隔格子)
SLTTIrrHerIntQui1DNonUni :1次元変則エルミート5次補間(不等間隔格子)
SLTTIrrHerIntQui1DUniLon :1次元変則エルミート5次補間(等間隔:経度方向専用)
SLTTHerIntCub1D :1次元エルミート3次補間
SLTTHerIntCub2D :2次元エルミート3次補間
——————— :————
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. SOLA, 4, 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. Quarterly Journal of the Royal Meteorological Society, 122(533), 1211-1226. doi:10.1002/qj.49712253310

種別型パラメタ Kind type parameter

Methods

Included Modules

dc_types dc_message gridset composition axesset sltt_const mpi_wrapper

Public Instance methods

Function :
f1 :real(DP), intent(in)
f2 :real(DP), intent(in)
g1 :real(DP), intent(in)
g2 :real(DP), intent(in)
dx :real(DP), intent(in)

エルミート3次補間 1D Hermite Cubic Interpolation

   1-----x------2

f:関数値、g:微分値、dx:点1と点2の間隔、Xi:点1と補間する点xとの間隔 f:function value, g:derivative, dx:distance between 1 and 2, Xi:distance between 1 and x

[Source]

  function SLTTHerIntCub1D(f1, f2, g1, g2, dx, Xi) result (fout)
    !エルミート3次補間
    ! 1D Hermite Cubic Interpolation
    !    1-----x------2
    ! f:関数値、g:微分値、dx:点1と点2の間隔、Xi:点1と補間する点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 :
f :real(DP), dimension(4), intent(in)
fx :real(DP), dimension(4), intent(in)
fy :real(DP), dimension(4), intent(in)
fxy :real(DP), dimension(4), intent(in)
dx :real(DP), intent(in)
dy :real(DP), intent(in)
Xix :real(DP), intent(in)

2次元エルミート3次補間 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

[Source]

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

水平方向の2次元補間。Enomoto (2008, SOLA)で提案された「スペクトルで計算した微分値を用いた双3次補間」を 発展させた方法、スペクトル微分を用いた変則エルミート5次補間を方向分離して行う。5次補間のたびに 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.

[Source]

 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 )
    ! 水平方向の2次元補間。Enomoto (2008, SOLA)で提案された「スペクトルで計算した微分値を用いた双3次補間」を
    ! 発展させた方法、スペクトル微分を用いた変則エルミート5次補間を方向分離して行う。5次補間のたびに 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
            ! 2次元補間のための配列配置(ワイド)
            ! 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次元エルミート3次補間; 2D Hermite Cubic Interpolation
          do n = 1, ncmax
          ! 2次元補間のための配列配置
          !    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次元エルミート3次補間
            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 :
f1 :real(8), intent(in)
f2 :real(8), intent(in)
f3 :real(8), intent(in)
f4 :real(8), intent(in)
g2 :real(8), intent(in)
g3 :real(8), intent(in)
dx21 :real(8), intent(in)
dx23 :real(8), intent(in)
dx24 :real(8), intent(in)

変則エルミート5次補間(2階微分不使用) 1D Hermite Quintic Interpolation 不等間隔格子の場合 non-equal separation

   1---2--x-3---4

f:関数値、g:微分値、、dx:点の間隔、Xi:点2と補間する点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)

[Source]

  function SLTTIrrHerIntQui1DNonUni(f1, f2, f3, f4, g2, g3, dx21, dx23, dx24, Xi) result (fout)
    ! 変則エルミート5次補間(2階微分不使用)
    ! 1D Hermite Quintic Interpolation
    ! 不等間隔格子の場合
    ! non-equal separation
    !    1---2--x-3---4
    ! f:関数値、g:微分値、、dx:点の間隔、Xi:点2と補間する点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))
#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 :
f1 :real(DP), intent(in)
f2 :real(DP), intent(in)
f3 :real(DP), intent(in)
f4 :real(DP), intent(in)
g2 :real(DP), intent(in)
g3 :real(DP), intent(in)
dx :real(DP), intent(in)

変則エルミート5次補間(2階微分不使用) 1D Hermite Quintic Interpolation (without 2nd derivatives) 等間隔格子の場合 equal separation

   1---2--x-3---4

f:関数値、g:微分値、、dx:点の間隔、Xi:点2と補間する点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)

[Source]

  function SLTTIrrHerIntQui1DUni(f1, f2, f3, f4, g2, g3, dx, Xi) result (fout)
    ! 変則エルミート5次補間(2階微分不使用)
    ! 1D Hermite Quintic Interpolation (without 2nd derivatives)
    ! 等間隔格子の場合
    ! equal separation
    !    1---2--x-3---4
    ! f:関数値、g:微分値、、dx:点の間隔、Xi:点2と補間する点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))
#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
Subroutine :
iexmin :integer , intent(in )
iexmax :integer , intent(in )
jexmin :integer , intent(in )
jexmax :integer , 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) :real(DP), intent(in )
xyz_ii(0:imax-1, 1:jmax/2, 1:kmax) :integer , intent(out)
xyz_jj(0:imax-1, 1:jmax/2, 1:kmax) :integer , 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) :real(DP), intent(out)

水平2次元ラグランジュ3次補間の係数計算 Calculation of factors for 2D Lagrangian cubic interpolation

[Source]

  subroutine SLTTLagIntCubCalcFactHor( iexmin, iexmax, jexmin, jexmax, x_ExtLon, y_ExtLat, xyz_MPLon, xyz_MPLat, xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify )
  ! 水平2次元ラグランジュ3次補間の係数計算
  ! 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 :
xyz_DPSigma(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyza_lcifz(0:imax-1, 1:jmax, 1:kmax, -1:2) :real(DP), intent(out)
xyz_kk(0:imax-1, 1:jmax, 1:kmax) :integer , intent(out)

鉛直1次元ラグランジュ3次補間のための係数計算 Calculation of factors for 1D Lagrangian cubic interpolation

[Source]

  subroutine SLTTLagIntCubCalcFactVer( xyz_DPSigma, xyza_lcifz, xyz_kk )
  ! 鉛直1次元ラグランジュ3次補間のための係数計算
  ! 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 :
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) :integer , 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(in )
xyz_MPArr( 0:imax-1, 1:jmax/2, 1:kmax) :real(DP), intent(out)

水平2次元ラグランジュ3次補間 2D Lagrangian cubic interpolation

[Source]

  subroutine SLTTLagIntCubIntHor( iexmin, iexmax, jexmin, jexmax, xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify, xyz_ExtArr, xyz_MPArr )
    ! 水平2次元ラグランジュ3次補間
    ! 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 :
xyz_Arr(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(in )
xyza_lcifz(0:imax-1, 1:jmax, 1:kmax, -1:2) :real(DP), intent(in )
xyz_kk(0:imax-1, 1:jmax, 1:kmax) :integer , intent(in )
xyz_ArrA(0:imax-1, 1:jmax, 1:kmax) :real(DP), intent(out)

鉛直1次元ラグランジュ3次補間 1D Lagrangian cubic interpolation

[Source]

  subroutine SLTTLagIntCubIntVer( xyz_Arr, xyza_lcifz, xyz_kk, xyz_ArrA )
  ! 鉛直1次元ラグランジュ3次補間
  ! 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

Private Instance methods

Function :
f1 :real(DP), intent(in)
f2 :real(DP), intent(in)
f3 :real(DP), intent(in)
f4 :real(DP), intent(in)
g2 :real(DP), intent(in)
g3 :real(DP), intent(in)

変則エルミート5次補間(2階微分不使用) 1D Hermite Quintic Interpolation (without 2nd derivatives) 経度方向(等間隔)専用 equal separation

   1---2--x-3---4

f:関数値、g:微分値、、dx:点の間隔、Xi:点2と補間する点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)

[Source]

  function SLTTIrrHerIntQui1DUniLon(f1, f2, f3, f4, g2, g3, Xi) result (fout)
    ! 変則エルミート5次補間(2階微分不使用)
    ! 1D Hermite Quintic Interpolation (without 2nd derivatives)
    ! 経度方向(等間隔)専用
    ! equal separation
    !    1---2--x-3---4
    ! f:関数値、g:微分値、、dx:点の間隔、Xi:点2と補間する点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))
#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 :
f :real(DP), dimension(16), intent(in)
fx :real(DP), dimension(16), intent(in)
fy :real(DP), dimension(16), intent(in)
fxy :real(DP), dimension(16), intent(in)
dy21 :real(DP), intent(in)
dy23 :real(DP), intent(in)
dy24 :real(DP), intent(in)
Xix :real(DP), intent(in)

2次元変則エルミート5次補間(2階微分不使用) 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)

[Source]

  function SLTTIrrHerIntQui2DHor(f, fx, fy, fxy, dy21, dy23, dy24, Xix, Xiy) result (fout)
    ! 2次元変則エルミート5次補間(2階微分不使用)
    ! 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 :
SLTTIntHor :character(*), intent(in )
iexmin :integer , intent(in )
iexmax :integer , intent(in )
jexmin :integer , intent(in )
jexmax :integer , intent(in )
y_ExtLat(jexmin:jexmax) :real(DP) , intent(in )
: 緯度の拡張配列 Extended array of Lat
xyz_DPLat(0:imax-1, 1:jmax/2, 1:kmax) :real(DP) , intent(in )
: 上流点の緯度 Lat of departure point

[Source]

  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
module_name
Constant :
module_name = ‘sltt_lagint :character(*), parameter
: モジュールの名称. Module name
version
Constant :
version = ’$Name: dcpam5-20150211 $’ // ’$Id: sltt_lagint.F90,v 1.4 2013/09/21 14:42:08 yot Exp $’ :character(*), parameter
: モジュールのバージョン Module version