!c Description: !c HE-VI 法を用いたエクスナー関数の計算. 陰解法で計算 !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) !--- モジュールの読み込み 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 !--- 暗黙の型宣言禁止 implicit none !--- 入出力変数 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 !共通に用いる変数 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) ) !--- 行列計算のための係数 !----- 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) ) & & * (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 !--- 連立一次方程式の解を求める do i = RegXMin + 1, RegXMax call LinSolv( A, B, C, D(i, RegZMin+1: RegZMax) ) end do !--- 計算結果を入力 ss_Exner_A = D !--- 境界条件 call boundary(ss_Bc, ss_Exner_A) end subroutine ExnerHeVi