!= Module Turbulence_3d ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA, Masatsugu ! Version:: $Id: turbulence_kw1978.f90,v 1.14 2012-01-12 11:53:10 sugiyama Exp $ ! Tag Name:: $Name: arare5-20120301 $ ! 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_kw1978 ! !モデルの物理過程を計算するために必要となる関数群を束ねたモジュール !具体的には以下の項を計算するための関数を格納する. ! * 乱流拡散項 ! * 乱流エネルギーの時間発展方程式に含まれる各項 ! * 散逸加熱項 ! !モジュール読み込み use dc_types, only: DP, STRING use dc_iounit, only: FileOpen use dc_message, only: MessageNotify use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut use mpi_wrapper,only: myrank use timeset, only: DelTimeLong, TimeN use gridset, only: imin, &! x 方向の配列の下限 & imax, &! x 方向の配列の上限 & jmin, &! y 方向の配列の下限 & jmax, &! y 方向の配列の上限 & kmin, &! z 方向の配列の下限 & kmax, &! z 方向の配列の上限 & nx,ny,nz,ncmax ! use basicset, only: xyz_PTempBZ, &!基本場の温位 & xyzf_QMixBZ, &!基本場の混合比 & xyz_QMixBZ, &!基本場の混合比 & xyz_QMixBZPerMolWt, &!基本場の混合比 & xyz_EffMolWtBZ, &!基本場の混合比 & xyz_VelSoundBZ, &! 音速 & xyz_VPTempBZ, &! 仮温位 & xyz_ExnerBZ !エクスナー関数の基本場 use constants,only: Grav, & & MolWtDry, & & CpDry use composition, only: MolWtWet, & & SpcWetSymbol, & & SpcWetID, & & CondNum, &!凝結過程の数 & IdxCG, &!凝結過程(蒸気)の配列添え字 & IdxCC, &!凝結過程(雲)の配列添え字 ! & IdxNH3, &!NH3(蒸気)の配列添え字 ! & IdxH2S, &!H2S(蒸気)の配列添え字 & GasNum, &!気体の数 & IdxG !気体の配列添え字 use ChemCalc, only: xyz_LatentHeat !潜熱 ! & ReactHeatNH4SH !NH4SH の反応熱 use axesset, only: & & x_dx, &! x 方向の格子点間隔 & y_dy, &! y 方向の格子点間隔 & z_dz, &! z 方向の格子点間隔 & 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_c4_module, only: pyz_c4dx_xyz, xqz_c4dy_xyz, xyr_c4dz_xyz 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 namelist_util, only: namelist_filename use setmargin, only: SetMargin_xyz use DExnerDt, only: xyz_DExnerDt_xyz !暗黙の型宣言禁止 implicit none !属性の指定 private !関数を public に設定 public Turbulence_KW1978_Init public Turbulence_KW1978_Forcing !変数定義 real(DP), save :: Cm = 2.0d-1 !乱流エネルギー診断式の係数 real(DP), save :: MixLen = 0.0d0 !平均混合距離 real(DP), save :: KmMax = 0.0d0 !最大値 contains !!!------------------------------------------------------------------------!!! subroutine turbulence_kw1978_init ! ! Turbulence モジュールの初期化ルーチン ! !暗黙の型宣言禁止 implicit none real(DP):: MinDelX real(DP):: MinDelY real(DP):: MinDelZ integer :: l integer :: unit ! NAMELIST から情報を取得 NAMELIST /turbulence_kw1978_nml/ Cm, KmMax call FileOpen(unit, file=namelist_filename, mode='r') read(unit, NML=turbulence_kw1978_nml) close(unit) !混合距離 MinDelX = minval(x_dx) MinDelY = minval(y_dy) MinDelZ = minval(z_dz) MixLen = min( MinDelZ, min( MinDelX, MinDelY ) ) ! KmMax が設定されていない場合. ! 安定性解析では, dt / l**2 < 0.5 を満たす必要がある. ここでは 0.1 にしておいた. ! if (KmMax == 0.0d0) then KmMax = 0.1 * (MixLen ** 2.0d0) / (DelTimeLong * 2.0d0) !LeapFrog end if ! Output ! if (myrank == 0) then call MessageNotify( "M", & & "Turbulence", "Cm = %f", d=(/Cm/)) call MessageNotify( "M", & & "Turbulence", "KmMax = %f", d=(/KmMax/)) call MessageNotify( "M", & & "Turbulence", "MixLen= %f", d=(/MixLen/)) end if call HistoryAutoAddVariable( & & varname='VelXTurb',& & dims=(/'x','y','z','t'/), & & longname='Turbulence term of velocity (x)', & & units='K.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelYTurb',& & dims=(/'x','y','z','t'/), & & longname='Turbulence term of velocity (y)', & & units='K.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelZTurb',& & dims=(/'x','y','z','t'/), & & longname='Turbulence term of velocity (z)', & & units='K.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='PTempTurb',& & dims=(/'x','y','z','t'/), & & longname='Turbulence term of potential temperature', & & units='K.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='PTempDisp',& & dims=(/'x','y','z','t'/), & & longname='Dissipation term of potential temperature', & & units='K.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='ExnerDisp',& & dims=(/'x','y','z','t'/), & & longname='Dissipation term of exner function', & & units='K.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='CDensTurb',& & dims=(/'x','y','z','t'/), & & longname='Turbulence term of cloud density', & & units='kg.m-3.s-1', & & xtype='float') do l = 1, ncmax call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(l))//'_Turb', & & dims=(/'x','y','z','t'/), & & longname='Turbulence term of ' & & //trim(SpcWetSymbol(l))//' mixing ratio', & & units='kg.kg-1.s-1', & & xtype='float') end do end subroutine turbulence_kw1978_init !!!------------------------------------------------------------------------!!! subroutine turbulence_KW1978_forcing( & & pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, & & xyz_PTempBl, xyz_ExnerBl, xyzf_QMixBl, & & xyz_KmBl, xyz_KhBl, xyz_CDensBl, & & pyz_DVelXDt, xqz_DVelYDt, xyr_DVelZDt, & & xyz_DPTempDt,xyz_DExnerDt, xyzf_DQMixDt, & & xyz_DKmDt, xyz_DCDensDt, & & xyz_KmAl, xyz_KhAl & ) implicit none real(DP),allocatable,intent(in) :: pyz_VelXBl(:,:,:) !水平速度 real(DP),allocatable,intent(in) :: xqz_VelYBl(:,:,:) !水平速度 real(DP),allocatable,intent(in) :: xyr_VelZBl(:,:,:) !鉛直速度 real(DP),allocatable,intent(in) :: pyz_PTempBl(:,:,:) !温位 real(DP),allocatable,intent(in) :: xyz_ExnerBl(:,:,:) !無次元圧力 real(DP),allocatable,intent(in) :: xyzf_QMixBl(:,:,:,:) !凝縮成分の混合比 real(DP),allocatable,intent(in) :: pyz_KmBl(:,:,:) !乱流拡散係数 real(DP),allocatable,intent(in) :: pyz_KhBl(:,:,:) !乱流拡散係数 real(DP),allocatable,intent(in) :: xyz_CDensBl(:,:,:) real(DP),allocatable,intent(inout):: pyz_DVelXDt(:,:,:) !スカラー量の水平乱流拡散 real(DP),allocatable,intent(inout):: xqz_DVelYDt(:,:,:) !スカラー量の水平乱流拡散 real(DP),allocatable,intent(inout):: xyr_DVelZDt(:,:,:) !スカラー量の水平乱流拡散 real(DP),allocatable,intent(inout):: xyz_DPTempDt(:,:,:) real(DP),allocatable,intent(inout):: xyz_DExnerDt(:,:,:) real(DP),allocatable,intent(inout):: xyzf_DQMixDt(:,:,:,:) real(DP),allocatable,intent(inout):: xyz_DKmDt(:,:,:) real(DP),allocatable,intent(inout) :: xyz_DCDensDt(:,:,:) real(DP),allocatable,intent(out):: xyz_KmAl(:,:,:) !乱流拡散係数 real(DP),allocatable,intent(out):: xyz_KhAl(:,:,:) !乱流拡散係数 real(DP),allocatable :: xyz_Buoy(:,:,:) !渦粘性係数の real(DP),allocatable :: xyz_Shear(:,:,:) !渦粘性係数の real(DP),allocatable :: xyz_Disp(:,:,:) !乱流エネルギーの消散 real(DP),allocatable :: xyz_DispPI(:,:,:) !乱流エネルギーの消散 real(DP),allocatable :: xyz_DispHeat(:,:,:) !乱流エネルギーの消散 real(DP),allocatable :: xyz_Turb(:,:,:) ! real(DP),allocatable :: xyzf_Turb(:,:,:,:) !スカラー量の水平乱流拡散 real(DP),allocatable :: pyz_Turb(:,:,:) real(DP),allocatable :: xqz_Turb(:,:,:) real(DP),allocatable :: xyr_Turb(:,:,:) real(DP),allocatable :: pyz_DVelXDt0(:,:,:) !スカラー量の水平乱流拡散 real(DP),allocatable :: xqz_DVelYDt0(:,:,:) !スカラー量の水平乱流拡散 real(DP),allocatable :: xyr_DVelZDt0(:,:,:) !スカラー量の水平乱流拡散 real(DP),allocatable :: xyz_DPTempDt0(:,:,:) real(DP),allocatable :: xyz_DExnerDt0(:,:,:) real(DP),allocatable :: xyzf_DQMixDt0(:,:,:,:) real(DP),allocatable :: xyz_DKmDt0(:,:,:) real(DP),allocatable :: xyz_DCDensDt0(:,:,:) real(DP),allocatable :: pyz_PTempBlAll(:,:,:) real(DP),allocatable :: xyzf_QMixBlAll(:,:,:,:) integer :: f !---------------------------------- ! 初期化 ! allocate( & & pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax), & & xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax), & & xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax), & & pyz_PTempBl(imin:imax,jmin:jmax,kmin:kmax), & & xyz_ExnerBl(imin:imax,jmin:jmax,kmin:kmax), & & xyzf_QMixBl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax), & & xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax), & & pyz_KhBl(imin:imax,jmin:jmax,kmin:kmax), & & xyz_CDensBl(imin:imax,jmin:jmax,kmin:kmax), & & pyz_DVelXDt(imin:imax,jmin:jmax,kmin:kmax), & & xqz_DVelYDt(imin:imax,jmin:jmax,kmin:kmax), & & xyr_DVelZDt(imin:imax,jmin:jmax,kmin:kmax), & & xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax), & & xyz_DExnerDt(imin:imax,jmin:jmax,kmin:kmax), & & xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax), & & xyz_DKmDt(imin:imax,jmin:jmax,kmin:kmax), & & xyz_DCDensDt(imin:imax,jmin:jmax,kmin:kmax), & & xyz_KmAl(imin:imax,jmin:jmax,kmin:kmax), & & xyz_KhAl(imin:imax,jmin:jmax,kmin:kmax), & & xyz_Buoy(imin:imax,jmin:jmax,kmin:kmax), & & xyz_Shear(imin:imax,jmin:jmax,kmin:kmax), & & xyz_Disp(imin:imax,jmin:jmax,kmin:kmax), & & xyz_DispPI(imin:imax,jmin:jmax,kmin:kmax), & & xyz_DispHeat(imin:imax,jmin:jmax,kmin:kmax), & & xyz_Turb(imin:imax,jmin:jmax,kmin:kmax), & & xyzf_Turb(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax), & & pyz_Turb(imin:imax,jmin:jmax,kmin:kmax), & & xqz_Turb(imin:imax,jmin:jmax,kmin:kmax), & & xyr_Turb(imin:imax,jmin:jmax,kmin:kmax), & & pyz_DVelXDt0(imin:imax,jmin:jmax,kmin:kmax), & & xqz_DVelYDt0(imin:imax,jmin:jmax,kmin:kmax), & & xyr_DVelZDt0(imin:imax,jmin:jmax,kmin:kmax), & & xyz_DPTempDt0(imin:imax,jmin:jmax,kmin:kmax), & & xyz_DExnerDt0(imin:imax,jmin:jmax,kmin:kmax), & & xyzf_DQMixDt0(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax), & & xyz_DKmDt0(imin:imax,jmin:jmax,kmin:kmax), & & xyz_DCDensDt0(imin:imax,jmin:jmax,kmin:kmax), & & pyz_PTempBlAll(imin:imax,jmin:jmax,kmin:kmax), & & xyzf_QMixBlAll(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) & & ) pyz_VelXBl = 0.0d0 xqz_VelYBl = 0.0d0 xyr_VelZBl = 0.0d0 pyz_PTempBl = 0.0d0 xyz_ExnerBl = 0.0d0 xyzf_QMixBl = 0.0d0 xyz_KmBl = 0.0d0 xyz_KhBl = 0.0d0 xyz_CDensBl = 0.0d0 pyz_DVelXDt = 0.0d0 xqz_DVelYDt = 0.0d0 xyr_DVelZDt = 0.0d0 xyz_DPTempDt = 0.0d0 xyz_DExnerDt = 0.0d0 xyzf_DQMixDt = 0.0d0 xyz_DKmDt = 0.0d0 xyz_DCDensDt = 0.0d0 xyz_KmAl = 0.0d0 xyz_KhAl = 0.0d0 xyz_Buoy = 0.0d0 xyz_Shear = 0.0d0 xyz_Disp = 0.0d0 xyz_DispPI = 0.0d0 xyz_DispHeat = 0.0d0 xyz_Turb = 0.0d0 xyzf_Turb = 0.0d0 pyz_Turb = 0.0d0 xqz_Turb = 0.0d0 xyr_Turb = 0.0d0 pyz_DVelXDt0 = 0.0d0 xqz_DVelYDt0 = 0.0d0 xyr_DVelZDt0 = 0.0d0 xyz_DPTempDt0 = 0.0d0 xyz_DExnerDt0 = 0.0d0 xyzf_DQMixDt0 = 0.0d0 xyz_DKmDt0 = 0.0d0 xyz_DCDensDt0 = 0.0d0 pyz_PTempBlAll = 0.0d0 xyzf_QMixBlAll = 0.0d0 pyz_PTempBlAll = pyz_PTempBl ! xyz_PTempBlAll = xyz_PTempBl + xyz_PTempBZ xyzf_QMixBlAll = xyzf_QMixBl + xyzf_QMixBZ pyz_DVelXDt0 = pyz_DVelXDt xqz_DVelYDt0 = xqz_DVelYDt xyr_DVelZDt0 = xyr_DVelZDt xyz_DPTempDt0 = xyz_DPTempDt xyz_DExnerDt0 = xyz_DExnerDt xyzf_DQMixDt0 = xyzf_DQMixDt xyz_DKmDt0 = xyz_DKmDt xyz_DCDensDt0 = xyz_DCDensDt !---------------------------------- ! 拡散係数の時間発展 (エネルギー方程式を Km の式に変形したもの) ! ! Buoyancy term ! xyz_Buoy = xyz_BuoyMoistKm(xyz_PTempBl, xyz_ExnerBl, xyzf_QMixBl) ! xyz_Buoy = & ! & - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) & ! & * xyz_avr_xyr( xyr_dz_xyz( xyz_PTempBlAll ) ) & ! & / ( 2.0d0 * xyz_PTempBZ ) xyz_Shear = & & ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) & & * ( & & ( xyz_dx_pyz( pyz_VelXBl ) ) ** 2.0d0 & & + ( xyz_dy_xqz( xqz_VelYBl ) ) ** 2.0d0 & & + ( xyz_dz_xyr( xyr_VelZBl ) ) ** 2.0d0 & & + 5.0d-1 & & * ( & & ( & & xyz_avr_pyr( pyr_dz_pyz( pyz_VelXBl ) ) & & + xyz_avr_pyr( pyr_dx_xyr( xyr_VelZBl ) ) & & ) ** 2.0d0 & & + ( & & xyz_avr_xqr( xqr_dy_xyr( xyr_VelZBl ) ) & & + xyz_avr_xqr( xqr_dz_xqz( xqz_VelYBl ) ) & & ) ** 2.0d0 & & + ( & & xyz_avr_pqz( pqz_dx_xqz( xqz_VelYBl ) ) & & + xyz_avr_pqz( pqz_dy_pyz( pyz_VelXBl ) ) & & ) ** 2.0d0 & & ) & & ) & & - xyz_KmBl * ( xyz_dx_pyz( pyz_VelXBl ) & & + xyz_dy_xqz( xqz_VelYBl ) & & + xyz_dz_xyr( xyr_VelZBl ) ) / 3.0d0 ! t - \Delta t で評価 ! xyz_Disp = - (xyz_KmBl ** 2.0d0) * 5.0d-1 / (MixLen ** 2.0d0) ! tendency ! xyz_DKmDt = xyz_DKmDt0 + xyz_Buoy + xyz_Shear + xyz_Disp ! 時間積分 ! ! xyz_KmAl = xyz_KmBl + (2.0d0 * DelTimeLong) * xyz_DKmDt xyz_KmAl = 155.0d0 ! xyz_KmAl = 0.0d0 ! 上限値の設定 ! xyz_KmAl = max( 0.0d0, min( xyz_KmAl, KmMax ) ) ! Kh の設定 ! ! xyz_KhAl = 3.0d0 * xyz_KmAl xyz_KhAl = xyz_KmAl !-------------------------------- ! 雲密度の tendency ! xyz_Turb = & & xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyz_CDensBl ) ) & & + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyz_CDensBl ) ) & & + xyz_dz_xyr( xyr_avr_xyz( xyz_KhBl ) * xyr_dz_xyz( xyz_CDensBl ) ) xyz_DCDensDt = xyz_DCDensDt0 + xyz_Turb ! output ! call HistoryAutoPut(TimeN, 'CDensTurb', xyz_Turb(1:nx,1:ny,1:nz)) !-------------------------------- ! 温位の tendency ! xyz_Turb = & & xyz_dx_pyz( pyz_KhBl * pyz_PTempBlAll ) ! & + xyz_dy_xqz( xqz_KhBl * xqz_PTempBlAll ) & ! & + xyz_dz_xyr( xyr_KhBl * xyr_PTempBlAll ) ! xyz_Turb = & ! & xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyz_PTempBlAll ) ) & ! & + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyz_PTempBlAll ) ) & ! & + xyz_dz_xyr( xyr_avr_xyz( xyz_KhBl ) * xyr_dz_xyz( xyz_PTempBlAll ) ) ! 散逸加熱項 ! xyz_DispHeat = (xyz_KmBl ** 3.0d0) & ! & / (xyz_ExnerBZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0)) xyz_DispHeat = 0.0d0 xyz_DPTempDt = xyz_DPTempDt0 + xyz_Turb + xyz_DispHeat call HistoryAutoPut(TimeN, 'PTempDisp', xyz_DispHeat(1:nx, 1:ny, 1:nz)) call HistoryAutoPut(TimeN, 'PTempTurb', xyz_Turb(1:nx, 1:ny, 1:nz)) !-------------------------------- ! 混合比の tendency ! do f = 1, ncmax xyzf_Turb(:,:,:,f) = & & xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyzf_QMixBlAll(:,:,:,f) ) ) & & + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyzf_QMixBlAll(:,:,:,f) ) ) & & + xyz_dz_xyr( xyr_avr_xyz( xyz_KhBl ) * xyr_dz_xyz( xyzf_QMixBlAll(:,:,:,f) ) ) call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_Turb', xyzf_Turb(1:nx,1:ny,1:nz,f)) end do xyzf_DQMixDt = xyzf_DQMixDt0 + xyzf_Turb !-------------------------------- ! 各速度成分の tendency ! pyz_Turb = & & 2.0d0 * pyz_dx_xyz( xyz_KmBl * xyz_dx_pyz( pyz_VelXBl ) )& & + pyz_dy_pqz( & & pqz_avr_xyz( xyz_KmBl ) * pqz_dx_xqz( xqz_VelYBl ) & & + pqz_avr_xyz( xyz_KmBl ) * pqz_dy_pyz( pyz_VelXBl ) & & ) & & + pyz_dz_pyr( & & pyr_avr_xyz( xyz_KmBl ) * pyr_dx_xyr( xyr_VelZBl ) & & + pyr_avr_xyz( xyz_KmBl ) * pyr_dz_pyz( pyz_VelXBl ) & & ) & & - 2.0d0 * pyz_dx_xyz( ( xyz_KmBl ** 2.0d0 ) ) & & / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) ) pyz_DVelXDt = pyz_DVelXDt0 + pyz_Turb call HistoryAutoPut(TimeN, 'VelXTurb', pyz_Turb(1:nx, 1:ny, 1:nz)) xqz_Turb = & & 2.0d0 * xqz_dy_xyz( xyz_KmBl * xyz_dy_xqz( xqz_VelYBl ) ) & & + xqz_dx_pqz( & & pqz_avr_xyz( xyz_KmBl ) * pqz_dy_pyz( pyz_VelXBl ) & & + pqz_avr_xyz( xyz_KmBl ) * pqz_dx_xqz( xqz_VelYBl ) & & ) & & + xqz_dz_xqr( & & xqr_avr_xyz( xyz_KmBl ) * xqr_dy_xyr( xyr_VelZBl ) & & + xqr_avr_xyz( xyz_KmBl ) * xqr_dz_xqz( xqz_VelYBl ) & & ) & & - 2.0d0 * xqz_dy_xyz( ( xyz_KmBl ** 2.0d0 ) ) & & / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) ) xqz_DVelYDt = xqz_DVelYDt0 + xqz_Turb call HistoryAutoPut(TimeN, 'VelYTurb', xqz_Turb(1:nx, 1:ny, 1:nz)) xyr_Turb = & & + 2.0d0 * xyr_dz_xyz( xyz_KmBl * xyz_dz_xyr( xyr_VelZBl ) ) & & + xyr_dx_pyr( & & pyr_avr_xyz( xyz_KmBl ) * pyr_dz_pyz( pyz_VelXBl ) & & + pyr_avr_xyz( xyz_KmBl ) * pyr_dx_xyr( xyr_VelZBl ) & & ) & & + xyr_dy_xqr( & & xqr_avr_xyz( xyz_KmBl ) * xqr_dz_xqz( xqz_VelYBl ) & & + xqr_avr_xyz( xyz_KmBl ) * xqr_dy_xyr( xyr_VelZBl ) & & ) & & - 2.0d0 * xyr_dz_xyz( ( xyz_KmBl ** 2.0d0 ) ) & & / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) ) xyr_DVelZDt = xyr_DVelZDt0 + xyr_Turb call HistoryAutoPut(TimeN, 'VelZTurb', xyr_Turb(1:nx, 1:ny, 1:nz)) !-------------------- ! Exner function ! xyz_DispPI = xyz_DExnerDt_xyz( xyz_DispHeat ) xyz_DExnerDt = xyz_DExnerDt0 + xyz_DispPI call HistoryAutoPut(TimeN, 'ExnerDisp', xyz_DispPI(1:nx, 1:ny, 1:nz)) ! Set Margin ! call SetMargin_xyz(xyz_KmAl) call SetMargin_xyz(xyz_KhAl) ! deallocate( & ! & pyz_VelXBl, & ! & xqz_VelYBl, & ! & xyr_VelZBl, & ! & xyz_PTempBl, & ! & xyz_ExnerBl, & ! & xyzf_QMixBl, & ! & xyz_KmBl, & ! & xyz_KhBl, & ! & xyz_CDensBl, & ! & pyz_DVelXDt, & ! & xqz_DVelYDt, & ! & xyr_DVelZDt, & ! & xyz_DPTempDt, & ! & xyz_DExnerDt, & ! & xyzf_DQMixDt, & ! & xyz_DKmDt, & ! & xyz_DCDensDt, & ! & xyz_KmAl, & ! & xyz_KhAl, & ! & xyz_Buoy, & ! & xyz_Shear, & ! & xyz_Disp, & ! & xyz_DispPI, & ! & xyz_DispHeat, & ! & xyz_Turb, & ! & xyzf_Turb, & ! & pyz_Turb, & ! & xqz_Turb, & ! & xyr_Turb, & ! & pyz_DVelXDt0, & ! & xqz_DVelYDt0, & ! & xyr_DVelZDt0, & ! & xyz_DPTempDt0, & ! & xyz_DExnerDt0, & ! & xyzf_DQMixDt0, & ! & xyz_DKmDt0, & ! & xyz_DCDensDt0, & ! & xyz_PTempBlAll, & ! & xyzf_QMixBlAll & ! & ) end subroutine Turbulence_KW1978_forcing !!!------------------------------------------------------------------------!!! function xyz_BuoyMoistKm(xyz_PTemp, xyz_Exner, xyzf_QMix) ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax) !温位 real(DP), intent(in) :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax) !無次元圧力 real(DP), intent(in) :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax) !凝縮成分の混合比 real(DP) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax) !凝縮成分の混合比 real(DP) :: xyzf_QMixAll2(imin:imax,jmin:jmax,kmin:kmax, ncmax) !凝縮成分の混合比 real(DP) :: xyz_BuoyMoistKm(imin:imax,jmin:jmax,kmin:kmax) ! real(DP) :: xyzf_LatentHeat(imin:imax,jmin:jmax,kmin:kmax, GasNum) !潜熱 real(DP) :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax) !温度 real(DP) :: xyz_EffHeat(imin:imax,jmin:jmax,kmin:kmax) ! real(DP) :: xyz_EffPTemp(imin:imax,jmin:jmax,kmin:kmax) !浮力に対する温度差の寄与 real(DP) :: xyz_EffMolWt(imin:imax,jmin:jmax,kmin:kmax) !浮力に対する分子量差の寄与 real(DP) :: xyzf_QMixPerMolWt(imin:imax,jmin:jmax,kmin:kmax, GasNum) !混合比/分子量 integer :: s !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xyz_TempAll = ( xyz_PTemp + xyz_PTempBZ ) * ( xyz_Exner + xyz_ExnerBZ ) xyzf_QMixAll = xyzf_QMixBZ + xyzf_QMix xyzf_LatentHeat = 0.0d0 !作業配列の初期化. 気体のみ利用 do s = 1, GasNum xyzf_QMixPerMolWt(:,:,:,s) = xyzf_QMix(:,:,:,IdxG(s)) / MolWtWet(IdxG(s)) end do !温度の効果 xyz_EffPTemp = xyz_PTemp / xyz_PTempBZ !分子量効果 + 引きづりの効果 xyz_EffMolWt = & & + sum(xyzf_QMixPerMolWt, 4) & & / ( 1.0d0 / MolWtDry + xyz_QMixBZPerMolWt ) & & - sum(xyzf_QMix, 4) / ( 1.0d0 + xyz_QMixBZ ) !蒸気が蒸発する場合の潜熱を計算 ! 分子量の部分はいつでも効くが潜熱は飽和していないと効かないので, ! 雲の混合比がゼロの時には, 潜熱の寄与はゼロとなるように調節している xyzf_QMixAll2 = xyzf_QMixAll ! xzya_QMixAll2(:,:,:,IdxNH3) = & ! & xyzf_QMixAll(:,:,:,IdxNH3) - xyzf_QMixAll(:,:,:,IdxH2S) do s = 1, CondNum xyzf_LatentHeat(:,:,:,s) = & & xyz_LatentHeat( SpcWetID(IdxCC(s)), xyz_TempAll ) & & * xyzf_QMixAll2(:,:,:,IdxCG(s)) & & * ( 5.0d-1 + sign( 5.0d-1, (xyzf_QMixAll2(:,:,:,IdxCC(s)) - 1.0d-4) ) ) end do xyz_EffHeat = ( sum( xyzf_LatentHeat, 4 ) * xyz_EffMolWtBZ & ! & + ReactHeatNH4SH * & ! & (xyzf_QMixAll(:,:,:,IdxH2S) + xyzf_QMixAll(:,:,:,IdxH2S)) & & ) / ( CpDry * xyz_ExnerBZ ) !乱流拡散係数の時間発展式の浮力項を決める xyz_BuoyMoistKm = & & - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) & & * xyz_avr_xyr( & & xyr_dz_xyz( & & xyz_EffHeat & & + xyz_PTempBZ / xyz_EffMolWtBZ & & * ( 1.0d0 + xyz_EffPTemp + xyz_EffMolWt ) & & ) & & ) & & / ( 2.0d0 * xyz_PTempBZ / xyz_EffMolWtBZ) end function xyz_BuoyMoistKm end module Turbulence_kw1978