    Copyright (C) GFD Dennou Club, 2004. All rights reserved.

begin

Subroutine EddyViscosity

  * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
  * Version: $Id: eddyviscosity.f90,v 1.2 2005/04/22 15:01:52 sugiyama Exp $
  * Tag Name: $Name:  $
  * Change History:

Overview

2.5 次の乱流クロージャーを利用し, 次の時刻の渦粘性係数と渦拡散係数を求めるサブルーチン.

Error Handling

Known Bugs

Note

 * 定式化は CReSS のマニュアルを参照した
 * 移流項の計算には ***_adv となる変数を利用する.
 * 生成項, 拡散項, 消散項の計算には ***_dif となる変数を利用する.

Future Plans

end

Methods

EddyViscosity

Public Instance methods

EddyViscosity( ss_Km_in, DelTime, fs_VelX_adv, sf_VelZ_adv, ss_Km_adv, fs_VelX_dif, sf_VelZ_dif, ss_PotTemp_dif, ss_Km_dif, ss_Km_out, ss_Kh_out )

ss_Km_in(DimXMin:DimXMax, DimZMin:DimZMax) : real(8), intent(in)

DelTime :

real(8), intent(in)

:	end 暗黙の型宣言禁止 begin Input

fs_VelX_adv(DimXMin:DimXMax, DimZMin:DimZMax) : real(8), intent(in)

sf_VelZ_adv(DimXMin:DimXMax, DimZMin:DimZMax) : real(8), intent(in)

ss_Km_adv(DimXMin:DimXMax, DimZMin:DimZMax) : real(8), intent(in)

fs_VelX_dif(DimXMin:DimXMax, DimZMin:DimZMax) : real(8), intent(in)

sf_VelZ_dif(DimXMin:DimXMax, DimZMin:DimZMax) : real(8), intent(in)

ss_PotTemp_dif(DimXMin:DimXMax, DimZMin:DimZMax) : real(8), intent(in)

ss_Km_dif(DimXMin:DimXMax, DimZMin:DimZMax) : real(8), intent(in)

ss_Km_out(DimXMin:DimXMax, DimZMin:DimZMax) :

real(8), intent(out)

Output

ss_Kh_out(DimXMin:DimXMax, DimZMin:DimZMax) : real(8), intent(out)

[Source]

subroutine EddyViscosity(                                    ss_Km_in,                                             DelTime,                                              fs_VelX_adv, sf_VelZ_adv, ss_Km_adv,                  fs_VelX_dif, sf_VelZ_dif, ss_PotTemp_dif, ss_Km_dif,  ss_Km_out, ss_Kh_out )
                                                                 !=begin
  !== Dependency

                                                                 !=end
  !== 暗黙の型宣言禁止
  implicit none
                                                                 !=begin
  !== Input
  real(8), intent(in)  :: DelTime
  real(8), intent(in)  :: ss_Km_in(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: fs_VelX_adv(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: sf_VelZ_adv(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_Km_adv(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: fs_VelX_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: sf_VelZ_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_PotTemp_dif(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)  :: ss_Km_dif(DimXMin:DimXMax, DimZMin:DimZMax)

  !== Output
  real(8), intent(out) :: ss_Km_out(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(out) :: ss_Kh_out(DimXMin:DimXMax, DimZMin:DimZMax)
                                                                 !=end
  !== Work
  real(8)  :: ss_AdvKm(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_BuoyKm(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_ShearKm(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_DiffKm(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_DispKm(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)  :: ss_TendKm(DimXMin:DimXMax, DimZMin:DimZMax)
  

  call BeginSub("EddyViscosity",         fmt="%c",            c1="calculation of eddy viscosity whith 2.5 order closure.")


  !=== スカラー量の移流. 中心差分. 
  ss_AdvKm =    ss_avr_fs( fs_VelX_adv * fs_dx_ss( ss_Km_adv ) )  + ss_avr_sf( sf_VelZ_adv * sf_dz_ss( ss_Km_adv ) ) 


  !=== 浮力による乱流エネルギー生成
  ss_BuoyKm =     3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( DelXZ ** 2.0d0 )        * ss_avr_sf( sf_dz_ss( ss_PotTemp_dif ) )        / ( 2.0d0 * ss_PotTempBasicZ )

  
  !=== 速度シアーによる乱流エネルギー生成
  ss_ShearKm =   ( Cm ** 2.0d0 ) * ( DelXZ ** 2.0d0 )                   * (                                                        ( ss_dx_fs( fs_VelX_dif ) ) ** 2.0d0               + ( ss_dz_sf( sf_VelZ_dif ) ) ** 2.0d0                + ( (   ss_avr_ff( ff_dz_fs( fs_VelX_dif ) )                + ss_avr_ff( ff_dx_sf( sf_VelZ_dif ) )              ) ** 2.0d0 ) / 2.0d0                            )                                               - ss_Km_dif                                              * ( ss_dx_fs( fs_VelX_dif ) + ss_dz_sf( sf_VelZ_dif ) ) / 3.0d0

  
  !=== 乱流エネルギーの拡散
  ss_DiffKm =                                                      ss_dx_fs( fs_dx_ss( ss_Km_dif ) ) / 2.0d0           + ss_dz_sf( sf_dz_ss( ss_Km_dif ) ) / 2.0d0           + ( ss_avr_fs( fs_dx_ss( ss_Km_dif ) ) ) ** 2.0d0     + ( ss_avr_sf( sf_dz_ss( ss_Km_dif ) ) ) ** 2.0d0
  

  !=== 乱流エネルギーの消散
  ss_DispKm = ( ss_Km_dif ** 2.0d0 ) / ( 2.0d0 * ( DelXZ ** 2.0d0 ) )
  
  
  !=== 項をまとめる
  ss_TendKm = - ss_AdvKm  - ss_BuoyKm  +  ss_ShearKm       + ss_DiffKm - ss_DispKm 

  
  !=== 時間積分
  ss_Km_out = ss_Km_in + DelTime * ss_TendKm


  !=== 境界条件
  call boundary( ss_BC, ss_Km_out )
  

  !=== 負の領域を零に固定
  where ( ss_Km_out < 0.0d0 )
     ss_Km_out = 0.0d0
  end where

  
  !=== 渦拡散定数を決める
  ss_Kh_out = 3.0d0 * ss_Km_out


  call EndSub("EddyViscosity")
  
end subroutine EddyViscosity

eddyviscosity.f90