exner.f90

Path: src/dynamics/exner.f90
Last Update: Tue Aug 30 21:24:53 JST 2005 
    Copyright (C) GFD Dennou Club, 2004. All rights reserved.

Begin

Subroutine Exner

  * Developer: SUGIYAMA Ko-ichiro
  * Version: $Id: exner.f90,v 1.3 2005/08/30 12:24:53 kitamo Exp $
  * Tag Name: $Name:  $
  * Change History:

Overview

陰解法を用いたエクスナー関数の計算. deepconv/arare では時間積分として He-VI 法を利用しているので, エクスナー関数は陰解法で解く.

Error Handling

Known Bugs

Note

 * 離散化する際, 上下境界条件として鉛直速度が零を与えている.
 * 空間方向に 2 次精度の離散化を陽に利用しているため, differentiate_center4 モジュールを指定することはできないので注意.

Future Plans

end

Methods

Exner  

Included Modules

dc_trace gridset timeset basicset bcset arareset linlib average differentiate_center2

Public Instance methods

Exner( sf_Fz_nl, fs_VelX_ns, fs_VelX_as, sf_VelZ_ns, ss_Exner_ns, ss_Exner_as )
sf_Fz_nl(DimXMin:DimXMax, DimZMin:DimZMax) :real(8), intent(in)
: Z 方向の外力項
fs_VelX_ns(DimXMin:DimXMax, DimZMin:DimZMax) :real(8), intent(in)
: end 
 暗黙の型宣言禁止

begin 

 Input

速度 u

fs_VelX_as(DimXMin:DimXMax, DimZMin:DimZMax) :real(8), intent(in)
: 速度 u
sf_VelZ_ns(DimXMin:DimXMax, DimZMin:DimZMax) :real(8), intent(in)
: 速度 w
ss_Exner_ns(DimXMin:DimXMax, DimZMin:DimZMax) :real(8), intent(in)
: 無次元圧力
ss_Exner_as(DimXMin:DimXMax, DimZMin:DimZMax) :real(8), intent(out)
: Output 無次元圧力 (1 ステップ後)

[Source]

subroutine Exner(  sf_Fz_nl, fs_VelX_ns, fs_VelX_as, sf_VelZ_ns, ss_Exner_ns,  ss_Exner_as )
                                                                 !=begin
  !== Dependency

                                                                 !=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

[Validate]