!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2004, 2006. All rights reserved. !--------------------------------------------------------------------- != Subroutine DynImpFunc ! ! * Developer: SUGIYAMA Ko-ichiro ! * Version: $Id: exner.f90,v 1.2.4.1 2005/08/30 13:41:00 kitamo Exp $ ! * Tag Name: $Name: arare3j $ ! * Change History: ! !== Overview ! !陰解法を用いたエクスナー関数の計算. !deepconv/arare では時間積分として He-VI 法を利用しているので, !エクスナー関数は陰解法で解く. ! !== Error Handling ! !== Known Bugs ! !== Note ! ! * 離散化する際, 上下境界条件として鉛直速度が零を与えている. ! * 空間方向に 2 次精度の離散化を陽に利用しているため, differentiate_center4 モジュールを指定することはできないので注意. ! !== Future Plans ! module DynImpFunc ! !陰解法を用いた力学過程の各項の計算モジュール. !エクスナー関数を陰解法で解く際に必要となる, 係数行列の要素の決定, !係数行列の LU 分解, LAPACK の関数のラッパーも提供. ! !モジュール読み込み use dc_message, only: MessageNotify !メッセージ出力 use gridset, only: DelZ, &! z 方向の格子点間隔 & DimXMin, &! x 方向の配列の下限 & DimXMax, &! x 方向の配列の上限 & DimZMin, &! z 方向の配列の下限 & DimZMax, &! z 方向の配列の上限 & RegXMin, &! x 方向の物理領域の下限 & RegXMax, &! x 方向の物理領域の上限 & RegZMin, &! z 方向の物理領域の下限 & RegZMax ! z 方向の物理領域の上限 use timeset, only: DelTimeShort !短い時間ステップ use damping, only: DampSound !音波の減衰係数 use basicset, only: CpDry, &!乾燥成分の比熱 & xz_VelSoundBasicZ, &!基本場の音速 & xz_DensBasicZ, &!基本場の密度 & xz_PotTempBasicZ, &!基本場の温位 & xz_EffMolWtBasicZ !基本場の分子量効果 use average, only: xr_avr_xz use differentiate_center2, only: xr_dz_xz, xz_dz_xr, xz_dx_pz !暗黙の型宣言禁止 implicit none !属性の指定 private !関数を public に設定 public DynImpFunc_init !初期化ルーチン public xz_Exner !エクスナー関数の計算 public xr_GradPi real(8) :: beta = 5.0d-1 !クランクニコルソン法なら 0.5 !完全陰解法なら 1 integer :: LN !配列のサイズ real(8), allocatable :: A(:) !係数行列 real(8), allocatable :: B(:) !係数行列 real(8), allocatable :: C(:) !係数行列 real(8), allocatable :: B2(:) !係数行列(LU 分解後) real(8), allocatable :: IPIV(:) !LU 分解の... real(8), allocatable :: xz_F1BasicZ(:,:) !係数行列の計算に用いる配列 real(8), allocatable :: xr_F2BasicZ(:,:) !係数行列の計算に用いる配列 real(8), allocatable :: xz_VPotTempBasicZ(:,:) !基本場の仮温位 !値の保存 save beta, LN, A, B, C, B2, IPIV save xz_F1BasicZ, xr_F2BasicZ, xz_VPotTempBasicZ contains !!!--------------------------------------------------------------------!!! subroutine DynImpFunc_init() ! !エクスナー関数を陰解法で解く際に必要となる, 係数行列の要素を決め, !LU 分解を行う. ! !暗黙の型宣言禁止 implicit none !配列の割り付け allocate( & & A(RegZMin+1:RegZMax), & & B(RegZMin+2:RegZMax), & & C(RegZMin+1:RegZMax-1), & & B2(RegZMin+3:RegZMax), & & xz_F1BasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & xr_F2BasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & xz_VPotTempBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & IPIV(RegZMin+1:RegZMax) ) !仮温位の定義 xz_VPotTempBasicZ = xz_PotTempBasicZ / xz_EffMolWtBasicZ !係数行列の計算 xz_F1BasicZ = ( xz_VelSoundBasicZ ** 2.0d0 ) & & / (CpDry * xz_DensBasicZ *(xz_VPotTempBasicZ ** 2.0d0)) xr_F2BasicZ = xr_avr_xz( & & CpDry * xz_DensBasicZ * (xz_VPotTempBasicZ ** 2.0d0) ) A(RegZMin+2: RegZMax-1) = & & 1.0d0 + (beta ** 2.0d0) & & * xz_F1BasicZ(RegXMax, RegZMin+2: RegZMax-1) & & * ( xr_F2BasicZ(RegXMax, RegZMin+2: RegZMax-1) & & + xr_F2BasicZ(RegXMax, RegZMin+1: RegZMax-2) ) & & * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0) A(RegZMin+1) = 1.0d0 & & + (beta ** 2.0d0) & & * xz_F1BasicZ(RegXMax, RegZMin+1) & & * xr_F2BasicZ(RegXMax, RegZMin+1) & & * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0) A(RegZMax) = 1.0d0 & & + (beta ** 2.0d0) & & * xz_F1BasicZ(RegXMax, RegZMax) & & * xr_F2BasicZ(RegXMax, RegZMax-1) & & * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0) B(RegZMin+2: RegZMax) = & & - ( beta ** 2.0d0 ) & & * xz_F1BasicZ(RegXMax, RegZMin+1: RegZMax-1) & & * xr_F2BasicZ(RegXMax, RegZMin+1: RegZMax-1) & & * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0 ) C(RegZMin+1: RegZMax-1) = & & - ( beta ** 2.0d0 ) & & * xz_F1BasicZ(RegXMax, RegZMin+2: RegZMax) & & * xr_F2BasicZ(RegXMax, RegZMin+1: RegZMax-1) & & * (DelTimeShort ** 2.0d0) / ( DelZ ** 2.0d0 ) !配列の大きさを保管 LN = RegZMax - RegZMin !係数行列を LU 分解 (LAPACK) call LinSolvLU(A, B, C, B2, IPIV) end subroutine DynImpFunc_init !!!--------------------------------------------------------------------!!! function xz_Exner(xr_FzNl, pz_VelXNs, pz_VelXAs, xr_VelZNs, xz_ExnerNs) ! !陰解法を用いたエクスナー関数の計算. ! !暗黙の型宣言禁止 implicit none !入出力変数 real(8), intent(in) :: pz_VelXNs(DimXMin:DimXMax, DimZMin:DimZMax) !速度 u [τ] real(8), intent(in) :: pz_VelXAs(DimXMin:DimXMax, DimZMin:DimZMax) !速度 u [τ+Δτ] real(8), intent(in) :: xr_VelZNs(DimXMin:DimXMax, DimZMin:DimZMax) !速度 w [τ] real(8), intent(in) :: xr_FzNl(DimXMin:DimXMax, DimZMin:DimZMax) !Z 方向の外力項[t] real(8), intent(in) :: xz_ExnerNs(DimXMin:DimXMax, DimZMin:DimZMax) !無次元圧力 real(8) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) !無次元圧力[τ+Δτ] !変数定義 real(8) :: D(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: E(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: F(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_DivVelNs(DimXMin:DimXMax, DimZMin:DimZMax) integer :: i !速度の収束を計算 xz_DivVelNs = xz_dx_pz( pz_VelXNs ) + xz_dz_xr( xr_VelZNs ) !行列計算のための係数 E = xr_dz_xz( DampSound * xz_DivVelNs ) & & - ( 1.0d0 - beta ) * xr_dz_xz( xz_ExnerNs ) & & + xr_FzNl / xr_avr_xz( CpDry * xz_VPotTempBasicZ ) F = - beta * xz_F1BasicZ * DelTimeShort & & * xz_dz_xr( xr_avr_xz( xz_DensBasicZ * xz_VPotTempBasicZ) & & * ( & & xr_VelZNs & & - xr_avr_xz(CpDry * xz_VPotTempBasicZ) * DelTimeShort & & * ( (1.0d0 - beta) * xr_dz_xz( xz_ExnerNs ) & & - xr_dz_xz( DampSound * xz_DivVelNs ) ) & & + xr_FzNl * DelTimeShort & & ) & & ) D = xz_ExnerNs & & - (1.0d0 - beta) & & * xz_F1BasicZ * DelTimeShort & & * xz_dz_xr( & & xr_avr_xz(xz_DensBasicZ * xz_VPotTempBasicZ) * xr_VelZNs & & ) & & - (xz_VelSoundBasicZ ** 2.0d0) * DelTimeShort & & / (CpDry * xz_VPotTempBasicZ) & & * xz_dx_pz( pz_VelXAs ) & & + F D(:, RegZMin+1) = D(:, RegZMin+1) & & - beta * xz_F1BasicZ(:, RegZMin+1) * ( DelTimeShort ** 2.0d0 ) & & * xr_F2BasicZ(:, RegZMin) * E(:,RegZMin) & & / DelZ D(:, RegZMax) = D(:, RegZMax) & & + beta * xz_F1BasicZ(:, RegZMax) * ( DelTimeShort ** 2.0d0 ) & & * xr_F2BasicZ(:, RegZMax) * E(:, RegZMax) & & / DelZ !連立一次方程式の解を求める do i = RegXMin + 1, RegXMax call LinSolv( D(i, RegZMin+1: RegZMax) ) end do !D の戻り値を出力 xz_Exner = D end function xz_Exner !!!--------------------------------------------------------------------!!! subroutine LinSolv(D) ! !LU 分解されている実 3 項行列の連立 1 次方程式(倍精度). LAPACK 利用 ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(inout) :: D(1,LN) !定数/解行列 integer :: NRHS integer :: INFO integer :: LDB character(1),parameter :: TRANS = 'N' !変数の初期化 NRHS = 0 INFO = 0 LDB = LN !解行列の計算. LAPACK を使用. call DGTTRS(TRANS, LN, NRHS, C, A, B, B2, IPIV, D, LDB, INFO) ! !解のコンディションをチェック. ! if (INFO /= 0) then ! call MessageNotify("Error", "lapack_linear", "INFO is not 0") ! stop ! end if end subroutine LinSolv !!!--------------------------------------------------------------------!!! subroutine LinSolvLU(A, B, C, B2, IPIV) ! !実 3 項行列の LU 分解(倍精度). LAPACK 利用 ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(inout) :: A(LN) !係数行列 real(8), intent(inout) :: B(LN-1) !係数行列 real(8), intent(inout) :: C(LN-1) !係数行列 real(8), intent(out) :: B2(LN-2) !係数行列 real(8), intent(out) :: IPIV(LN) !部分ピボット交換の情報を格納 integer :: INFO !変数の初期化 INFO = 0 !解行列の計算. LAPACK を使用. call DGTTRF(LN, C, A, B, B2, IPIV, INFO) ! !解のコンディションをチェック. ! if (INFO /= 0) then ! call MessageNotify("Error", "lapack_linear", "INFO is not 0") ! stop ! end if end subroutine LinSolvLU !!!--------------------------------------------------------------------!!! function xr_GradPi(xz_ExnerA, xz_ExnerN, pz_VelX, xr_VelZ) ! ! x 方向に半格子ずれた点での圧力傾度力項の計算. ! クランク・ニコルソン法を用いるために, 時刻 τ と τ+Δτでの ! エクスナー関数の値を引数として取る. ! 音波減衰項を含めた形式で定式化してあることに注意. ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: xz_ExnerA(DimXMin:DimXMax, DimZMin:DimZMax) !エクスナー関数の擾乱 real(8), intent(in) :: xz_ExnerN(DimXMin:DimXMax, DimZMin:DimZMax) !エクスナー関数の擾乱 real(8), intent(in) :: pz_VelX(DimXMin:DimXMax, DimZMin:DimZMax) !水平速度 real(8), intent(in) :: xr_VelZ(DimXMin:DimXMax, DimZMin:DimZMax) !鉛直速度 real(8) :: xr_GradPi(DimXMin:DimXMax, DimZMin:DimZMax) !圧力傾度力 real(8) :: xz_DivVel(DimXMin:DimXMax, DimZMin:DimZMax) !速度の収束 !速度の収束 xz_DivVel = xz_dx_pz( pz_VelX ) + xz_dz_xr( xr_VelZ ) !速度 w の圧力勾配 xr_GradPi = 0.0d0 xr_GradPi = & & xr_avr_xz(CpDry * xz_PotTempBasicZ / xz_EffMolWtBasicZ ) & & * ( & & beta * xr_dz_xz( xz_ExnerA ) & & + (1.0d0 - beta) * xr_dz_xz( xz_ExnerN ) & & - xr_dz_xz( DampSound * xz_DivVel ) & & ) end function xr_GradPi end module DynImpFunc