!= 物質移流計算 (セミラグランジュ法）で用いる補間法 ! != Interpolation methods for Semi-Lagrangian method ! ! Authors:: Hiroki KASHIMURA, Yoshiyuki O. TAKAHASHI, Shin-ichi Takehiro ! Version:: ! Tag Name:: $Name: $ ! Copyright:: Copyright (C) GFD Dennou Club, 2013-2026. All rights reserved. ! License:: See COPYRIGHT[link:../../../COPYRIGHT] ! module sltt_lagint ! != セミラグランジュ法で用いる補間法 (for ISPACK3/spml2) ! != Interpolation methods for Semi-Lagrangian method (for ISPACK3/spml2) ! ! 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 :: 水平２次元ラグランジュ３次補間用の係数計算 ! 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. ! 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 ! use dc_types, only: DP, & ! 倍精度実数型. Double precision. & TOKEN ! キーワード. Keywords. ! メッセージ出力 ! Message output ! use dc_message, only: MessageNotify ! 格子点設定 ! Grid points settings ! #ifdef LIB_MPI use gridset, only: & & imax, & ! 経度格子点数. ! Number of grid points in longitude & jmax, & ! 緯度格子点数. ! Number of grid points in latitude & kmax, & ! 鉛直層数. ! Number of vertical level & imax_global ! 経度格子点数 (全球). ! Number of grid points in longitude on whole globe use ua_mpi_module_base, only : kc, pva_xvb #else use gridset, only: & & imax, & ! 経度格子点数. ! Number of grid points in longitude & jmax, & ! 緯度格子点数. ! Number of grid points in latitude & kmax, & ! 鉛直層数. ! Number of vertical level & kc => kmax,& ! 並列鉛直層数. ! Number of vertical level for parallel compuation & imax_global ! 経度格子点数 (全球). ! Number of grid points in longitude on whole globe #endif 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 2026/05/12 12:00:00 takepiro Exp $' ! モジュールのバージョン ! Module version contains !-------------------------------------------------------------------------------------- subroutine SLTTLagIntCubCalcFactHor( & & 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 ) :: jexmin integer , intent(in ) :: jexmax real(DP), intent(in ) :: x_ExtLon(0:imax_global-1) real(DP), intent(in ) :: y_ExtLat(jexmin:jexmax) real(DP), intent(in ) :: xyz_MPLon(0:imax_global-1, 1:jmax, 1:kc) real(DP), intent(in ) :: xyz_MPLat(0:imax_global-1, 1:jmax, 1:kc) integer , intent(out) :: xyz_ii(0:imax_global-1, 1:jmax, 1:kc) integer , intent(out) :: xyz_jj(0:imax_global-1, 1:jmax, 1:kc) real(DP), intent(out) :: xyza_lcifx(0:imax_global-1, 1:jmax, 1:kc, -1:2) real(DP), intent(out) :: xyza_lcify(0:imax_global-1, 1:jmax, 1:kc, -1:2) ! ! local variables ! integer :: ii integer :: jj integer :: i, j, k, j2 real(DP) :: ExtLonm1, ExtLonp1, ExtLonp2 real(DP) :: dlon if ( imax_global == 1 ) then dlon = 0.0D0 else dlon = PIx2/imax_global endif xyz_ii = int( xyz_MPLon / ( PIx2 / imax_global ) ) do k = 1, kc do j = 1, jmax do i = 0, imax_global-1 ! 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 end do end do end do do k = 1, kc do j = 1, jmax do i = 0, imax_global-1 jj = xyz_jj(i,j,k) if ( jj-1 < jexmin ) then call MessageNotify( 'E', module_name, & & 'Latitudinal point for interporation factor calculation ' & & // 'is 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 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, kc do j = 1, jmax do i = 0, imax_global-1 ii = xyz_ii(i,j,k) jj = xyz_jj(i,j,k) ExtLonm1 = x_ExtLon(ii) - dlon ExtLonp1 = x_ExtLon(ii) + dlon ExtLonp2 = x_ExtLon(ii) + dlon*2 xyza_lcifx(i,j,k,-1) = & & ( xyz_MPLon(i,j,k) - x_ExtLon(ii ) ) & & * ( xyz_MPLon(i,j,k) - ExtLonp1 ) & & * ( xyz_MPLon(i,j,k) - ExtLonp2 ) & & * (-InvDeltaLon**3)/6.0_DP ! economical xyza_lcifx(i,j,k, 0) = & & ( xyz_MPLon(i,j,k) - ExtLonm1 ) & & * ( xyz_MPLon(i,j,k) - ExtLonp1 ) & & * ( xyz_MPLon(i,j,k) - ExtLonp2 ) & & * (0.5_DP*InvDeltaLon**3) ! economical xyza_lcifx(i,j,k, 1) = & & ( xyz_MPLon(i,j,k) - ExtLonm1 ) & & * ( xyz_MPLon(i,j,k) - x_ExtLon(ii ) ) & & * ( xyz_MPLon(i,j,k) - ExtLonp2 ) & & * (-0.5_DP*InvDeltaLon**3) ! economical xyza_lcifx(i,j,k, 2) = & & ( xyz_MPLon(i,j,k) - ExtLonm1 ) & & * ( xyz_MPLon(i,j,k) - x_ExtLon(ii ) ) & & * ( xyz_MPLon(i,j,k) - ExtLonp1 ) & & * (InvDeltaLon**3)/6.0_DP ! economical !!$ xyza_lcifx(i,j,k,-1) = & !!$ & ( xyz_MPLon(i,j,k) - x_ExtLon(ii ) ) & !!$ & * ( xyz_MPLon(i,j,k) - x_ExtLon(iip1) ) & !!$ & * ( xyz_MPLon(i,j,k) - x_ExtLon(iip2) ) & !!$ & * (-InvDeltaLon**3)/6.0_DP ! economical !!$ ! & / ( ( x_ExtLon(iim1) - x_ExtLon(ii ) ) & !!$ ! & * ( x_ExtLon(iim1) - x_ExtLon(iip1) ) & !!$ ! & * ( x_ExtLon(iim1) - x_ExtLon(iip2) ) ) !!$ xyza_lcifx(i,j,k, 0) = & !!$ & ( xyz_MPLon(i,j,k) - x_ExtLon(iim1) ) & !!$ & * ( xyz_MPLon(i,j,k) - x_ExtLon(iip1) ) & !!$ & * ( xyz_MPLon(i,j,k) - x_ExtLon(iip2) ) & !!$ & * (0.5_DP*InvDeltaLon**3) ! economical !!$ ! & / ( ( x_ExtLon(ii ) - x_ExtLon(iim1) ) & !!$ ! & * ( x_ExtLon(ii ) - x_ExtLon(iip1) ) & !!$ ! & * ( x_ExtLon(ii ) - x_ExtLon(iip2) ) ) !!$ xyza_lcifx(i,j,k, 1) = & !!$ & ( xyz_MPLon(i,j,k) - x_ExtLon(iim1) ) & !!$ & * ( xyz_MPLon(i,j,k) - x_ExtLon(ii ) ) & !!$ & * ( xyz_MPLon(i,j,k) - x_ExtLon(iip2) ) & !!$ & * (-0.5_DP*InvDeltaLon**3) ! economical !!$ ! & / ( ( x_ExtLon(iip1) - x_ExtLon(iim1) ) & !!$ ! & * ( x_ExtLon(iip1) - x_ExtLon(ii ) ) & !!$ ! & * ( x_ExtLon(iip1) - x_ExtLon(iip2) ) ) !!$ xyza_lcifx(i,j,k, 2) = & !!$ & ( xyz_MPLon(i,j,k) - x_ExtLon(iim1) ) & !!$ & * ( xyz_MPLon(i,j,k) - x_ExtLon(ii ) ) & !!$ & * ( xyz_MPLon(i,j,k) - x_ExtLon(iip1) ) & !!$ & * (InvDeltaLon**3)/6.0_DP ! economical !!$ ! & / ( ( x_ExtLon(iip2) - x_ExtLon(iim1) ) & !!$ ! & * ( x_ExtLon(iip2) - x_ExtLon(ii ) ) & !!$ ! & * ( x_ExtLon(iip2) - x_ExtLon(iip1) ) ) !!$ 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( & & jexmin, jexmax, & ! (in) & xyz_ii, xyz_jj, xyza_lcifx, xyza_lcify, xyz_ExtArr, & ! (in) & xyz_MPArr & ! (out) & ) ! 水平２次元ラグランジュ３次補間 ! 2D Lagrangian cubic interpolation integer , intent(in ) :: jexmin integer , intent(in ) :: jexmax integer , intent(in ) :: xyz_ii(0:imax_global-1, 1:jmax, 1:kc) integer , intent(in ) :: xyz_jj(0:imax_global-1, 1:jmax, 1:kc) real(DP), intent(in ) :: xyza_lcifx(0:imax_global-1, 1:jmax, 1:kc, -1:2) real(DP), intent(in ) :: xyza_lcify(0:imax_global-1, 1:jmax, 1:kc, -1:2) real(DP), intent(in ) :: xyz_ExtArr(0:imax_global-1, jexmin:jexmax, 1:kc) real(DP), intent(out) :: xyz_MPArr (0:imax_global-1, 1:jmax, 1:kc) ! ! local variables ! integer :: ii integer :: jj integer :: kk integer :: i, j, k integer :: iim1, iip1, iip2 do k = 1, kc do j = 1, jmax do i = 0, imax_global-1 ii = xyz_ii(i,j,k) jj = xyz_jj(i,j,k) kk = k iim1 = modulo(ii-1,imax_global) iip1 = modulo(ii+1,imax_global) iip2 = modulo(ii+2,imax_global) xyz_MPArr(i,j,k) = & & xyza_lcify(i,j,k,-1) & & * ( xyza_lcifx(i,j,k,-1) * xyz_ExtArr(iim1,jj-1,kk) & & + xyza_lcifx(i,j,k, 0) * xyz_ExtArr(ii ,jj-1,kk) & & + xyza_lcifx(i,j,k, 1) * xyz_ExtArr(iip1,jj-1,kk) & & + xyza_lcifx(i,j,k, 2) * xyz_ExtArr(iip2,jj-1,kk) )& & + xyza_lcify(i,j,k, 0) & & * ( xyza_lcifx(i,j,k,-1) * xyz_ExtArr(iim1,jj ,kk) & & + xyza_lcifx(i,j,k, 0) * xyz_ExtArr(ii ,jj ,kk) & & + xyza_lcifx(i,j,k, 1) * xyz_ExtArr(iip1,jj ,kk) & & + xyza_lcifx(i,j,k, 2) * xyz_ExtArr(iip2,jj ,kk) )& & + xyza_lcify(i,j,k, 1) & & * ( xyza_lcifx(i,j,k,-1) * xyz_ExtArr(iim1,jj+1,kk) & & + xyza_lcifx(i,j,k, 0) * xyz_ExtArr(ii ,jj+1,kk) & & + xyza_lcifx(i,j,k, 1) * xyz_ExtArr(iip1,jj+1,kk) & & + xyza_lcifx(i,j,k, 2) * xyz_ExtArr(iip2,jj+1,kk) )& & + xyza_lcify(i,j,k, 2) & & * ( xyza_lcifx(i,j,k,-1) * xyz_ExtArr(iim1,jj+2,kk) & & + xyza_lcifx(i,j,k, 0) * xyz_ExtArr(ii ,jj+2,kk) & & + xyza_lcifx(i,j,k, 1) * xyz_ExtArr(iip1,jj+2,kk) & & + xyza_lcifx(i,j,k, 2) * xyz_ExtArr(iip2,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( & & 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 ) :: jexmin integer , intent(in ) :: jexmax real(DP), intent(in ) :: x_ExtLon(0:imax_global-1) ! 経度の拡張配列 ! Extended array of Lon real(DP), intent(in ) :: y_ExtLat(jexmin:jexmax) ! 緯度の拡張配列 ! Extended array of Lat real(DP), intent(in ) :: xyz_DPLon(0:imax_global-1, 1:jmax, 1:kc) ! 上流点の経度 ! Lon at departure point real(DP), intent(in ) :: xyz_DPLat(0:imax_global-1, 1:jmax, 1:kc) ! 上流点の緯度 ! Lat at departure point real(DP), intent(in ) :: xyzf_ExtQMixB(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax) ! 物質混合比の拡張配列 ! Extended array of mix ration of tracers real(DP), intent(out) :: xyzf_QMixMinA(0:imax_global-1, 1:jmax, 1:kc, 1:ncmax) ! 上流点を囲む 4 格子点の最小値 ! Minimum mixing ratio of tracers at 4 grid points ! surrounding a departure point real(DP), intent(out) :: xyzf_QMixMaxA(0:imax_global-1, 1:jmax, 1:kc, 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:: iis 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, kc do j = 1, jmax do i = 0, imax_global-1 ! 上流点を囲む4点を探す ! Determine four grid points around the departure point ! isten = int(xyz_DPLon(i,j,k)*InvDeltaLon) do ii = 0, imax_global-1 if (x_ExtLon(ii) > xyz_DPLon(i, j, k)) then isten = ii-1 exit endif isten = imax_global-1 enddo do jj = jexmin, jexmax if (y_ExtLat(jj) > xyz_DPLat(i, j, k)) then jsten = jj-1 exit endif enddo ! 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 iis = modulo(isten+ii,imax_global) xy_Z(ii,jj) = xyzf_ExtQMixB(iis,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( & & 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 ) :: jexmin integer , intent(in ) :: jexmax real(DP), intent(in ) :: x_ExtLon(0:imax_global-1) ! 経度の拡張配列 ! Extended array of Lon real(DP), intent(in ) :: y_ExtLat(jexmin:jexmax) ! 緯度の拡張配列 ! Extended array of Lat real(DP), intent(in ) :: xyz_DPLon(0:imax_global-1, 1:jmax, 1:kc) ! 上流点の経度 ! Lon at departure point real(DP), intent(in ) :: xyz_DPLat(0:imax_global-1, 1:jmax, 1:kc) ! 上流点の緯度 ! Lat at departure point real(DP), intent(in ) :: xyzf_ExtQMixB(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax) ! 物質混合比の拡張配列 ! Extended array of mix ration of tracers real(DP), intent(out) :: xyzf_QMixA (0:imax-1, 1:jmax, 1:kmax, 1:ncmax) ! 次ステップの物質混合比 ! Mix ration of tracers at next time-step real(DP) :: xyzf_ExtQMixA (0:imax_global-1, 1:jmax, 1:kc, 1:ncmax) ! 次ステップの物質混合比 ! Mix ration of tracers at next time-step real(DP) :: Lon_isten ! 基準経度座標 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 格子点上の混合比を格納する作v業変数 ! 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 integer :: isp1 ! 配列配置のための作業変数 do k = 1, kc do j = 1, jmax do i = 0, imax_global-1 ! 上流点を囲む4点を探す ! Determine four grid points around the departure point ! isten = int(xyz_DPLon(i,j,k)*InvDeltaLon) do ii = 0, imax_global-1 if (x_ExtLon(ii) > xyz_DPLon(i, j, k)) then isten = ii-1 Lon_isten = x_ExtLon(ii)-DeltaLon exit endif isten = imax_global-1 Lon_isten = x_ExtLon(imax_global-1) enddo 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 isp1 = modulo(isten+ii,imax_global) xy_Z(ii,jj) = xyzf_ExtQMixB(isp1,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) - Lon_isten ) + xy_Z(0,jj) !!$ 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_ExtQMixA(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 do n = 1, ncmax #ifdef LIB_MPI xyzf_QMixA(:,:,:,n) = pva_xvb(xyzf_ExtQMixA(:,1:jmax,:,n)) #else xyzf_QMixA(:,:,:,n) = xyzf_ExtQMixA(:,1:jmax,:,n) #endif end do end subroutine SLTTIrrLinInt !-------------------------------------------------------------------------------------- subroutine SLTTIrrHerIntK13( & & jexmin, jexmax, & ! (in) & x_ExtLon, y_ExtLat, xyz_DPLon, xyz_DPLat, xyz_ExtQMixB, & ! (in) & xyz_ExtQMixB_dlon, xyz_ExtQMixB_dlat, xyz_ExtQMixB_dlonlat,& ! (in) & 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 ) :: jexmin integer , intent(in ) :: jexmax real(DP), intent(in ) :: x_ExtLon(0:imax_global-1) ! 経度の拡張配列 ! Extended array of Lon real(DP), intent(in ) :: y_ExtLat(jexmin:jexmax) ! 緯度の拡張配列 ! Extended array of Lat real(DP), intent(in ) :: xyz_DPLon(0:imax_global-1, 1:jmax, 1:kc) ! 上流点の経度 ! Lon at departure point real(DP), intent(in ) :: xyz_DPLat(0:imax_global-1, 1:jmax, 1:kc) ! 上流点の緯度 ! Lat at departure point real(DP), intent(in ) :: xyz_ExtQMixB(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax) ! 物質混合比の拡張配列 ! Extended array of mix ration of tracers real(DP), intent(in) :: xyz_ExtQMixB_dlon(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax) ! 物質混合比の経度微分の拡張配列 ! Extended array of zonal derivative of the mix ration real(DP), intent(in) :: xyz_ExtQMixB_dlat(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax) ! 物質混合比の緯度微分の拡張配列 ! Extended array of meridional derivative of the mix ration real(DP), intent(in) :: xyz_ExtQMixB_dlonlat(0:imax_global-1, jexmin:jexmax, 1:kc, 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, 1:kmax, 1:ncmax) ! 次ステップの物質混合比 ! Mix ration of tracers at next time-step real(DP) :: xyz_ExtQMixA (0:imax_global-1, 1:jmax, 1:kc, 1:ncmax) ! 次ステップの物質混合比 !---fxx, fyy, fxxy, fxyy, fxxyy ! real(DP), intent(inout) :: xyz_ExtQMixB_dlon2(-2+0:imax-1+3, -jew+1:jmax+jew, 1:kmax) ! real(DP), intent(inout) :: xyz_ExtQMixB_dlat2(-2+0:imax-1+3, -jew+1:jmax+jew, 1:kmax) ! real(DP), intent(inout) :: xyz_ExtQMixB_dlon2lat(-2+0:imax-1+3, -jew+1:jmax+jew, 1:kmax) ! real(DP), intent(inout) :: xyz_ExtQMixB_dlonlat2(-2+0:imax-1+3, -jew+1:jmax+jew, 1:kmax) ! real(DP), intent(inout) :: xyz_ExtQMixB_dlon2lat2(-2+0:imax-1+3, -jew+1:jmax+jew, 1:kmax) real(DP) :: Lon_isten ! 基準経度座標 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:: iise 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 integer :: isp1 ! 配列配置のための作業変数 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) !!$ & jexmin, jexmax, & ! (in) !!$ & x_ExtLon, xyz_DPLon & ! (in) !!$ & ) !!$ ! Check whether a latitude of departure point is within an extended array !!$ call SLTTLagIntChkDPLat( & !!$ & SLTTIntHor, & ! (in) !!$ & jexmin, jexmax, & ! (in) !!$ & y_ExtLat, xyz_DPLat & ! (in) !!$ & ) do k = 1, kc do j = 1, jmax do i = 0, imax_global-1 ! 上流点を囲む4点を探す ! Determine four grid points around the departure point do ii = 0, imax_global-1 if (x_ExtLon(ii) > xyz_DPLon(i, j, k)) then isten = ii-1 Lon_isten = x_ExtLon(ii)-DeltaLon exit endif isten = imax_global-1 Lon_isten = x_ExtLon(imax_global-1) enddo do jj = jexmin, jexmax if (y_ExtLat(jj) > xyz_DPLat(i, j, k)) then jsten = jj-1 exit endif enddo ! 水平方向の補間方法の選択 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 iise = modulo(isten+ii,imax_global) a_z(ii+num) = xyz_ExtQMixB(iise, jsten+jj, k, n) a_zx(ii+num) = xyz_ExtQMixB_dlon(iise, jsten+jj, k, n) a_zy(ii+num) = xyz_ExtQMixB_dlat(iise, jsten+jj, k, n) a_zxy(ii+num) = xyz_ExtQMixB_dlonlat(iise, jsten+jj, k, n) enddo enddo ! 方向分離型2次元変則エルミート5次補間 ! 2D Hermite quintic interpolation xyz_ExtQMixA(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)-Lon_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 isp1 = modulo(isten+1,imax_global) 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(isp1, jsten, k, n) a_zx(2) = xyz_ExtQMixB_dlon(isp1, jsten, k, n) a_zy(2) = xyz_ExtQMixB_dlat(isp1, jsten, k, n) a_zxy(2) = xyz_ExtQMixB_dlonlat(isp1, jsten, k, n) ! a_zxx(2) = xyz_ExtQMixB_dlon2(isp1, jsten, k, n) ! a_zyy(2) = xyz_ExtQMixB_dlat2(isp1, jsten, k, n) ! a_zxxy(2) = xyz_ExtQMixB_dlon2lat(isp1, jsten, k, n) ! a_zxyy(2) = xyz_ExtQMixB_dlonlat2(isp1, jsten, k, n) ! a_zxxyy(2) = xyz_ExtQMixB_dlon2lat2(isp1, jsten, k, n) a_z(3) = xyz_ExtQMixB(isp1, jsten+1, k, n) a_zx(3) = xyz_ExtQMixB_dlon(isp1, jsten+1, k, n) a_zy(3) = xyz_ExtQMixB_dlat(isp1, jsten+1, k, n) a_zxy(3) = xyz_ExtQMixB_dlonlat(isp1, jsten+1, k, n) ! a_zxx(3) = xyz_ExtQMixB_dlon2(isp1, jsten+1, k, n) ! a_zyy(3) = xyz_ExtQMixB_dlat2(isp1, jsten+1, k, n) ! a_zxxy(3) = xyz_ExtQMixB_dlon2lat(isp1, jsten+1, k, n) ! a_zxyy(3) = xyz_ExtQMixB_dlonlat2(isp1, jsten+1, k, n) ! a_zxxyy(3) = xyz_ExtQMixB_dlon2lat2(isp1, 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_ExtQMixA(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 IN NAMELIST.' ) end select enddo enddo enddo do n = 1, ncmax #ifdef LIB_MPI xyz_QMixA(:,:,:,n) = pva_xvb(xyz_ExtQMixA(:,1:jmax,:,n)) #else xyz_QMixA(:,:,:,n) = xyz_ExtQMixA(:,1:jmax,:,n) #endif end do 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) & 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 ) :: jexmin integer , intent(in ) :: jexmax real(DP) , intent(in ) :: x_ExtLon(0:imax_global-1) ! 緯度の拡張配列 ! Extended array of Lat real(DP) , intent(in ) :: xyz_DPLon(0:imax_global-1, 1:jmax, 1:kc) ! 上流点の緯度 ! 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, kc do j = 1, jmax do i = 0, imax_global-1 iedge = 0 !!$ iedge = imax_global-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 = 0 iedge = imax_global-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) & 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 ) :: 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_global-1, 1:jmax, 1:kc) ! 上流点の緯度 ! 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, kc do j = 1, jmax do i = 0, imax_global-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