!= 物質移流過程 (セミラグランジュ法)
!
!= Semi-Lagrangian Tracer Transport scheme
!
! Authors:: Hiroki KASHIMURA, Yoshiyuki O. TAKAHASHI, Shin-ichi Takehiro
! Version:: $Id: sltt.F90,v 1.3 2026/03/02 22:00:00 takepiro Exp $
! Tag Name:: $Name: $
! Copyright:: Copyright (C) GFD Dennou Club, 2013-2026. All rights reserved.
! License:: See COPYRIGHT[link:../../../COPYRIGHT]
!
module sltt
!
!= 物質移流 (セミラグランジュ法, Enomoto (2008) modified, ISPACK3/spml2 用)
!
!= Tracer Transport (Semi-Lagrangian method, Enomoto (2008) modified)
!
! Note that Japanese and English are described in parallel.
!
! 物質移流を非保存型のセミラグランジュ法で演算するモジュールです.
! 上流点探索には Williamson and Rasch (1989, MWR) を
! 補間には Enomoto (2008) を応用した方法を用いています。
! すなわちスペクトルから求めた１階微分の値を利用した５次精度の変則エルミート補間です。
! 非負を保証するために arcsine 変換フィルタを用いています。
! スペクトル変換・高精度補間に由来する人工的な短波を除去するために Sun et al. (1996) の
! 単調フィルタを応用したものを部分的に用いている。
!
! This is a tracer transport module. Semi-Lagrangian method (Enomoto 2008 modified, for ISPACK3/spml2)
! Arcsine transformation filter is used to avoid negative values.
! Monotonicity filter (Sun et al 1996) is partly used.
!
!== Procedures List
!
! SLTTMain :: 移流計算
! SLTTInit :: 初期化
! SLTTTest :: 移流テスト用
! --------------------- :: ------------
! SLTTMain :: Main subroutine for SLTT
! SLTTInit :: Initialization for SLTT
! SLTTTest :: Generate velocity for SLTT Test
!
!== NAMELIST
!
! NAMELIST#
!
!== References
!
! * Kashimura, H., T. Enomoto, Y. O. Takahashi, 2013:
! Non-negative filter using arcsine transformation for tracer advection with semi-Lagrangian scheme.
! NCTAM, 62.
!
! * Enomoto, T., 2008:
! Bicubic Interpolation with Spectral Derivatives.
! SOLA, 4, 5-8. doi:10.2151/sola.2008-002
!
! * Williamson, D. L., and Rasch, P. J., 1989:
! Two-dimensional semi-Lagrangian transport with shape-preserving interpolation.
! Mon. Wea. Rev., 117, 102-129.
!
! * 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
! モジュール引用 ; USE statements
!
! 種別型パラメタ
! Kind type parameter
!
use dc_types, only: DP, & ! 倍精度実数型. Double precision.
& TOKEN ! キーワード. Keywords.
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
!
! MPI
!
use mpi_wrapper, only : MPIWrapperFindMaxVal, myrank
! 時刻管理
! Time control
!
use timeset, only: &
& DelTime
! 格子点設定
! 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
& lmax, & ! スペクトルデータの配列サイズ
! Size of array for spectral data
& imax_global ! 経度格子点数 (全球).
! Number of grid points in longitude on whole globe
use ua_mpi_module_base, only : kc
#else
use gridset, only: &
& imax, & ! 経度格子点数.
! Number of grid points in longitude
& jmax, & ! 緯度格子点数.
! Number of grid points in latitude
& kmax, & ! 鉛直層数.
! Number of vertical level
& kc => kmax,& ! 並列鉛直層数.
! Number of vertical level for parallel compuation
& lmax, & ! スペクトルデータの配列サイズ
! Size of array for spectral data
& imax_global ! 経度格子点数 (全球).
! Number of grid points in longitude on whole globe
#endif
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: &
& ncmax, &
! 成分の数
! Number of composition
& CompositionInqFlagAdv
! 質量の補正
! Mass fixer
!
use mass_fixer, only: &
& MassFixerBC02, MassFixerBC02Layer, MassFixerBC02Column, &
& MassFixer, MassFixerR95, MassFixerWO94, MassFixerColumn!, MassFixerLayer
! 宣言文 ; Declaration statements
!
implicit none
private
! 公開手続き
! Public procedure
!
public :: SLTTInit
public :: SLTTMain
public :: SLTTHorAdv
public :: SLTTVerAdv
public :: SLTTCheckInit
! 公開変数
! Public variables
!
! 非公開変数
! Private variables
!
logical, save :: sltt_inited = .false.
! 初期設定フラグ.
! Initialization flag
logical, save :: sltt_check_inited = .false.
! 初期設定フラグ.
! Initialization flag
real(DP) , save, allocatable :: x_SinLon(:)
! $\sin\lambda_S$
real(DP) , save, allocatable :: x_CosLon(:)
! $\cos\lambda_S$
real(DP) , save, allocatable :: y_SinLat(:)
! $\sin\varphai_S$
real(DP) , save, allocatable :: y_CosLat(:)
! $\cos\varphai_S$
real(DP) , save, allocatable :: x_ExtLon(:)
! $ x_LonSの拡張配列。
!Extended array of x_LonS.
real(DP) , save, allocatable :: y_ExtLat(:)
! $ x_LatSの拡張配列。
!Extended array of x_LatS.
logical, save :: FlagSLTTArcsineHor
logical, save :: FlagSLTTArcsineVer
! Arcsine変換の非負フィルタフラグ
! Flag for non-negative filter using arcsine trasformation
real(DP), save :: SLTTArcSineFactor
character(TOKEN), save :: SLTTIntHor
! 水平方向の補間方法を指定するキーワード
! Keyword for Interpolation Method for Horizontal direction
character(TOKEN), save :: SLTTIntVer
! 鉛直方向の補間方法を指定するキーワード
! Keyword for Interpolation Method for Vertical direction
character(*), parameter:: module_name = 'sltt'
! モジュールの名称.
! Module name
character(*), parameter:: version = &
& '$Name: $' // &
& '$Id: sltt.F90,v 1.8 2026/03/02 22:00:00 takepiro Exp $'
! モジュールのバージョン
! Module version
!--------------------------------------------------------------------------------------
contains
!--------------------------------------------------------------------------------------
subroutine SLTTMain( &
& xyr_PressB, xyr_PressA, & !(in )
& xyz_UN, xyz_VN, xyr_SigDotN, & !(in )
& xyzf_DQMixDtPhy, & !(in )
& xyzf_QMixB, & !(in )
& xyzf_QMixA & !(out)
& )
! セミラグランジュ法による物質移流計算を行う。
! Calculates tracer transports by Semi-Lagrangian method
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
use timeset , only : &
& TimeN, &
& DelTime
! $\Delta t$
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: &
& ncmax, &
! 成分の数
! Number of composition
& a_QMixName, &
! 成分の変数名
! Name of variables for composition
& CompositionInqFlagAdv
real(DP), intent(in ) :: xyr_PressB(0:imax-1, 1:jmax, 0:kmax)
!
! Pressure at current time step
real(DP), intent(in ) :: xyr_PressA(0:imax-1, 1:jmax, 0:kmax)
!
! Pressure at next time step
real(DP), intent(in ) :: xyz_UN (0:imax-1, 1:jmax, 1:kmax)
! 東西風速
! Zonal Wind
real(DP), intent(in ) :: xyz_VN (0:imax-1, 1:jmax, 1:kmax)
! 南北風速
! Meridional Wind
real(DP), intent(in ) :: xyr_SigDotN(0:imax-1, 1:jmax, 0:kmax)
! 鉛直流速（SigmaDot）
real(DP), intent(in ):: xyzf_DQMixDtPhy(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! $ \left(\DP{q}{t}\right)^{phy} $ .
! 外力項 (物理過程) による比湿変化.
! Temperature tendency by external force terms (physical processes)
real(DP), intent(in ) :: xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 物質混合比
! Mix ratio of the tracers
real(DP), intent(out) :: xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 物質混合比
! Mix ratio of the tracers
! 作業変数
! Work variables
!
real(DP) :: f_QMixMax(1:ncmax)
! 各物質混合比の最大値
! Maximum of each mix ratio of the tracers
real(DP) :: f_QMixProcMax(1:ncmax)
! 各物質混合比のプロセス内最大値
! Maximum of each mix ratio of the tracers in each process
real(DP) :: f_QMixLinMax(1:ncmax)
real(DP) :: f_QMixLinProcMax(1:ncmax)
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
real(DP) :: xyz_UTest (0:imax-1, 1:jmax, 1:kmax)
! 東西風速（テスト用）
! Zonal Wind (for test)
real(DP) :: xyz_VTest (0:imax-1, 1:jmax, 1:kmax)
! 南北風速（テスト用）
! Meridional Wind (for test)
real(DP) :: xyr_SigDotTest(0:imax-1, 1:jmax, 0:kmax)
! 鉛直流速（テスト用）;SigmaDot (for test)
real(DP) :: xyzf_QMixSave(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP) :: xyzf_QMixLinATentative(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP) :: xyzf_QMixLinA (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! Variables for monotone limiter
real(DP) :: xyzf_QMixMinA (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP) :: xyzf_QMixMaxA (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP) :: xyzf_QMixSaveMassFix (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP) :: xyzf_DQMixDtHorMassFix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP) :: xyzf_DQMixDtVerMassFix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP) :: xyzf_DQMixDtTotMassFix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! セミラグランジュ法による物質移流計算
! Semi-Lagrangian method for tracer transport
xyzf_QMixA = xyzf_QMixB + xyzf_DQMixDtPhy * 2.0_DP * DelTime
! Save a variable for mass fixer
xyzf_QMixSave = xyzf_QMixA
! Mass fixer
! Constituents
!
call MassFixerColumn( &
& xyr_PressB, & ! (in)
& xyzf_QMixA, & ! (inout)
& xyr_PressRef = xyr_PressB, & ! (in) optional
& xyzf_QMixRef = xyzf_QMixSave & ! (in) optional
& )
! Save a variable for mass fixer
xyzf_QMixSave = xyzf_QMixA
! Variable for linear interpolation
xyzf_QMixLinA = xyzf_QMixA
if ( FlagSLTTArcsineHor ) then
! 非負を保証するための arcsine変換フィルタ
! Arcsine transformation for non-negative filter
do n = 1, ncmax
f_QMixProcMax(n) = maxval( xyzf_QMixA(:,:,:,n) )
end do
call MPIWrapperFindMaxVal( &
& ncmax, f_QMixProcMax, & ! (in)
& f_QMixMax & ! (out)
& )
f_QMixMax = f_QMixMax * SLTTArcSineFactor + 1.0e-14_DP
do n = 1, ncmax
xyzf_QMixA(:,:,:,n) = &
& 0.5_DP*(asin(2.0_DP*xyzf_QMixA(:,:,:,n)/f_QMixMax(n) - 1.0_DP))
end do
! arcsine transformed variable is used for linear interpolation too
xyzf_QMixLinA = xyzf_QMixA
f_QMixLinProcMax = f_QMixProcMax
f_QMixLinMax = f_QMixMax
end if
! 水平セミラグ
! Horizontal
xyzf_QMixA = SLTTHorAdv( xyzf_QMixA, xyz_UN, xyz_VN, & ! (in)
& xyzf_QMixLinA = xyzf_QMixLinA, & ! (inout) optional
& xyzf_QMixMinA = xyzf_QMixMinA, & ! (out) optional
& xyzf_QMixMaxA = xyzf_QMixMaxA ) ! (out) optional
! Monotonic filter
! see Diamantakis and Flemming (2014) for BS limiter
! but limiter is applied separately in horizontal and vertical directions
#ifdef SLTT2D1DMONOTONIC
xyzf_QMixA = max( min( xyzf_QMixA, xyzf_QMixMaxA ), xyzf_QMixMinA )
#endif
!==================================================
! Calculation in a case in which mass fixer applied in horizontal and
! vertical directions separately
!
!!$ if (FlagSLTTArcsine) then
!!$ ! 非負を保証するための arcsine変換フィルタ（逆変換）
!!$ ! Arcsine transformation for non-negative filter
!!$ do n = 1, ncmax
!!$ xyzf_QMixA(:,:,:,n) = &
!!$ & f_QMixMax(n)*(0.5_DP)*(sin(2.0_DP*xyzf_QMixA(:,:,:,n))+1.0_DP)
!!$! xyzf_QMixLinA(:,:,:,n) = &
!!$! & f_QMixMax(n)*(0.5_DP)*(sin(2.0_DP*xyzf_QMixLinA(:,:,:,n))+1.0_DP)
!!$ enddo
!!$ endif
!!$ !
!!$ xyzf_QMixSaveMassFix = xyzf_QMixA
!!$ !
!!$ call MassFixerBC02Layer( &
!!$ & xyr_PressA, & ! (in)
!!$ & xyzf_QMixA, & ! (inout)
!!$ & xyzf_QMixLinA, & ! (in)
!!$ & xyr_PressB, & ! (in)
!!$ & xyzf_QMixSave & ! (in)
!!$ & )
!!$ !
!!$ xyzf_DQMixDtHorMassFix = &
!!$ & ( xyzf_QMixA - xyzf_QMixSaveMassFix ) / ( 2.0_DP * DelTime )
!!$ !
!!$ ! Save a variable for mass fixer
!!$ xyzf_QMixSave = xyzf_QMixA
!!$ !
!!$ ! Variable for linear interpolation
!!$ xyzf_QMixLinATentative = xyzf_QMixA
!!$ !
!!$ if (FlagSLTTArcsine) then
!!$ ! 非負を保証するための arcsine変換フィルタ
!!$ ! Arcsine transformation for non-negative filter
!!$ !
!!$ do n = 1, ncmax
!!$ f_QMixProcMax(n) = maxval( xyzf_QMixA(:,:,:,n) )
!!$ end do
!!$ call MPIWrapperFindMaxVal( &
!!$ & ncmax, f_QMixProcMax, & ! (in)
!!$ & f_QMixMax & ! (out)
!!$ & )
!!$ f_QMixMax = f_QMixMax * SLTTArcSineFactor + 1.0e-14_DP
!!$ do n = 1, ncmax
!!$ xyzf_QMixA(:,:,:,n) = &
!!$ & 0.5_DP*(asin(2.0_DP*xyzf_QMixA(:,:,:,n)/f_QMixMax(n) - 1.0_DP))
!!$ end do
!!$ end if
!==================================================
! Calculation in a case in which mass fixer applied in horizontal and
! vertical directions in a same time
!
if ( ( .not. FlagSLTTArcsineHor ) .and. ( FlagSLTTArcsineVer ) ) then
! 非負を保証するための arcsine変換フィルタ
! Arcsine transformation for non-negative filter
do n = 1, ncmax
f_QMixProcMax(n) = maxval( xyzf_QMixA(:,:,:,n) )
end do
call MPIWrapperFindMaxVal( &
& ncmax, f_QMixProcMax, & ! (in)
& f_QMixMax & ! (out)
& )
f_QMixMax = f_QMixMax * SLTTArcSineFactor + 1.0e-14_DP
do n = 1, ncmax
xyzf_QMixA(:,:,:,n) = &
& 0.5_DP*(asin(2.0_DP*xyzf_QMixA(:,:,:,n)/f_QMixMax(n) - 1.0_DP))
end do
do n = 1, ncmax
f_QMixLinProcMax(n) = maxval( xyzf_QMixLinA(:,:,:,n) )
end do
call MPIWrapperFindMaxVal( &
& ncmax, f_QMixLinProcMax, & ! (in)
& f_QMixLinMax & ! (out)
& )
f_QMixLinMax = f_QMixLinMax * SLTTArcSineFactor + 1.0e-14_DP
do n = 1, ncmax
xyzf_QMixLinA(:,:,:,n) = &
& 0.5_DP*(asin(2.0_DP*xyzf_QMixLinA(:,:,:,n)/f_QMixLinMax(n) - 1.0_DP))
end do
else if ( ( FlagSLTTArcsineHor ) .and. ( .not. FlagSLTTArcsineVer ) ) then
! 非負を保証するための arcsine変換フィルタ（逆変換）
! Arcsine transformation for non-negative filter
do n = 1, ncmax
xyzf_QMixA(:,:,:,n) = &
& f_QMixMax(n)*(0.5_DP)*(sin(2.0_DP*xyzf_QMixA(:,:,:,n))+1.0_DP)
end do
do n = 1, ncmax
xyzf_QMixLinA(:,:,:,n) = &
& f_QMixLinMax(n)*(0.5_DP)*(sin(2.0_DP*xyzf_QMixLinA(:,:,:,n))+1.0_DP)
end do
end if
xyzf_DQMixDtHorMassFix = 0.0_DP
xyzf_QMixLinATentative = xyzf_QMixLinA
!==================================================
! 鉛直セミラグ
! Vertical
xyzf_QMixA = SLTTVerAdv( xyr_SigDotN, xyzf_QMixA, &
& xyzf_QMixLin = xyzf_QMixLinATentative, & ! (in ) optional
& xyzf_QMixLinA = xyzf_QMixLinA, & ! (out ) optional
& xyzf_QMixMinA = xyzf_QMixMinA, & ! (inout) optional
& xyzf_QMixMaxA = xyzf_QMixMaxA ) ! (inout) optional
! Monotonic filter
! see Diamantakis and Flemming (2014) for BS limiter
! but limiter is applied separately in horizontal and vertical directions
#ifdef SLTT2D1DMONOTONIC
xyzf_QMixA = max( min( xyzf_QMixA, xyzf_QMixMaxA ), xyzf_QMixMinA )
#endif
! Vertical advection by finite difference method
!
!!$ do n = 1, ncmax
!!$ k = 1
!!$ xyrf_QMixA(:,:,k,n) = 1.0e100_DP
!!$ do k = 1, kmax-1
!!$ xyrf_QMixA(:,:,k,n) = &
!!$ & ( xyzf_QMixA(:,:,k,n) + xyzf_QMixA(:,:,k+1,n) ) / 2.0_DP
!!$ end do
!!$ k = kmax
!!$ xyrf_QMixA(:,:,k,n) = 1.0e100_DP
!!$ end do
!!$ do n = 1, ncmax
!!$ do k = 1, kmax
!!$ xyzf_QMixA(:,:,k,n) = xyzf_QMixA(:,:,k,n) &
!!$ & + ( &
!!$ & - ( xyr_SigDotN(:,:,k-1) * xyrf_QMixA(:,:,k-1,n) &
!!$ & - xyr_SigDotN(:,:,k ) * xyrf_QMixA(:,:,k ,n) ) &
!!$ & / z_DelSigma(k) &
!!$ & + xyzf_QMixA(:,:,k,n) &
!!$ & * ( xyr_SigDotN(:,:,k-1) - xyr_SigDotN(:,:,k ) ) &
!!$ & / z_DelSigma(k) &
!!$ & ) * 2.0_DP * DelTime
!!$ end do
!!$ end do
if ( FlagSLTTArcsineVer ) then
! 非負を保証するための arcsine変換フィルタ（逆変換）
! Arcsine transformation for non-negative filter
do n = 1, ncmax
xyzf_QMixA(:,:,:,n) = &
& f_QMixMax(n)*(0.5_DP)*(sin(2.0_DP*xyzf_QMixA(:,:,:,n))+1.0_DP)
end do
do n = 1, ncmax
xyzf_QMixLinA(:,:,:,n) = &
& f_QMixLinMax(n)*(0.5_DP)*(sin(2.0_DP*xyzf_QMixLinA(:,:,:,n))+1.0_DP)
end do
end if
! Mass fixer
!!$ call MassFixerColumn( &
!!$ & xyr_PressA, & ! (in)
!!$ & xyzf_QMixA, & ! (inout)
!!$! & xyr_PressRef = xyr_PressB, & ! (in) optional
!!$ & xyr_PressRef = xyr_PressA, & ! (in) optional
!!$ & xyzf_QMixRef = xyzf_QMixSave & ! (in) optional
!!$ & )
!==================================================
! Calculation in a case in which other type of mass fixer is applied
!
!!$ xyzf_QMixSaveMassFix = xyzf_QMixA
!!$ !
!!$! call MassFixer( &
!!$! call MassFixerWO94( &
!!$ call MassFixerR95( &
!!$ & xyr_PressA, & ! (in)
!!$ & xyzf_QMixA, & ! (inout)
!!$ & xyr_PressRef = xyr_PressB, & ! (in) optional
!!$ & xyzf_QMixRef = xyzf_QMixSave & ! (in) optional
!!$ & )
!!$ !
!!$ xyzf_DQMixDtVerMassFix = 0.0_DP
!==================================================
! Calculation in a case in which mass fixer applied in horizontal and
! vertical directions separately
!
!!$ xyzf_QMixSaveMassFix = xyzf_QMixA
!!$ !
!!$ call MassFixerBC02Column( &
!!$ & xyr_PressA, & ! (in)
!!$ & xyzf_QMixA, & ! (inout)
!!$ & xyzf_QMixLinA, & ! (in)
!!$ & xyr_PressA, & ! (in)
!!$ & xyzf_QMixSave & ! (in)
!!$ & )
!!$ !
!!$ xyzf_DQMixDtVerMassFix = &
!!$ & ( xyzf_QMixA - xyzf_QMixSaveMassFix ) / ( 2.0_DP * DelTime )
!!$ !
!!$ xyzf_DQMixDtTotMassFix = &
!!$ & xyzf_DQMixDtHorMassFix + xyzf_DQMixDtVerMassFix
!==================================================
! Calculation in a case in which mass fixer applied in horizontal and
! vertical directions in a same time
!
xyzf_QMixSaveMassFix = xyzf_QMixA
!
call MassFixerBC02( &
& xyr_PressA, & ! (in)
& xyzf_QMixA, & ! (inout)
& xyzf_QMixLinA, & ! (in)
& xyr_PressB, & ! (in)
& xyzf_QMixSave & ! (in)
& )
!
xyzf_DQMixDtVerMassFix = 0.0_DP
!==================================================
xyzf_DQMixDtTotMassFix = &
& + ( xyzf_QMixA - xyzf_QMixSaveMassFix ) / ( 2.0_DP * DelTime )
do n = 1, ncmax
call HistoryAutoPut( TimeN, 'SLD'//trim(a_QMixName(n))//'DtHorMassFix', &
& xyzf_DQMixDtHorMassFix(:,:,:,n) )
call HistoryAutoPut( TimeN, 'SLD'//trim(a_QMixName(n))//'DtVerMassFix', &
& xyzf_DQMixDtVerMassFix(:,:,:,n) )
call HistoryAutoPut( TimeN, 'SLD'//trim(a_QMixName(n))//'DtTotMassFix', &
& xyzf_DQMixDtTotMassFix(:,:,:,n) )
end do
end subroutine SLTTMain
!----------------------------------------------------------------------------
function SLTTHorAdv( xyzf_QMix, xyz_U, xyz_V, & ! (in)
& xyzf_QMixLinA, & ! (inout) optional
& xyzf_QMixMinA, xyzf_QMixMaxA & ! (out) optional
& ) result( xyzf_QMixA )
! セミラグランジュ法による水平移流の計算
! Calculates tracer transports by Semi-Lagrangian method for horizontal direction
use timeset , only : DelTime, TimeN
! $\Delta t$
use axesset , only : x_Lon, y_Lat
! $\lambda, \varphai$ lon and lat
use sltt_const , only : dtjw, jexmin, jexmax
use sltt_extarr, only : SLTTExtArr, SLTTExtArrf
! 配列拡張ルーチン
! Expansion of arrays
use sltt_dp , only : SLTTDPHor
! 水平上流点探索
! Finding departure point in horizontal
use sltt_lagint, only : SLTTIrrHerIntK13, SLTTIrrLinInt, SLTTLagIntHorMaxMin
! 水平２次元の補間
! 2D Interpolation in horizontal
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
use gtool_history, only: HistoryPut
! SPMODEL ライブラリ, 球面上の問題を球面調和函数変換により解く(多層対応)
! SPMODEL library, problems on sphere are solved with spherical harmonics (multi layer is supported)
!
#ifdef LIB_MPI
!! modify for bug of Intel fortran (2026/01/24 S.Takehiro)
! use ua_mpi_module, only: &
! & wa_xya => ua_pva, &
! & xya_wa => pva_ua, &
! & wa_DLon_wa => ua_DLon_ua, &
! & xya_GradLat_wa => pva_GradLat_ua
use ua_mpi_module_base, only: &
& wa_xya => ua_pva, &
& xya_wa => pva_ua, &
& pva_xvb
use ua_mpi_module_deriv, only: &
& wa_DLon_wa => ua_DLon_ua, &
& xya_GradLat_wa => pva_GradLat_ua
#elif AXISYMMETRY
use wa_zonal_module, only: &
& wa_xya, xya_wa , &
& wa_DLon_wa, xya_GradLat_wa
#else
use wa_module, only: &
& wa_xya, xya_wa , &
& wa_DLon_wa, xya_GradLat_wa
#endif
real(DP), intent(in ) :: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 現在時刻の物質混合比
! Present mix ratio of the tracers
real(DP), intent(in ) :: xyz_U (0:imax-1, 1:jmax, 1:kmax)
! 東西風速
! Zonal Wind
real(DP), intent(in ) :: xyz_V (0:imax-1, 1:jmax, 1:kmax)
! 南北風速
! Meridional Wind
real(DP), intent(inout), optional :: xyzf_QMixLinA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 次ステップの物質混合比
! Next mix ratio of the tracers estimated by linear interpolation
real(DP), intent(out), optional :: xyzf_QMixMinA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP), intent(out), optional :: xyzf_QMixMaxA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP) :: xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 次ステップの物質混合比
! Next mix ratio of the tracers
!
! local variables
!
real(DP) :: xyzf_ExtQMix(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax)
! 現在時刻の物質混合比の拡張配列（南半球）
! Extended array (SH) of present mix ratio of the tracers.
real(DP) :: xyzf_ExtQMixLinA(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax)
! 現在時刻の物質混合比の拡張配列（南半球）
! Extended array (SH) of present mix ratio of the tracers.
real(DP) :: xyz_ExtU (0:imax_global-1, jexmin:jexmax, 1:kc)
! 東西風速の拡張配列（南半球）
! Extended array (SH) of Zonal Wind
real(DP) :: xyz_ExtV (0:imax_global-1, jexmin:jexmax, 1:kc)
! 南北風速の拡張配列（南半球）
! Extended array (SH) of Meridional Wind
integer:: i, ii ! 東西方向に回る 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 dimension of constituents
real(DP) :: xyz_DPLon(0:imax_global-1, 1:jmax, 1:kc)
! 上流点経度（南半球）
! Lon of the departure point (SH)
real(DP) :: xyz_DPLat(0:imax_global-1, 1:jmax, 1:kc)
! 上流点緯度（南半球）
! Lat of the departure point (SH)
!---fx, fy, fxy
real(DP) :: xyzf_QMix_dlon(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 物質混合比の経度微分（グリッド）
! Zonal derivative of the mix ratio (on grid)
real(DP) :: xyzf_QMix_dlat(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 物質混合比の緯度微分（グリッド）
! Meridional derivative of the mix ratio (on grid)
real(DP) :: xyzf_QMix_dlonlat(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 物質混合比の緯度経度微分（グリッド）
! Zonal and meridional derivative of the mix ratio (on grid)
real(DP) :: xyzf_ExtQMix_dlon(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax)
! 物質混合比の経度微分の拡張配列（南半球）
! Extended array (SH) of zonal derivative of the mix ratio
real(DP) :: xyzf_ExtQMix_dlat(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax)
! 物質混合比の緯度微分の拡張配列（南半球）
! Extended array (SH) of meridional derivative of the mix ratio
real(DP) :: xyzf_ExtQMix_dlonlat(0:imax_global-1, jexmin:jexmax, 1:kc, 1:ncmax)
! 物質混合比の緯度経度微分の拡張配列（南半球）
! Extended array (SH) of zonal and meridional derivative of the mix ratio
real(DP) :: wzf_QMix(1:lmax, 1:kmax, 1:ncmax)
! 物質混合比の経度微分（スペクトル）
! Zonal derivative of the mix ratio (on grid)
real(DP) :: wzf_QMix_dlon(1:lmax, 1:kmax, 1:ncmax)
! 物質混合比の経度微分（スペクトル）
! Zonal derivative of the mix ratio (on grid)
real(DP) :: PM ! 配列拡張する際、極ごえ後に符号が変わる場合は -1.0を与える。そうでない場合は1.0を与える。
! Sign change flag for array extension; -1.0 for sign change over the pole, 1.0 for no sign change
!---fxx, fyy, fxxyy
! real(DP) :: xyzf_QMix_dlon2(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_QMix_dlat2(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_QMix_dlon2lat2(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_ExtQMixS_dlon2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_ExtQMixN_dlon2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_ExtQMixS_dlat2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_ExtQMixN_dlat2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_ExtQMixS_dlon2lat2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_ExtQMixN_dlon2lat2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax, 1:ncmax)
!----fxxy
! real(DP) :: xyzf_QMix_dlon2lat(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_ExtQMixS_dlon2lat(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_ExtQMixN_dlon2lat(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax, 1:ncmax)
!----fxyy
! real(DP) :: xyzf_QMix_dlonlat2(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_ExtQMixS_dlonlat2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax, 1:ncmax)
! real(DP) :: xyzf_ExtQMixN_dlonlat2(-2+0:imax-1+3, -jew+1:jmax/2+jew, 1:kmax, 1:ncmax)
!----
! real(DP) :: wzf_QMix_dlon2(1:lmax, 1:kmax, 1:ncmax)
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. sltt_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
! QMixの微分計算（スペクトル変換利用）
! Derivatives of QMix
do n = 1, ncmax
wzf_QMix(:,:,n) = wa_xya(xyzf_QMix(:,:,:,n))
! グリッド→スペクトル
! grid -> spectrum
xyzf_QMix_dlat(:,:,:,n) = xya_GradLat_wa(wzf_QMix(:,:,n))
! スペクトル→グリッド緯度微分
! spectrum -> grid (dQ/dlat)
wzf_QMix_dlon(:,:,n) = wa_Dlon_wa(wzf_QMix(:,:,n))
! スペクトル→スペクトル経度微分
! spectrum -> spectrum (dQ/dlon)
xyzf_QMix_dlon(:,:,:,n) = xya_wa(wzf_QMix_dlon(:,:,n))
! スペクトル経度微分→グリッド経度微分
! spectrum (dQ/dlon) -> grid (dQ/dlon)
xyzf_QMix_dlonlat(:,:,:,n) = xya_GradLat_wa(wzf_QMix_dlon(:,:,n))
! スペクトル経度微分→グリッド緯度経度微分
! spectrum (dQ/dlon) -> grid (d^2Q/dlon dlat)
!---fxx, fyy, fxxy, fxyy, fxxyy を計算
!xyzf_QMix_dlon2(:,:,:,n) = xya_wa(wa_Dlon_wa(wzf_QMix_dlon(:,:,n)))
!xyzf_QMix_dlat2(:,:,:,n) = xya_GradLat_wa(wa_xya(xyzf_QMix_dlat(:,:,:,n)))
!xyzf_QMix_dlon2lat(:,:,:,n) = xya_GradLat_wa(wa_xya(xyzf_QMix_dlon2(:,:,:,n)))
!xyzf_QMix_dlonlat2(:,:,:,n) = xya_GradLat_wa(wa_xya(xyzf_QMix_dlonlat(:,:,:,n)))
!xyzf_QMix_dlon2lat2(:,:,:,n) = xya_GradLat_wa(wa_xya(xyzf_QMix_dlon2lat(:,:,:,n)))
enddo
! 配列の分割と拡張
! Division and extension of arrays
!
! 配列の分割と拡張
! Division and extension of arrays
! 配列拡張する際、極ごえ後に符号が変わる場合は -1.0を与える。そうでない場合は1.0を与える。
!!$ pm = 1.0_DP
pm = -1.0_DP
call SLTTExtArrf( &
& xyzf_QMix_dlon, PM, & ! (in)
& xyzf_ExtQMix_dlon & ! (out)
& )
pm = -1.0_DP
call SLTTExtArrf( &
& xyzf_QMix_dlat, PM, & ! (in)
& xyzf_ExtQMix_dlat & ! (out)
& )
pm = -1.0_DP
call SLTTExtArrf( &
& xyzf_QMix_dlonlat, PM, & ! (in)
& xyzf_ExtQMix_dlonlat & ! (out)
& )
PM = 1.0_DP
call SLTTExtArrf( &
& xyzf_QMix, PM, & ! (in)
& xyzf_ExtQMix & ! (out)
& )
#if defined(AXISYMMETRY)
xyz_ExtU = 0.0D0
#else
PM = -1.0_DP
call SLTTExtArr( &
& xyz_U, PM, & ! (in)
& xyz_ExtU & ! (out)
& )
#endif
PM = -1.0_DP
call SLTTExtArr( &
& xyz_V, PM, & ! (in)
& xyz_ExtV & ! (out)
& )
if ( present( xyzf_QMixLinA ) ) then
! Extention of array for linear interpolation
PM = 1.0_DP
! call SLTTExtArrExt2( &
! & x_SinLon, x_CosLon, & ! (in)
! & xyzf_QMixLinA, PM, & ! (in)
! & xyzf_ExtQMixLinA & ! (out)
! & )
call SLTTExtArrf( &
& xyzf_QMixLinA, PM, & ! (in)
& xyzf_ExtQMixLinA & ! (out)
& )
end if
if ( sltt_check_inited ) then
call HistoryPut('Q',xyzf_ExtQMix(:,:,:,1))
call HistoryPut('Qlon',xyzf_ExtQMix_dlon(:,:,:,1))
call HistoryPut('Qlat',xyzf_ExtQMix_dlat(:,:,:,1))
call HistoryPut('Qlonlat',xyzf_ExtQMix_dlonlat(:,:,:,1))
call HistoryPut('U',xyz_ExtU)
call HistoryPut('V',xyz_ExtV)
end if
! 上流点の計算
! estimation of departure point
call SLTTDPHor( &
& DelTime, x_Lon, y_Lat, y_SinLat, y_CosLat, & ! (in)
& jexmin, jexmax, & ! (in)
& x_ExtLon, y_ExtLat, xyz_ExtU, xyz_ExtV, & ! (in)
& xyz_DPLon, xyz_DPLat & ! (out)
& )
call HistoryAutoPut( TimeN, 'Qlat', xyzf_QMix_dlat(:,:,:,1) )
#ifdef LIB_MPI
call HistoryAutoPut( TimeN, 'DPLon', pva_xvb(xyz_DPLon) )
call HistoryAutoPut( TimeN, 'DPLat', pva_xvb(xyz_DPLat) )
#else
call HistoryAutoPut( TimeN, 'DPLon', xyz_DPLon )
call HistoryAutoPut( TimeN, 'DPLat', xyz_DPLat )
#endif
!=========================== ここまで ====================================
! 補間
! Interpolation
! do n = 1, ncmax
call SLTTIrrHerIntK13( &
& jexmin, jexmax, & ! (in)
& x_ExtLon, y_ExtLat, xyz_DPLon, xyz_DPLat, & ! (in)
& xyzf_ExtQMix(:,:,:,:), xyzf_ExtQMix_dlon(:,:,:,:), & ! (in)
& xyzf_ExtQMix_dlat(:,:,:,:), xyzf_ExtQMix_dlonlat(:,:,:,:), & ! (in)
& SLTTIntHor, & ! (in)
& xyzf_QMixA(:,:,:,:) & ! (out)
& )
! enddo
if ( present( xyzf_QMixLinA ) ) then
call SLTTIrrLinInt( &
& jexmin, jexmax, & ! (in)
& x_ExtLon, y_ExtLat, xyz_DPLon, xyz_DPLat, xyzf_ExtQMixLinA, & ! (in)
& xyzf_QMixLinA & ! (out)
& )
end if
if ( ( ( present( xyzf_QMixMinA ) ) .and. &
& ( .not. present( xyzf_QMixMaxA ) ) ) .or. &
& ( ( .not. present( xyzf_QMixMinA ) ) .and. &
& ( present( xyzf_QMixMaxA ) ) ) ) then
call MessageNotify( 'E', module_name, &
& 'QMixMinA has to be present when QMixMaxA is present, and vice versa.' )
end if
if ( present( xyzf_QMixMinA ) ) then
call SLTTLagIntHorMaxMin( &
& jexmin, jexmax, & ! (in)
& x_ExtLon, y_ExtLat, xyz_DPLon, xyz_DPLat, xyzf_ExtQMix, & ! (in)
& xyzf_QMixMinA, xyzf_QMixMaxA & ! (out)
& )
end if
end function SLTTHorAdv
!--------------------------------------------------------------------------------------
function SLTTVerAdv( xyr_SigmaDot, xyzf_QMix, &
& xyzf_QMixLin, & ! (in ) optional
& xyzf_QMixLinA, & ! (out) optional
& xyzf_QMixMinA, xyzf_QMixMaxA & ! (out) optional
& ) result( xyzf_QMixA )
! セミラグランジュ法による鉛直移流の計算
! Calculates tracer transports by Semi-Lagrangian method for vertical direction
use axesset, only : z_Sigma ! 鉛直座標; Sigma coordinate
use timeset, only : DelTime ! $\Delta t$
use sltt_dp, only : SLTTDPVer ! 鉛直上流点探索; Finding departure point in vertical
use sltt_lagint, only : &
& SLTTIrrHerIntQui1DNonUni, & ! 不当間隔格子の五次補間; Quintic Interpolation for non-uniform grids
& SLTTHerIntCub1D
real(DP), intent(in ) :: xyr_SigmaDot(0:imax-1, 1:jmax, 0:kmax)
! 鉛直流速（SigmaDot）
real(DP), intent(in ) :: xyzf_QMix (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 現在時刻の物質混合比
! Present mix ratio of the tracers
real(DP), intent(in ), optional :: xyzf_QMixLin (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP), intent(out), optional :: xyzf_QMixLinA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP), intent(out), optional :: xyzf_QMixMinA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP), intent(out), optional :: xyzf_QMixMaxA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP) :: xyzf_QMixA (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 次ステップの物質混合比
! Next mix ratio of the tracers
!
! local variables
!
real(DP) :: xyz_DPSigma(0:imax-1, 1:jmax, 1:kmax)
! 上流点高度
! Sigma of the departure point
integer:: i ! 東西方向に回る DO ループ用作業変数
! Work variables for DO loop in zonal direction
integer:: j ! 南北方向に回る DO ループ用作業変数
! Work variables for DO loop in meridional direction
integer:: k, kk ! 鉛直方向に回る DO ループ用作業変数
! Work variables for DO loop in vertical direction
integer:: n ! 組成方向に回る DO ループ用作業変数
! Work variables for DO loop in dimension of constituents
integer:: xy_kk(0:imax-1, 1:jmax)
! 上流点の上下のグリッドを探索するための作業変数
! Work variable for finding the grid just above the departure point
real(DP) :: xyzf_QMix_dz(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
! 物質混合比の鉛直微分
! Vertical derivative of the mix ratio
real(DP) :: xyzf_ExtQMix(0:imax-1, 1:jmax, 1-2:kmax+2, 1:ncmax)
! 物質混合比の拡張配列
! Extended array of the mix ratio
real(DP) :: z_ExtSigma(1-2:kmax+2)
! σ座標の拡張配列
! Extended array of the sigma coordinate
real(DP) :: xyf_F11(0:imax-1, 1:jmax, 1:ncmax)
! 微分計算時に用いる作業変数
! work variable for the derivative calculation
real(DP) :: xyf_F22(0:imax-1, 1:jmax, 1:ncmax)
! 微分計算時に用いる作業変数
! work variable for the derivative calculation
real(DP) :: xyf_F12(0:imax-1, 1:jmax, 1:ncmax)
! 微分計算時に用いる作業変数
! work variable for the derivative calculation
real(DP) :: xyf_F21(0:imax-1, 1:jmax, 1:ncmax)
! 微分計算時に用いる作業変数
! work variable for the derivative calculation
real(DP) :: s1, t1, s2, t2, r1, r2
! 微分計算時に用いる作業変数
! work variable for the derivative calculation
real(DP) :: xyzf_QMixLinLV (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
real(DP) :: xyzf_ExtQMixLinLV(0:imax-1, 1:jmax, 1-2:kmax+2, 1:ncmax)
! 実行文 ; Executable statement
!
! 初期化確認
! Initialization check
!
if ( .not. sltt_inited ) then
call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
end if
if ( ( present( xyzf_QMixLin ) ) .and. ( .not. present( xyzf_QMixLinA ) ) ) then
call MessageNotify( 'E', module_name, &
& 'If xyzf_QMixLinA has to be present when xyzf_QMixLin ise present.' )
end if
if ( ( ( present( xyzf_QMixMinA ) ) .and. &
& ( .not. present( xyzf_QMixMaxA ) ) ) .or. &
& ( ( .not. present( xyzf_QMixMinA ) ) .and. &
& ( present( xyzf_QMixMaxA ) ) ) ) then
call MessageNotify( 'E', module_name, &
& 'QMixMinA has to be present when QMixMaxA is present, and vice versa.' )
end if
if ( present( xyzf_QMixLin ) ) then
xyzf_QMixLinLV = xyzf_QMixLin
else
xyzf_QMixLinLV = xyzf_QMix
end if
! 上流点探索
! estimation of departure point
!
call SLTTDPVer( &
& DelTime, xyr_SigmaDot, & ! (in )
& xyz_DPSigma & ! (out)
& )
! 配列拡張（z_Sigma）
! Array extension for z_Sigma
z_ExtSigma(-1) = 2.0_DP - z_Sigma(2)
z_ExtSigma(0) = 2.0_DP - z_Sigma(1)
z_ExtSigma(1:kmax) = z_Sigma(1:kmax)
z_ExtSigma(kmax+1) = -z_Sigma(kmax)
z_ExtSigma(kmax+2) = -z_Sigma(kmax-1)
! 配列拡張（xyzf_QMix）
! Array extension for Q_Mix
xyzf_ExtQMix(:,:,-1,:) = xyzf_QMix(:,:,2,:)
xyzf_ExtQMix(:,:,0,:) = xyzf_QMix(:,:,1,:)
xyzf_ExtQMix(:,:,1:kmax,:) = xyzf_QMix(:,:,1:kmax,:)
xyzf_ExtQMix(:,:,kmax+1,:) = xyzf_QMix(:,:,kmax,:)
xyzf_ExtQMix(:,:,kmax+2,:) = xyzf_QMix(:,:,kmax-1,:)
xyzf_ExtQMixLinLV(:,:,-1,:) = xyzf_QMixLinLV(:,:,2,:)
xyzf_ExtQMixLinLV(:,:,0,:) = xyzf_QMixLinLV(:,:,1,:)
xyzf_ExtQMixLinLV(:,:,1:kmax,:) = xyzf_QMixLinLV(:,:,1:kmax,:)
xyzf_ExtQMixLinLV(:,:,kmax+1,:) = xyzf_QMixLinLV(:,:,kmax,:)
xyzf_ExtQMixLinLV(:,:,kmax+2,:) = xyzf_QMixLinLV(:,:,kmax-1,:)
! xyzf_QMix_dz（微分）を求める
! calculate xyzf_QMix_dz
do k = 1 , kmax
s1 = z_ExtSigma(k) - z_ExtSigma(k-1)
t1 = z_ExtSigma(k+1) - z_ExtSigma(k)
s2 = z_ExtSigma(k) - z_ExtSigma(k-2)
t2 = z_ExtSigma(k+2) - z_ExtSigma(k)
if (s1 == t1 .and. s2 == t2 .and. s1 + s1 == s2) then
! 格子が等間隔の場合
! Uniform depth
! 4次精度
! 4th order
xyzf_QMix_dz(:,:,k,:) = ( 8.0_DP*( xyzf_ExtQMix(:,:,k+1,:) - xyzf_ExtQMix(:,:,k-1,:)) &
& - ( xyzf_ExtQMix(:,:,k+2,:) - xyzf_ExtQMix(:,:,k-2,:) ) )/12.0_DP
else
! 格子が不当間隔の場合
! Non-uniform depth
xyf_F11 = (s1*s1*xyzf_ExtQMix(:,:,k+1,:) +(t1*t1 - s1*s1)*xyzf_ExtQMix(:,:,k,:) - t1*t1*xyzf_ExtQMix(:,:,k-1,:))&
& /(s1*t1*(s1+t1))
xyf_F22 = (s2*s2*xyzf_ExtQMix(:,:,k+2,:) +(t2*t2 - s2*s2)*xyzf_ExtQMix(:,:,k,:) - t2*t2*xyzf_ExtQMix(:,:,k-2,:))&
& /(s2*t2*(s2+t2))
xyf_F21 = (s2*s2*xyzf_ExtQMix(:,:,k+1,:) +(t1*t1 - s2*s2)*xyzf_ExtQMix(:,:,k,:) - t1*t1*xyzf_ExtQMix(:,:,k-2,:))&
& /(s2*t1*(s2+t1))
xyf_F12 = (s1*s1*xyzf_ExtQMix(:,:,k+2,:) +(t2*t2 - s1*s1)*xyzf_ExtQMix(:,:,k,:) - t2*t2*xyzf_ExtQMix(:,:,k-1,:))&
& /(s1*t2*(s1+t2))
r1 = t1 - s1 - t2 + s2
r2 = t1 - s2 - t2 + s1
!４次精度
! 4th order
xyzf_QMix_dz(:,:,k,:) = ( (xyf_F11*s2*t2 - xyf_F22*s1*t1)*r2 - (xyf_F21*s1*t2 - xyf_F12*s2*t1)*r1 ) &
& / ( (s2*t2-s1*t1)*r2 - (s1*t2-s2*t1)*r1 )
!3次精度
! 3rd order
! xyzf_QMix_dz(:,:,k,:) = (xyf_F11*s2*t2 - xyf_F22(:,:,:)*s1*t1)/(s2*t2 - s1*t1)
!2次精度
! 2nd order
! xyzf_QMix_dz(:,:,k,:) = xyf_F11
end if
end do
xy_kk = 2
do k = 1, kmax
do j = 1, jmax
do i = 0, imax-1
if ( xyz_DPSigma(i,j,k) >= z_Sigma(1) ) then ! DPが z_Sigma(1) と 地表面(sigma = 1.0)の間の場合
! if DP is between z_Sigma(1) and the ground (sigma = 1.0)
xyzf_QMixA(i,j,k,:) = xyzf_QMix(i,j,1,:) ! Q_1で一定とする。
! use Q_1 for interpolated value
if ( present( xyzf_QMixLinA ) ) then
xyzf_QMixLinA(i,j,k,:) = xyzf_QMixLinLV(i,j,1,:)
end if
if ( present( xyzf_QMixMinA ) ) then
xyzf_QMixMinA(i,j,k,:) = xyzf_QMix(i,j,1,:)
xyzf_QMixMaxA(i,j,k,:) = xyzf_QMix(i,j,1,:)
end if
elseif (xyz_DPSigma(i,j,k) <= z_Sigma(kmax)) then! DPが z_Sigma(kmax) と 大気上端(sigma = 0.0)の間
! if DP is between z_Sigma(kmax) and the upper boundary (sigma = 0.0)
xyzf_QMixA(i,j,k,:) = xyzf_QMix(i,j,kmax,:) ! Q_kmaxで一定とする。
! use Q_kmax for interpolated value
if ( present( xyzf_QMixLinA ) ) then
xyzf_QMixLinA(i,j,k,:) = xyzf_QMixLinLV(i,j,kmax,:)
end if
if ( present( xyzf_QMixMinA ) ) then
xyzf_QMixMinA(i,j,k,:) = xyzf_QMix(i,j,kmax,:)
xyzf_QMixMaxA(i,j,k,:) = xyzf_QMix(i,j,kmax,:)
end if
else
do kk = xy_kk(i,j), kmax
if ( xyz_DPSigma(i,j,k) > z_Sigma(kk) ) then
select case (SLTTIntVer)
case("HQ") ! 変則エルミート５次補間; Irregular Hermite Quintic interpolation
do n = 1, ncmax
xyzf_QMixA(i,j,k,n) = SLTTIrrHerIntQui1DNonUni(xyzf_ExtQMix(i,j,kk-2,n), xyzf_ExtQMix(i,j,kk-1,n), &
& xyzf_ExtQMix(i,j,kk,n), xyzf_ExtQMix(i,j,kk+1,n), &
& xyzf_QMix_dz(i,j,kk-1,n), xyzf_QMix_dz(i,j,kk,n), &
& z_ExtSigma(kk-2)-z_ExtSigma(kk-1), z_ExtSigma(kk)-z_ExtSigma(kk-1), &
& z_ExtSigma(kk+1)-z_ExtSigma(kk-1), xyz_DPSigma(i,j,k)-z_ExtSigma(kk-1))
end do
case("HC") ! エルミート３次補間; Hermitian Cubic interpolation
do n = 1, ncmax
xyzf_QMixA(i,j,k,n) = SLTTHerIntCub1D( xyzf_ExtQMix(i,j,kk-1,n), xyzf_ExtQMix(i,j,kk,n),&
& xyzf_QMix_dz(i,j,kk-1,n), xyzf_QMix_dz(i,j,kk,n),&
& z_ExtSigma(kk)-z_ExtSigma(kk-1), &
& xyz_DPSigma(i,j,k)-z_ExtSigma(kk-1))
end do
case default
call MessageNotify( 'E', module_name, 'GIVE CORRECT KEYWORD FOR IN NAMELIST.' )
end select
if ( present( xyzf_QMixLinA ) ) then
! Linear interporation
do n = 1, ncmax
xyzf_QMixLinA(i,j,k,n) = &
& ( xyzf_ExtQMixLinLV(i,j,kk,n) - xyzf_ExtQMixLinLV(i,j,kk-1,n) ) &
& / ( z_ExtSigma(kk) - z_ExtSigma(kk-1) ) &
& * ( xyz_DPSigma(i,j,k) - z_ExtSigma(kk-1) ) &
& + xyzf_ExtQMixLinLV(i,j,kk-1,n)
end do
end if
if ( present( xyzf_QMixMinA ) ) then
do n = 1, ncmax
xyzf_QMixMinA(i,j,k,n) = &
& min( xyzf_QMix(i,j,kk-1,n), &
& xyzf_QMix(i,j,kk,n) )
xyzf_QMixMaxA(i,j,k,n) = &
& max( xyzf_QMix(i,j,kk-1,n), &
& xyzf_QMix(i,j,kk,n) )
end do
end if
xy_kk(i,j) = kk
exit
end if
end do
end if
end do
end do
end do
end function SLTTVerAdv
!-------------------------------------------------
subroutine SLTTInit
! セミラグランジュ法の初期化処理
! Initialization for Semi-Lagrangian method
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoAddVariable
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: &
& ncmax, &
! 成分の数
! Number of composition
& a_QMixName
! 成分の変数名
! Name of variables for composition
! 座標データ設定
! Axes data settings
!
use axesset, only: &
& r_Sigma, &
! $ \sigma $ レベル (半整数).
! Half $ \sigma $ level
& z_Sigma, & ! $ \sigma $ レベル (整数).
! Full $ \sigma $ level
& x_Lon, y_Lat, &
& AxNameX, AxNameY, AxNameZ, AxNameT
use sltt_const , only : SLTTConstInit
use sltt_extarr, only : SLTTExtArrInit
! NAMELIST ファイル入力に関するユーティリティ
! Utilities for NAMELIST file input
!
use namelist_util, only: namelist_filename, NmlutilMsg
! メッセージ制御
! Message control
!
use mpi_messagecntl, only : DoesOutputMPIMessage
! 種別型パラメタ
! Kind type parameter
!
use dc_types, only: &
& STDOUT ! 標準出力の装置番号. Unit number of standard output
! ファイル入出力補助
! File I/O support
!
use dc_iounit, only: FileOpen
use sltt_const , only : jexmin, jexmax
! メッセージ制御
! Message control
!
use mpi_messagecntl, only : DoesOutputMPIMessage
!
! local variables
!
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
integer:: unit_nml ! NAMELIST ファイルオープン用装置番号.
! Unit number for NAMELIST file open
integer:: iostat_nml ! NAMELIST 読み込み時の IOSTAT.
! IOSTAT of NAMELIST read
! NAMELIST 変数群
! NAMELIST group name
!
namelist /sltt_nml/ &
& FlagSLTTArcsineHor, FlagSLTTArcsineVer, SLTTIntHor, SLTTIntVer, SLTTArcSineFactor
! 実行文 ; Executable statement
!
if ( sltt_inited ) return
if ( mod( jmax, 2 ) /= 0 ) then
stop 'jmax cannot be divided by 2.'
end if
call SLTTConstInit
! デフォルト値の設定
! Default values settings
!
FlagSLTTArcsineHor = .true.
FlagSLTTArcsineVer = .true.
SLTTArcSineFactor = 1.05_DP
SLTTIntHor = "HQ"
SLTTIntVer = "HQ"
! NAMELIST の読み込み
! NAMELIST is input
!
if ( trim(namelist_filename) /= '' ) then
call FileOpen( unit_nml, & ! (out)
& namelist_filename, mode = 'r' ) ! (in)
rewind( unit_nml )
read( unit_nml, & ! (in)
& nml = sltt_nml, & ! (out)
& iostat = iostat_nml ) ! (out)
close( unit_nml )
call NmlutilMsg( iostat_nml, module_name ) ! (in)
if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = sltt_nml )
end if
allocate( x_SinLon(0:imax-1) )
allocate( x_CosLon(0:imax-1) )
allocate( y_SinLat(1:jmax) )
allocate( y_CosLat(1:jmax) )
do i = 0, imax-1
x_SinLon(i) = sin( x_Lon(i) )
x_CosLon(i) = cos( x_Lon(i) )
end do
do j = 1, jmax
y_SinLat(j) = sin( y_Lat(j) )
y_CosLat(j) = cos( y_Lat(j) )
end do
allocate( x_ExtLon( 0:imax_global-1 ) )
allocate( y_ExtLat( jexmin:jexmax ) )
call SLTTExtArrInit( &
& x_Lon, y_Lat, & ! (in )
& x_ExtLon, y_ExtLat & ! (out)
& )
! ヒストリデータ出力のためのへの変数登録
! Register of variables for history data output
!
do n = 1, ncmax
call HistoryAutoAddVariable( 'SLD'//trim(a_QMixName(n))//'DtHorMassFix', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'tendency of horizontal mass fix of '//trim(a_QMixName(n)), 's-1' )
call HistoryAutoAddVariable( 'SLD'//trim(a_QMixName(n))//'DtVerMassFix', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'tendency of vertical mass fix of '//trim(a_QMixName(n)), 's-1' )
call HistoryAutoAddVariable( 'SLD'//trim(a_QMixName(n))//'DtTotMassFix', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'tendency of mass fix of '//trim(a_QMixName(n)), 's-1' )
end do
call HistoryAutoAddVariable( 'Qlat', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'Qlat', '1' )
call HistoryAutoAddVariable( 'DPLon', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'DPLon', '1' )
call HistoryAutoAddVariable( 'DPLat', &
& (/ AxNameX, AxNameY, AxNameZ, AxNameT /), &
& 'DPLat', '1' )
! 印字 ; Print
!
call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
call MessageNotify( 'M', module_name, ' FlagSLTTArcsineHor = %b', l = (/ FlagSLTTArcsineHor /) )
call MessageNotify( 'M', module_name, ' FlagSLTTArcsineVer = %b', l = (/ FlagSLTTArcsineVer /) )
call MessageNotify( 'M', module_name, ' SLTTArcsineFactor = %f', d = (/ SLTTArcsineFactor /) )
call MessageNotify( 'M', module_name, ' SLTTIntHor = %c', c1 = trim( SLTTIntHor ) )
call MessageNotify( 'M', module_name, ' SLTTIntVer = %c', c1 = trim( SLTTIntVer ) )
call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )
sltt_inited = .true.
end subroutine SLTTInit
!
!----------------------------------------------------------------------
!
subroutine SLTTCheckInit
use dc_string, only : cprintf
use dc_types, only: &
& STRING ! 文字列. Strings.
use sltt_const , only : jexmin, jexmax, PIx2
use gtool_history, only: HistoryCreate, HistoryAddVariable, HistoryPut
use axesset, only: &
& z_Sigma ! $ \sigma $ (整数).
! $ \sigma $ (Full)
character(STRING) :: testfilename
!
! チェック用出力
!
testfilename = cprintf("sltt_test_%04d.nc",(/myrank/))
call HistoryCreate( & ! ヒストリー作成
& file=trim(adjustl(testfilename)), title='sltt test', &
& source='Sample program of gtool_history/gtool5', &
& institution='GFD_Dennou Club davis project', &
& dims=(/'lon ','lat ','sig ','time'/), &
& dimsizes=(/imax_global,jexmax-jexmin+1,kc,0/),&
& longnames=(/'longitude','latitude ','sigma ','time '/),&
& units=(/'deg','deg','1 ','s '/), &
& origin=0.0, interval=real(DelTime) )
call HistoryPut('lon',x_ExtLon*360.0D0/PIx2) ! 次元変数出力
call HistoryPut('lat',y_ExtLat*360.0D0/PIx2) ! 次元変数出力
! call HistoryPut('sig',z_Sigma) ! 次元変数出力
call HistoryAddVariable( & ! 変数定義
& varname='Q', dims=(/'lon ','lat ','sig ','time'/), &
& longname='QMix', units='1', xtype='double')
call HistoryAddVariable( & ! 変数定義
& varname='Qlon', dims=(/'lon ','lat ','sig ','time'/), &
& longname='QMix (lon.derivative)', units='1', xtype='double')
call HistoryAddVariable( & ! 変数定義
& varname='Qlat', dims=(/'lon ','lat ','sig ','time'/), &
& longname='QMix (lat.derivative)', units='1', xtype='double')
call HistoryAddVariable( & ! 変数定義
& varname='Qlonlat', dims=(/'lon ','lat ','sig ','time'/), &
& longname='QMix (lon.lat.derivative)', units='1', xtype='double')
call HistoryAddVariable( & ! 変数定義
& varname='U', dims=(/'lon ','lat ','sig ','time'/), &
& longname='U', units='m/s', xtype='double')
call HistoryAddVariable( & ! 変数定義
& varname='V', dims=(/'lon ','lat ','sig ','time'/), &
& longname='V', units='m/s', xtype='double')
sltt_check_inited = .true.
end subroutine SLTTCheckInit
!--------------------------------------------------------------------------------------
end module sltt