!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine Scalar
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: scalar.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!スカラー量の時間発展を解く
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * 移流項の計算には ***_adv となる変数を利用する. 
!  * 生成項, 乱流拡散項, 数値粘性項の計算には ***_dif となる変数を利用する. 
!
!== Future Plans
!
                                                                 !=end
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  
  use dc_trace, only: BeginSub, EndSub
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax
  use bcset,    only: ss_BC
  use average,  only: ss_avr_sf
  use basicset, only: ss_PotTempBasicZ
  use arareset, only: NuH, NuV
  use average,  only: ss_avr_fs, ss_avr_sf, fs_avr_ss, sf_avr_ss
  use differentiate_center4, only: fs_dx_ss, ss_dx_fs, sf_dz_ss, ss_dz_sf
                                                                 !=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
  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 bcset,         only: ss_BC, sf_BC, fs_BC
  use differentiate_center4, only: fs_dx_ss, sf_dz_ss, ss_dx_fs, ss_dz_sf
                                                                 !=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

!!!
!!! 上流差分法を用いてスカラーのフラックスを計算
!!!
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 Scalar
