!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine Exner_HeVi
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: exner_hevi.f90,v 1.2 2005/01/28 12:41:21 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!HE-VI 法を用いたエクスナー関数の計算. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!離散化する際, 上下境界条件として鉛直速度が零を与えている. 
!
!== Future Plans
!
                                                                 !=end

subroutine Exner_HeVi(&
     & ss_DivVel_n, sf_Fz_n, fs_VelX_a, sf_VelZ_n, ss_Exner_n, ss_Exner_a)
                                                                 !=begin
  !== Dependency
  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 arareprm, 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)      :: ss_DivVel_n(DimXMin:DimXMax, DimZMin:DimZMax)  !速度の収束
  real(8), intent(in)      :: fs_VelX_a(DimXMin:DimXMax, DimZMin:DimZMax)    !速度 u
  real(8), intent(in)      :: sf_VelZ_n(DimXMin:DimXMax, DimZMin:DimZMax)    !速度 w
  real(8), intent(in)      :: sf_Fz_n(DimXMin:DimXMax, DimZMin:DimZMax)      !Z 方向の外力項
  real(8), intent(in)      :: ss_Exner_n(DimXMin:DimXMax, DimZMin:DimZMax)   !無次元圧力

  !== Output
  real(8), intent(out)     :: ss_Exner_a(DimXMin:DimXMax, DimZMin:DimZMax)   !無次元圧力 (1 ステップ後)
                                                                 !=end
  real(8)                  :: A(RegZMin+1:RegZMax)  
  real(8)                  :: B(RegZMin+2:RegZMax)  
  real(8)                  :: C(RegZMin+1:RegZMax-1)  
  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)  
  integer                  :: i


  !共通に用いる変数
  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 で代表させてある
  A(RegZMin+2: RegZMax-1) = &
       & 1.0d0 + (beta ** 2.0d0) &
       &   * ss_F1BasicZ(RegXMax, RegZMin+2: RegZMax-1)  &  
       &   * ( sf_F2BasicZ(RegXMax, RegZMin+2: RegZMax-1)  &
       &       +  sf_F2BasicZ(RegXMax, RegZMin+1: RegZMax-2) ) &
       &   * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0)
  
  A(RegZMin+1) = 1.0d0 &
       & + (beta ** 2.0d0) &
       &   * ss_F1BasicZ(RegXMax, RegZMin+1)  &
       &   * sf_F2BasicZ(RegXMax, RegZMin+1) &
       &   * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0)

  A(RegZMax) = 1.0d0 &
       & + (beta ** 2.0d0) &
       &   * ss_F1BasicZ(RegXMax, RegZMax)  &
       &   * sf_F2BasicZ(RegXMax, RegZMax-1) &
       &   * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0)
  
  B(RegZMin+2: RegZMax) = &
       & - ( beta ** 2.0d0 ) &
       &   * ss_F1BasicZ(RegXMax, RegZMin+1: RegZMax-1) &
       &   * sf_F2BasicZ(RegXMax, RegZMin+1: RegZMax-1) &
       &   * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0 )
  
  C(RegZMin+1: RegZMax-1) = &
       & - ( beta ** 2.0d0 ) &
       &   * ss_F1BasicZ(RegXMax, RegZMin+2: RegZMax) &
       &   * sf_F2BasicZ(RegXMax, RegZMin+1: RegZMax-1) &
       &   * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0 )
  
  E = sf_dz_ss( alpha * ss_DivVel_n ) &
       & - ( 1.0d0 - beta ) * sf_dz_ss( ss_Exner_n ) &
       & + sf_Fz_n / sf_avr_ss( ss_CpBasicZ * ss_PotTempBasicZ )

  D = ss_Exner_n &
       & - (1.0d0 - beta) &
       &   * ss_F1BasicZ  &
       &   * ss_dz_sf(    &
       &       sf_avr_ss( ss_DensBasicZ * ss_PotTempBasicZ ) * sf_VelZ_n ) &
       &   * DelTimeShort &
       & - ( ss_F1BasicZ * ss_DensBasicZ * ss_PotTempBasicZ ) &
       &   * ss_dx_fs( fs_VelX_a ) * DelTimeShort  &
       & - beta &
       &   * ss_F1BasicZ &
       &   * ss_dz_sf( sf_avr_ss( ss_DensBasicZ * ss_PotTempBasicZ) &
       &     * ( sf_VelZ_n &
       &         - sf_avr_ss(ss_CpBasicZ * ss_PotTempBasicZ) * DelTimeShort &
       &           * ( (1.0d0 - beta) * sf_dz_ss( ss_Exner_n ) &
       &               - sf_dz_ss( alpha * ss_DivVel_n ) )  &
       &         + sf_Fz_n * 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( A, B, C, D(i, RegZMin+1: RegZMax) )
  end do
  
  !--- 計算結果を入力
  ss_Exner_a = D
  
end subroutine Exner_HeVi
