!c Description: 
!c   HE-VI $BK!$rMQ$$$?%(%/%9%J!<4X?t$N7W;;(B. $B1"2rK!$G7W;;(B
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2004, All rights reserved

subroutine ExnerHeVi(&
     & ss_VelDiv_B, sf_Fz_B, fs_VelX_A, sf_VelZ_B, ss_Exner_B, ss_Exner_A)

  !--- $B%b%8%e!<%k$NFI$_9~$_(B
  use gridset
  use timeset,  only: beta, DelTShort
  use basicprm, only: ss_VelSoundBasicZ, ss_CpBasicZ, &
       & ss_DensBasicZ, ss_VThetaBasicZ, alpha
  use if_display
  use if_avr
  use if_diff
  use if_boundary
  use linlib,  only: LinSolv
  
  !--- $B0EL[$N7?@k8@6X;_(B
  implicit none
  
  !--- $BF~=PNOJQ?t(B
  real(8), intent(in)      :: ss_VelDiv_B(DimXMin:DimXMax, DimZMin:DimZMax)  
  real(8), intent(in)      :: fs_VelX_A(DimXMin:DimXMax, DimZMin:DimZMax)  
  real(8), intent(in)      :: sf_VelZ_B(DimXMin:DimXMax, DimZMin:DimZMax)  
  real(8), intent(in)      :: sf_Fz_B(DimXMin:DimXMax, DimZMin:DimZMax)  
  real(8), intent(in)      :: ss_Exner_B(DimXMin:DimXMax, DimZMin:DimZMax)  
  real(8), intent(out)     :: ss_Exner_A(DimXMin:DimXMax, DimZMin:DimZMax)  

  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


  !$B6&DL$KMQ$$$kJQ?t(B
  ss_F1BasicZ = ( ss_VelSoundBasicZ ** 2.0d0 ) &
       &        / (ss_CpBasicZ * ss_DensBasicZ * (ss_VthetaBasicZ ** 2.0d0) )

  sf_F2BasicZ = sf_avr_ss( &
       & ss_CpBasicZ * ss_DensBasicZ * (ss_VthetaBasicZ ** 2.0d0) )


  !--- $B9TNs7W;;$N$?$a$N78?t(B
  !----- BasicZ $B$JG[Ns$O(B X $BJ}8~$O0lMM$J$?$a(B, $B0l<!85G[Ns$H$7$F07$&$?$a$K(B
  !----- RegXMax $B$GBeI=$5$;$F$"$k(B
  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) ) &
       &   * (DelTShort ** 2.0d0) / ( DelZ ** 2.0d0)
  
  A(RegZMin+1) = 1.0d0 &
       & + (beta ** 2.0d0) &
       &   * ss_F1BasicZ(RegXMax, RegZMin+1)  &
       &   * sf_F2BasicZ(RegXMax, RegZMin+1) &
       &   * (DelTShort ** 2.0d0) / ( DelZ ** 2.0d0)

  A(RegZMax) = 1.0d0 &
       & + (beta ** 2.0d0) &
       &   * ss_F1BasicZ(RegXMax, RegZMax)  &
       &   * sf_F2BasicZ(RegXMax, RegZMax-1) &
       &   * (DelTShort ** 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) &
       &   * (DelTShort ** 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) &
       &   * (DelTShort ** 2.0d0) / ( DelZ ** 2.0d0 )
  
  E = sf_dz_ss( alpha * ss_VelDiv_B ) &
       & - ( 1.0d0 - beta ) * sf_dz_ss( ss_Exner_B ) &
       & + sf_Fz_B / sf_avr_ss( ss_CpBasicZ * ss_VThetaBasicZ )

  D = ss_Exner_B &
       & - (1.0d0 - beta) &
       &   * ss_F1BasicZ  &
       &   * ss_dz_sf(    &
       &       sf_avr_ss( ss_DensBasicZ * ss_VThetaBasicZ ) * sf_VelZ_B ) &
       &   * DelTShort &
       & - ( ss_F1BasicZ * ss_DensBasicZ * ss_VThetaBasicZ ) &
       &   * ss_dx_fs( fs_VelX_A ) * DelTShort  &
       & - beta &
       &   * ss_F1BasicZ &
       &   * ss_dz_sf( sf_avr_ss( ss_DensBasicZ * ss_VThetaBasicZ) &
       &     * ( sf_VelZ_B &
       &         - sf_avr_ss( ss_CpBasicZ * ss_VThetaBasicZ) * DelTShort &
       &           * ( (1.0d0 - beta) * sf_dz_ss( ss_Exner_B ) &
       &               - sf_dz_ss( alpha * ss_VelDiv_B ) )  &
       &         + sf_Fz_B * DelTShort &
       &        ) &
       &     ) * DelTShort
  
  D(:, RegZMin+1) = D(:, RegZMin+1) &
       & - beta * ss_F1BasicZ(:, RegZMin+1) &
       &   * sf_F2BasicZ(:, RegZMin) &
       &   * E(:,RegZMin) &
       &   * ( DelTShort ** 2.0d0 ) / DelZ

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


  !--- $BO"N)0l<!J}Dx<0$N2r$r5a$a$k(B
  do i = RegXMin + 1, RegXMax
     call LinSolv( A, B, C, D(i, RegZMin+1: RegZMax) )
  end do
  
  !--- $B7W;;7k2L$rF~NO(B
  ss_Exner_A = D
  
  !--- $B6-3&>r7o(B
  call boundary(ss_Bc, ss_Exner_A)
  
end subroutine ExnerHeVi
