!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine AdvectScalar_MPDATA
!
!   * Developer: ODAKA Masatsugu (odakker@gfd-dennou.org)
!   * Version: $Id: advectscalar_mpdata.f90,v 1.3 2005/02/18 06:51:38 odakker 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
! を時間積分する. 2005/02/18 現在, 本数値コードではこちらで実装している.
!
!       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
  use dc_trace,      only: BeginSub, EndSub 
  use gridset,       only: DimXMin, DimXMax, DimZMin, DimZMax, DelX, DelZ
  use average,       only: fs_avr_sf, sf_avr_fs, sf_avr_ss, fs_avr_ss
  use differentiate, only: fs_dx_ss, sf_dz_ss, ss_dx_fs, ss_dz_sf
  use bcset,         only: ss_BC, sf_BC, fs_BC
                                                                 !=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
 
  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
 

  !--- 上流差分法を用いてスカラーのフラックスを計算

  call FluxScalar(fs_VelX, sf_VelZ, ss_Scalar, fs_FluxX, sf_FluxZ)

  !--- 発散の計算

  call DivergenceVel(fs_VelX, sf_VelZ, ss_DivVel)

  !--- 移流による増分を計算

  ss_AdvScalar = ss_dx_fs(fs_FluxX) + ss_dz_sf(sf_FluxZ)  
!&                - ss_Scalar * ss_DivVel

  !--- スカラーの仮値を計算

  call Integral(- ss_AdvScalar, DelTime, ss_Scalar, ss_Scalar_tmp )
  call boundary(ss_BC, ss_Scalar_tmp)


  !--- MPDATA スキームによる補正 ----------------------------------------
  !--- 上流差分によるスカラーフラックスが 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")

  contains

  !--- 上流差分法を用いてスカラーのフラックスを計算

  subroutine FluxScalar(fs_VelX, sf_VelZ, ss_Scalar, &
&                       fs_FluxX, sf_FluxZ)

    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(out) :: fs_FluxX(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(out) :: sf_FluxZ(DimXMin:DimXMax, DimZMin:DimZMax)
   

    fs_FluxX(DimXMin:DimXMax-1, :) =  (&
&      (fs_VelX(DimXMin:DimXMax-1, :) + abs(fs_VelX(DimXMin:DimXMax-1, :))) &
&              *ss_Scalar(DimXMin:DimXMax-1, :)       &
&     +(fs_VelX(DimXMin:DimXMax-1, :) - abs(fs_VelX(DimXMin:DimXMax-1, :))) &
&              *ss_Scalar(DimXMin+1:DimXMax, :) )*0.5d0

    sf_FluxZ(:, DimZMin:DimZMax-1) =  (& 
&       (sf_VelZ(:, DimZMin:DimZMax-1) + abs(sf_VelZ(:, DimZMin:DimZMax-1))) &
&              *ss_Scalar(:, DimZMin:DimZMax-1)       &
&     + (sf_VelZ(:, DimZMin:DimZMax-1) - abs(sf_VelZ(:, DimZMin:DimZMax-1))) &
&              *ss_Scalar(:, DimZMin+1:DimZMax) )*0.5d0


  end subroutine FluxScalar
  
end subroutine AdvectScalar_MPDATA

