!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=Begin
!= Subroutine Exner
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: exner.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!陰解法を用いたエクスナー関数の計算. 
!deepconv/arare では時間積分として He-VI 法を利用しているので, 
!エクスナー関数は陰解法で解く.
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * 離散化する際, 上下境界条件として鉛直速度が零を与えている. 
!  * 空間方向に 2 次精度の離散化を陽に利用しているため, differentiate_center4 モジュールを指定することはできないので注意.
!
!== Future Plans
!
                                                                 !=end

subroutine Exner( &
     & sf_Fz_nl, fs_VelX_ns, fs_VelX_as, sf_VelZ_ns, ss_Exner_ns, &
     & ss_Exner_as )
                                                                 !=begin
  !== Dependency
  use dc_trace, only: BeginSub, EndSub  
  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, &
       &              RegXMax, RegXMin, RegZMin, RegZMax, DelZ
  use timeset,  only: DelTimeShort
  use basicset, only: ss_VelSoundBasicZ, ss_CpBasicZ, &
       &              ss_DensBasicZ, ss_PotTempBasicZ
  use bcset,    only: ss_BC
  use arareset, only: alpha, beta
  use linlib,   only: LinSolv
  use average,  only: sf_avr_ss
  use differentiate_center2, only: sf_dz_ss, ss_dz_sf, ss_dx_fs
                                                                 !=end  
  !== 暗黙の型宣言禁止
  implicit none
                                                                 !=begin  
  !== Input
  real(8), intent(in)      :: fs_VelX_ns(DimXMin:DimXMax, DimZMin:DimZMax)    !速度 u
  real(8), intent(in)      :: fs_VelX_as(DimXMin:DimXMax, DimZMin:DimZMax)    !速度 u
  real(8), intent(in)      :: sf_VelZ_ns(DimXMin:DimXMax, DimZMin:DimZMax)    !速度 w
  real(8), intent(in)      :: sf_Fz_nl(DimXMin:DimXMax, DimZMin:DimZMax)      !Z 方向の外力項
  real(8), intent(in)      :: ss_Exner_ns(DimXMin:DimXMax, DimZMin:DimZMax)   !無次元圧力

  !== Output
  real(8), intent(out)     :: ss_Exner_as(DimXMin:DimXMax, DimZMin:DimZMax)   !無次元圧力 (1 ステップ後)
                                                                 !=end
  !== Work
  real(8)                  :: D(DimXMin:DimXMax, DimZMin:DimZMax)  
  real(8)                  :: E(DimXMin:DimXMax, DimZMin:DimZMax)  
  real(8)                  :: ss_F1BasicZ(DimXMin:DimXMax, DimZMin:DimZMax)  
  real(8)                  :: sf_F2BasicZ(DimXMin:DimXMax, DimZMin:DimZMax)  
  real(8)                  :: ss_DivVel_ns(DimXMin:DimXMax, DimZMin:DimZMax)
  integer                  :: i


  call BeginSub("Exner", &
&           fmt="%c",        &
&           c1="Calculate Exner with HE-VI (ss_Exner_as).")


  !=== 速度の div を計算
  ss_DivVel_ns =  ss_dx_fs( fs_VelX_ns ) + ss_dz_sf( sf_VelZ_ns )

  !=== 共通に用いる変数
  ss_F1BasicZ = ( ss_VelSoundBasicZ ** 2.0d0 ) &
       &        / (ss_CpBasicZ * ss_DensBasicZ *(ss_PotTempBasicZ ** 2.0d0))
  
  sf_F2BasicZ = sf_avr_ss( &
       & ss_CpBasicZ * ss_DensBasicZ * (ss_PotTempBasicZ ** 2.0d0) )


  !=== 行列計算のための係数
  !==== BasicZ な配列は X 方向は一様なため, 一次元配列として扱うために
  !==== RegXMax で代表させてある
  
  E = sf_dz_ss( alpha * ss_DivVel_ns ) &
       & - ( 1.0d0 - beta ) * sf_dz_ss( ss_Exner_ns ) &
       & + sf_Fz_nl / sf_avr_ss( ss_CpBasicZ * ss_PotTempBasicZ )

  D = ss_Exner_ns &
       & - (1.0d0 - beta) &
       &   * ss_F1BasicZ  &
       &   * ss_dz_sf(    &
       &       sf_avr_ss( ss_DensBasicZ * ss_PotTempBasicZ ) * sf_VelZ_ns ) &
       &   * DelTimeShort &
       & - ( ss_F1BasicZ * ss_DensBasicZ * ss_PotTempBasicZ ) &
       &   * ss_dx_fs( fs_VelX_as ) * DelTimeShort  &
       & - beta &
       &   * ss_F1BasicZ &
       &   * ss_dz_sf( sf_avr_ss( ss_DensBasicZ * ss_PotTempBasicZ) &
       &     * ( sf_VelZ_ns &
       &         - sf_avr_ss(ss_CpBasicZ * ss_PotTempBasicZ) * DelTimeShort &
       &           * ( (1.0d0 - beta) * sf_dz_ss( ss_Exner_ns ) &
       &               - sf_dz_ss( alpha * ss_DivVel_ns ) )  &
       &         + sf_Fz_nl * DelTimeShort &
       &        ) &
       &     ) * DelTimeShort
  
  D(:, RegZMin+1) = D(:, RegZMin+1) &
       & - beta * ss_F1BasicZ(:, RegZMin+1) &
       &   * sf_F2BasicZ(:, RegZMin) &
       &   * E(:,RegZMin) &
       &   * ( DelTimeShort ** 2.0d0 ) / DelZ

  D(:, RegZMax) = D(:, RegZMax) &
       & + beta * ss_F1BasicZ(:, RegZMax) &
       &   * sf_F2BasicZ(:, RegZMax) &
       &   * E(:, RegZMax) &
       &   * ( DelTimeShort ** 2.0d0 ) / DelZ


  !=== 連立一次方程式の解を求める
  do i = RegXMin + 1, RegXMax
     call LinSolv( D(i, RegZMin+1: RegZMax) )
  end do
  
  !=== 計算結果を入力
  ss_Exner_as = D

  !=== 境界条件
  call boundary(ss_BC, ss_Exner_as)

  call EndSub("Exner")
  
end subroutine Exner
