!= Module Turbulence_3d ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA, Masatsugu ! Version:: $Id: turbulence_3d.f90,v 1.2 2007/10/01 08:57:59 odakker Exp $ ! Tag Name:: $Name: arare4-20071012 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! !モデルの物理過程を計算するために必要となる関数群を束ねたモジュール !具体的には以下の項を計算するための関数を格納する. ! * 乱流拡散項 ! * 乱流エネルギーの時間発展方程式に含まれる各項 ! * 散逸加熱項 ! !== Error Handling ! !== Bugs ! !== Note ! ! * 1.5 次のクロージャーを利用する ! * 定式化は CReSS のマニュアルを参照した ! !== Future Plans ! ! module Turbulence_3d ! !モデルの物理過程を計算するために必要となる関数群を束ねたモジュール !具体的には以下の項を計算するための関数を格納する. ! * 乱流拡散項 ! * 乱流エネルギーの時間発展方程式に含まれる各項 ! * 散逸加熱項 ! !モジュール読み込み use dc_types, only : DP use gridset_3d, only: DimXMin, &! x 方向の配列の下限 & DimXMax, &! x 方向の配列の上限 & DimYMin, &! y 方向の配列の下限 & DimYMax, &! y 方向の配列の上限 & DimZMin, &! z 方向の配列の下限 & DimZMax, &! z 方向の配列の上限 & SpcNum, &! z 方向の配列の上限 & x_dx, &! x 方向の格子点間隔 & y_dy, &! y 方向の格子点間隔 & z_dz ! z 方向の格子点間隔 use basicset_3d, only: CpDry, &!乾燥成分の比熱 & Grav, &!重力加速度 & xyz_PotTempBasicZ, &!基本場の温位 & xyz_ExnerBasicZ !エクスナー関数の基本場 use xyz_module, only: & & xyz_avr_pyz, xyr_avr_pyr, xqz_avr_pqz, & & pyz_avr_xyz, pyr_avr_xyr, pqz_avr_xqz, & & xyz_avr_xqz, pyz_avr_pqz, xyr_avr_xqr, & & xqz_avr_xyz, pqz_avr_pyz, xqr_avr_xyr, & & xyz_avr_xyr, pyz_avr_pyr, xqz_avr_xqr, & & xyr_avr_xyz, pyr_avr_pyz, xqr_avr_xqz, & & pqz_avr_xyz, pyr_avr_xyz, xqr_avr_xyz, & & xyz_avr_pqz, xyz_avr_pyr, xyz_avr_xqr use xyz_deriv_module, only: xyz_dx_pyz, xyz_dy_xqz, xyz_dz_xyr, & & pyz_dx_xyz, xqz_dy_xyz, xyr_dz_xyz, & & pqz_dx_xqz, pqz_dy_pyz, pyz_dy_pqz, & & pyr_dx_xyr, pyz_dz_pyr, pyr_dz_pyz, & & pyr_dz_pyz, pyr_dx_xyr, xyr_dx_pyr, & & xqr_dz_xqz, xqr_dy_xyr, xyr_dy_xqr, & & pqz_dy_pyz, pqz_dx_xqz, xqz_dx_pqz, & & xqr_dy_xyr, xqr_dz_xqz, xqz_dz_xqr use StorePotTemp_3d, only: StorePotTempDisp, StorePotTempTurb ! use StoreMixRt, only: StoreMixRtTurb !暗黙の型宣言禁止 implicit none !属性の指定 private !関数を public に設定 public Turbulence_Init public xyz_TurbScalar public xyza_TurbScalar public pyz_TurbVelX public xqz_TurbVelY public xyr_TurbVelZ public xyz_ShearKm public xyz_DispKm public xyz_DispHeat public xyz_BuoyKm public EddyViscosity !変数定義 real(DP) :: Cm = 2.0d-1 !乱流エネルギー診断式の係数 real(DP) :: MixLen = 0.0d0 !平均混合距離 !値を保存 save Cm, MixLen contains !!!------------------------------------------------------------------------!!! subroutine turbulence_init() ! ! Turbulence モジュールの初期化ルーチン ! !暗黙の型宣言禁止 implicit none real(DP):: MinDelX real(DP):: MinDelY real(DP):: MinDelZ !混合距離 MinDelX = minval(x_dx) MinDelY = minval(y_dy) MinDelZ = minval(z_dz) ! MixLen = ( MinDelX * MinDelY * MinDelZ ) ** (1.0d0/3.0d0) MixLen = min( MinDelZ, min( MinDelX, MinDelY ) ) end subroutine turbulence_init !!!------------------------------------------------------------------------!!! function xyz_TurbScalar(xyz_Var, xyz_Kh) ! ! x, y, z 方向に半格子ずれた点での乱流拡散 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP),intent(in) :: xyz_Var(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スカラー量 real(DP),intent(in) :: xyz_Kh(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !乱流拡散係数 real(DP) :: xyz_TurbScalar(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スカラー量の水平乱流拡散 ! xyz_TurbScalar = 0.0d0 !初期化 xyz_TurbScalar = & & xyz_dx_pyz( pyz_avr_xyz( xyz_Kh ) * pyz_dx_xyz( xyz_Var ) ) & & + xyz_dy_xqz( xqz_avr_xyz( xyz_Kh ) * xqz_dy_xyz( xyz_Var ) ) & & + xyz_dz_xyr( xyr_avr_xyz( xyz_Kh ) * xyr_dz_xyz( xyz_Var ) ) call StorePotTempTurb( xyz_TurbScalar ) end function xyz_TurbScalar !!!------------------------------------------------------------------------!!! function xyza_TurbScalar(xyza_Var, xyz_Kh) ! ! x, z 方向に半格子ずれた点での乱流拡散 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP),intent(in) :: xyza_Var(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !スカラー量 real(DP),intent(in) :: xyz_Kh(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !乱流拡散係数 real(DP) :: xyza_TurbScalar(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !スカラー量の水平乱流拡散 integer :: s do s = 1, SpcNum xyza_TurbScalar(:,:,:,s) = xyz_TurbScalar( xyza_Var(:,:,:,s), xyz_Kh ) end do ! call StoreMixRtTurb( xyza_TurbScalar ) end function xyza_TurbScalar !!!------------------------------------------------------------------------!!! function pyz_TurbVelX(xyz_Km, pyz_VelX, xqz_VelY, xyr_VelZ) ! ! 水平速度に対する乱流拡散 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP),intent(in) :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !水平速度 real(DP),intent(in) :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !水平速度 real(DP),intent(in) :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !鉛直速度 real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !乱流拡散係数 real(DP) :: pyz_TurbVelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スカラー量の水平乱流拡散 ! pyz_TurbVelX = 0.0d0 !初期化 pyz_TurbVelX = & & 2.0d0 * pyz_dx_xyz( xyz_Km * xyz_dx_pyz( pyz_VelX ) ) & & + pyz_dy_pqz( & & pqz_avr_xyz( xyz_Km ) * pqz_dx_xqz( xqz_VelY ) & & + pqz_avr_xyz( xyz_Km ) * pqz_dy_pyz( pyz_VelX ) & & ) & & + pyz_dz_pyr( & & pyr_avr_xyz( xyz_Km ) * pyr_dx_xyr( xyr_VelZ ) & & + pyr_avr_xyz( xyz_Km ) * pyr_dz_pyz( pyz_VelX ) & & ) & & - 2.0d0 * pyz_dx_xyz( ( xyz_Km ** 2.0d0 ) ) & & / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) ) end function pyz_TurbVelX !!!------------------------------------------------------------------------!!! function xqz_TurbVelY(xyz_Km, pyz_VelX, xqz_VelY, xyr_VelZ) ! ! 水平速度に対する乱流拡散 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP),intent(in) :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !水平速度 real(DP),intent(in) :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !水平速度 real(DP),intent(in) :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !鉛直速度 real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !乱流拡散係数 real(DP) :: xqz_TurbVelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スカラー量の水平乱流拡散 ! xqz_TurbVelY = 0.0d0 !初期化 xqz_TurbVelY = & & 2.0d0 * xqz_dy_xyz( xyz_Km * xyz_dy_xqz( xqz_VelY ) ) & & + xqz_dx_pqz( & & pqz_avr_xyz( xyz_Km ) * pqz_dy_pyz( pyz_VelX ) & & + pqz_avr_xyz( xyz_Km ) * pqz_dx_xqz( xqz_VelY ) & & ) & & + xqz_dz_xqr( & & xqr_avr_xyz( xyz_Km ) * xqr_dy_xyr( xyr_VelZ ) & & + xqr_avr_xyz( xyz_Km ) * xqr_dz_xqz( xqz_VelY ) & & ) & & - 2.0d0 * xqz_dy_xyz( ( xyz_Km ** 2.0d0 ) ) & & / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) ) end function xqz_TurbVelY !!!------------------------------------------------------------------------!!! function xyr_TurbVelZ(xyz_Km, pyz_VelX, xqz_VelY, xyr_VelZ) ! !鉛直速度に対する乱流拡散 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP),intent(in) :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !水平速度 real(DP),intent(in) :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !水平速度 real(DP),intent(in) :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !鉛直速度 real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !乱流拡散係数 real(DP) :: xyr_TurbVelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スカラー量の水平乱流拡散 ! xyr_TurbVelZ = 0.0d0 !初期化 xyr_TurbVelZ = & & + 2.0d0 * xyr_dz_xyz( xyz_Km * xyz_dz_xyr( xyr_VelZ ) ) & & + xyr_dx_pyr( & & pyr_avr_xyz( xyz_Km ) * pyr_dz_pyz( pyz_VelX ) & & + pyr_avr_xyz( xyz_Km ) * pyr_dx_xyr( xyr_VelZ ) & & ) & & + xyr_dy_xqr( & & xqr_avr_xyz( xyz_Km ) * xqr_dz_xqz( xqz_VelY ) & & + xqr_avr_xyz( xyz_Km ) * xqr_dy_xyr( xyr_VelZ ) & & ) & & - 2.0d0 * xyr_dz_xyz( ( xyz_Km ** 2.0d0 ) ) & & / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) ) end function xyr_TurbVelZ !!!------------------------------------------------------------------------!!! function xyz_ShearKm(xyz_Km, pyz_VelX, xqz_VelY, xyr_VelZ) ! !速度シアーによる乱流エネルギー生成項 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP) :: xyz_ShearKm(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !渦粘性係数の移流 real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !渦粘性係数 real(DP),intent(in) :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !水平速度 real(DP),intent(in) :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !水平速度 real(DP),intent(in) :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !鉛直速度 ! xyz_ShearKm = 0.0d0 xyz_ShearKm = & & ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) & & * ( & & ( xyz_dx_pyz( pyz_VelX ) ) ** 2.0d0 & & + ( xyz_dy_xqz( xqz_VelY ) ) ** 2.0d0 & & + ( xyz_dz_xyr( xyr_VelZ ) ) ** 2.0d0 & & + 5.0d-1 & & * ( & & ( & & xyz_avr_pyr( pyr_dz_pyz( pyz_VelX ) ) & & + xyz_avr_pyr( pyr_dx_xyr( xyr_VelZ ) ) & & ) ** 2.0d0 & & + ( & & xyz_avr_xqr( xqr_dy_xyr( xyr_VelZ ) ) & & + xyz_avr_xqr( xqr_dz_xqz( xqz_VelY ) ) & & ) ** 2.0d0 & & + ( & & xyz_avr_pqz( pqz_dx_xqz( xqz_VelY ) ) & & + xyz_avr_pqz( pqz_dy_pyz( pyz_VelX ) ) & & ) ** 2.0d0 & & ) & & ) & & - xyz_Km * ( xyz_dx_pyz( pyz_VelX ) & & + xyz_dy_xqz( xqz_VelY ) & & + xyz_dz_xyr( xyr_VelZ ) ) / 3.0d0 end function xyz_ShearKm !!!------------------------------------------------------------------------!!! function xyz_TurbKm(xyz_Km) ! ! 乱流エネルギーの拡散項 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP) :: xyz_TurbKm(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) ! 乱流エネルギーの拡散項 real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !渦粘性係数 ! xz_TurbKm = 0.0d0 xyz_TurbKm = & & xyz_TurbScalar( xyz_Km, xyz_Km ) * 5.0d-1 & & + (xyz_avr_pyz( pyz_dx_xyz( xyz_Km ) ) ) ** 2.0d0 & & + (xyz_avr_xqz( xqz_dy_xyz( xyz_Km ) ) ) ** 2.0d0 & & + (xyz_avr_xyr( xyr_dz_xyz( xyz_Km ) ) ) ** 2.0d0 end function xyz_TurbKm !!!------------------------------------------------------------------------!!! function xyz_DispKm(xyz_Km) ! !乱流エネルギーの消散項 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP) :: xyz_DispKm(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !乱流エネルギーの消散 real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !渦粘性係数 ! xz_DispKm = 0.0d0 xyz_DispKm = - (xyz_Km ** 2.0d0) * 5.0d-1 / (MixLen ** 2.0d0) end function xyz_DispKm !!!------------------------------------------------------------------------!!! function xyz_DispHeat(xyz_Km) ! !温位の散逸加熱項 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP) :: xyz_DispHeat(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !乱流エネルギーの消散 real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !渦粘性係数 ! xz_DispHeat = 0.0d0 ! xyz_DispHeat = (xyz_Km ** 3.0d0) * xyz_EffMolWtBasicZ & ! & / (xyz_ExnerBasicZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0)) xyz_DispHeat = (xyz_Km ** 3.0d0) & & / (xyz_ExnerBasicZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0)) !値の保管 call StorePotTempDisp( xyz_DispHeat ) end function xyz_DispHeat !!!------------------------------------------------------------------------!!! function xyz_BuoyKm(xyz_PotTemp) ! !乱流エネルギーの浮力項を計算. 乾燥大気版. ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP),intent(in) :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温位擾乱 real(DP) :: xyz_BuoyKm(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !浮力項 !浮力項の計算 ! xz_BuoyKm = 0.0d0 xyz_BuoyKm = & & - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) & & * xyz_avr_xyr( xyr_dz_xyz( xyz_PotTemp + xyz_PotTempBasicZ ) ) & & / ( 2.0d0 * xyz_PotTempBasicZ ) end function xyz_BuoyKm !!!------------------------------------------------------------------------!!! subroutine EddyViscosity(pyz_VelX, xqz_VelY, xyr_VelZ, xyz_PotTemp, xyz_Km, xyz_KmA) ! ! 乱流エネルギーを定常とした場合の渦拡散係数, 渦粘性係数を求める. ! この乱流パラメタリゼーションは, Mellor and Yamada (1974) の ! Level 1 Closure に対応するが, Level1 Closure に存在する ! \bar{\theta^2} は無視されている. !--- 暗黙の型宣言禁止 implicit none real(DP),intent(in) :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP),intent(in) :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP),intent(in) :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP),intent(in) :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP),intent(in) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP),intent(out) :: xyz_KmA(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) xyz_KmA = & & sqrt ( & & max( 0.0d0, & & ( xyz_ShearKm( xyz_Km, pyz_VelX, xqz_VelY, xyr_VelZ ) + & & xyz_BuoyKm( xyz_PotTemp ) ) * MixLen ** 2.0d0 )) end subroutine EddyViscosity end module Turbulence_3d