subroutine Scalar( ss_Scalar_in, DelTime, fs_VelX_adv, sf_VelZ_adv, ss_Scalar_adv, ss_Scalar_dif, ss_Kh_dif, ss_Scalar_out)
!=begin
!=end
!== 暗黙の型宣言禁止
implicit none
!=begin
!== Input
real(8), intent(in) :: DelTime
real(8), intent(in) :: ss_Scalar_in(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(in) :: fs_VelX_adv(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(in) :: sf_VelZ_adv(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(in) :: ss_Scalar_adv(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(in) :: ss_Scalar_dif(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(in) :: ss_Kh_dif(DimXMin:DimXMax, DimZMin:DimZMax)
!== Output
real(8), intent(out) :: ss_Scalar_out(DimXMin:DimXMax, DimZMin:DimZMax)
!=end
!== Work
real(8) :: ss_AdvScalar(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: ss_SrcScalar(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: ss_TurbScalar(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: ss_NumDiffScalar(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: ss_TendScalar(DimXMin:DimXMax, DimZMin:DimZMax)
call BeginSub("Scalar", fmt="%c", c1="Time integration of scalar component.")
!=== 移流項の計算. 中心差分を利用して計算.
ss_AdvScalar = ss_avr_fs( fs_VelX_adv * fs_dx_ss( ss_Scalar_adv ) ) + ss_avr_sf( sf_VelZ_adv * sf_dz_ss( ss_Scalar_adv ) )
! !=== 移流項の計算. MPDATA スキームを利用して計算
! call AdvectScalar_MPDATA( fs_VelX_adv, sf_VelZ_adv,
! ss_Scalar_adv, DelTime,
! ss_AdvScalar_adv )
!=== 生成項
ss_SrcScalar = ss_avr_sf( sf_VelZ_adv * sf_dz_ss( ss_PotTempBasicZ ) )
!=== 乱流拡散項 (時刻は t - Δt で評価)
ss_TurbScalar = ss_dx_fs( fs_avr_ss(ss_Kh_dif) * fs_dx_ss( ss_Scalar_dif ) ) + ss_dz_sf( sf_avr_ss(ss_Kh_dif) * sf_dz_ss( ss_Scalar_dif ) )
!=== 数値粘性項 (時刻は t - Δt で評価)
ss_NumDiffScalar = NuH * ( ss_dx_fs( fs_dx_ss( ss_Scalar_dif) ) ) + NuV * ( ss_dz_sf( sf_dz_ss( ss_Scalar_dif) ) )
!=== 項をまとめる
ss_TendScalar = - ss_AdvScalar - ss_SrcScalar + ss_TurbScalar + ss_NumDiffScalar
!=== 温位の時間積分
ss_Scalar_out = ss_Scalar_in + DelTime * ss_TendScalar
!=== 境界条件
call boundary(ss_BC, ss_Scalar_out)
call EndSub("Scalar")
contains
!=begin
!== Subroutine AdvectScalar_MPDATA
!
! * Developer: ODAKA Masatsugu (odakker@gfd-dennou.org)
! * Version: $Id: scalar.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $
! * Tag Name: $Name: $
! * Change History:
!
!=== Overview
!
! スカラー量の移流を MPDATA スキームを用いて計算する.
!
!=== Error Handling
!
!=== Known Bugs
!
!=== Note
!
! スキームの都合上, 時間積分は前進差分を用いて行う.
! 実引数に渡す変数配列の時間レベルに注意すること.
!
! 非定常な流れの場合, 引数 fs_VelX, fs_VelZ の時間レベルは ss_Scalar の
! 時間レベルから DelTime/2 だけ先行させたものにする.
!
! 速度場を leapfrog スキームで解く場合に簡単に実装するなら, 時刻 t における
! fs_VelX_n, fs_VelZ_n を用いて時間刻み 2*DelTime の前進差分として ss_Scala
! を時間積分する.
!
! t-Δt t t+Δt
! -----|-------------|-------------|-----
! [fs_VelX_n,sf_VelZ_n]
! ss_Scalar_b ss_Scalar_a
! |======== 2*DelTime =======>|
!
!
! 精度良く計算するなら時刻 t+Δt/2 における fs_VelX_mid, sf_VelZ_mid を
! 用いて時間刻み DelTime の前進差分として ss_Scalar を計算する.
!
! t-Δt t t+Δt
! -----|-------------|------|------|-----
! [fs_VelX_mid,sf_VelZ_mid]
! s_Scalar_n ss_Scalar_a
! |== DelTime =>|
!
!=== Future Plans
!
!=end
subroutine AdvectScalar_MPDATA(fs_VelX, sf_VelZ, ss_Scalar, DelTime, ss_AdvScalar)
!=begin
!=== Dependency
!=end
!=== 暗黙の型宣言禁止
implicit none
!=begin
!=== Input
real(8), intent(in) :: fs_VelX(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(in) :: sf_VelZ(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(in) :: ss_Scalar(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(in) :: DelTime
!=== Output
real(8), intent(out) :: ss_AdvScalar(DimXMin:DimXMax, DimZMin:DimZMax)
!=end
!=== Work
real(8) :: ss_Scalar_tmp(DimXMin:DimXMax, DimZMin:DimZMax) ! 仮の値
real(8) :: ss_AdvScalar_tmp(DimXMin:DimXMax, DimZMin:DimZMax)! 仮の値
real(8) :: ss_DivVel(DimXMin:DimXMax, DimZMin:DimZMax) ! 発散
real(8) :: fs_FluxX(DimXMin:DimXMax, DimZMin:DimZMax) ! フラックス
real(8) :: sf_FluxZ(DimXMin:DimXMax, DimZMin:DimZMax) ! フラックス
real(8) :: fs_AntiDiffVelX(DimXMin:DimXMax, DimZMin:DimZMax) ! 反拡散速度
real(8) :: sf_AntiDiffVelZ(DimXMin:DimXMax, DimZMin:DimZMax) ! 反拡散速度
real(8) :: fs_AntiDiffFluxX(DimXMin:DimXMax, DimZMin:DimZMax)! フラックス
real(8) :: sf_AntiDiffFluxZ(DimXMin:DimXMax, DimZMin:DimZMax)! フラックス
call BeginSub("AdvectScalar_MPDATA", fmt="%c", c1="Calculate advection term for scalar variables (ss_AdvScalar).")
!=== 初期化
fs_FluxX = 0.0d0
sf_FluxZ = 0.0d0
ss_AdvScalar = 0.0d0
!=== 上流差分法を用いてスカラーのフラックスを計算
call FluxScalar(fs_VelX, sf_VelZ, ss_Scalar, fs_FluxX, sf_FluxZ)
!=== 発散の計算
ss_DivVel = ss_dx_fs( fs_VelX ) + ss_dz_sf( sf_VelZ )
!=== 移流による増分を計算
ss_AdvScalar = ss_dx_fs(fs_FluxX) + ss_dz_sf(sf_FluxZ)
! - ss_Scalar * ss_DivVel
!=== スカラーの仮値を計算
ss_Scalar_tmp = ss_Scalar + DelTime * ( - ss_AdvScalar )
! call Integral(- ss_AdvScalar, DelTime, ss_Scalar, ss_Scalar_tmp )
call boundary(ss_BC, ss_Scalar_tmp)
!=== MPDATA スキームによる補正
!=== 初期化
fs_AntiDiffVelX = 0.0d0
sf_AntiDiffVelZ = 0.0d0
fs_AntiDiffFluxX = 0.0d0
sf_AntiDiffFluxZ = 0.0d0
ss_AdvScalar_tmp = 0.0d0
ss_DivVel = 0.0d0
!=== 上流差分によるスカラーフラックスが 0 でない場所だけ反拡散速度を計算
where(fs_FluxX /= 0.0d0)
fs_AntiDiffVelX = 0.5d0 * ( (abs(fs_VelX)*DelX - DelTime*fs_VelX**2)*fs_dx_ss(ss_Scalar_tmp) - DelTime*fs_VelX*fs_avr_sf(sf_VelZ)*fs_avr_sf(sf_dz_ss(ss_Scalar_tmp))) / (Max(fs_avr_ss(ss_Scalar_tmp), epsilon(0.0d0))) - 0.5d0*fs_VelX*fs_avr_ss(ss_DivVel)
end where
where (sf_FluxZ /= 0.0d0)
sf_AntiDiffVelZ = 0.5d0 * ( (abs(sf_VelZ)*DelZ - DelTime*sf_VelZ**2)*sf_dz_ss(ss_Scalar_tmp) - DelTime*sf_VelZ*sf_avr_fs(fs_VelX)*sf_avr_fs(fs_dx_ss(ss_Scalar_tmp))) / (MAX(sf_avr_ss(ss_Scalar_tmp), epsilon(0.0d0))) - 0.5d0*sf_VelZ*sf_avr_ss(ss_DivVel)
end where
!=== 反拡散速度による過剰な輸送を抑制
! 反拡散速度が実際の速度を越えている場合は値を 0 にする.
where(abs(fs_AntiDiffVelX) > abs(fs_VelX) )
fs_AntiDiffVelX = 0.0d0
end where
where(abs(sf_AntiDiffVelZ) > abs(sf_VelZ) )
sf_AntiDiffVelZ = 0.0d0
end where
!=== 反拡散速度と上流差分法を用いてスカラーのフラックスを計算
call FluxScalar(fs_AntiDiffVelX, sf_AntiDiffVelZ, ss_Scalar_tmp, fs_AntiDiffFluxX, sf_AntiDiffFluxZ)
!=== 反拡散速度の移流による増分を計算
ss_AdvScalar_tmp = ss_dx_fs(fs_AntiDiffFluxX) + ss_dz_sf(sf_AntiDiffFluxZ)
!=== 全増分を計算
ss_AdvScalar = ss_AdvScalar_tmp + ss_AdvScalar
call EndSub("AdvectScalar_MPDATA")
end subroutine AdvectScalar_MPDATA