!= Module Turbulence_kw1978 ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA, Masatsugu ! Version:: $Id: turbulence_kw1978_v2.f90,v 1.4 2015/02/19 02:17:24 sugiyama Exp $ ! Tag Name:: $Name: $ ! Copyright:: Copyright (C) GFD Dennou Club, 2014. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] module Turbulence_kw1978 ! ! Klemp and Wilhelmson (1978) の乱流過程 ! 乱流エネルギーの時間発展方程式を解くことで乱流拡散係数を決める ! !モジュール読み込み ! use dc_types, only: DP, STRING !暗黙の型宣言禁止 ! 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 !乱流拡散係数の最大値 real(DP), save :: Cm_Cm_MixLen_MixLen = 0.0d0 logical, save :: FlagDExnerDtTurb =.true. !圧力方程式に乱流拡散項を考慮するかのスイッチ integer, save :: IDTurbulence = 0 ! 短い時間ステップの計算方法 integer, parameter :: IDTurbulence_std = 1 ! 1: 微分平均モジュールの利用. integer, parameter :: IDTurbulence_2D = 2 ! 2: 2D 版計算ルーチン integer, parameter :: IDTurbulence_3D = 3 ! 3: 3D 版計算ルーチン character(STRING), parameter:: module_name = 'turbulence_kw1978_v3' contains !!!------------------------------------------------------------------------!!! subroutine turbulence_kw1978_init ! ! Turbulence モジュールの初期化ルーチン ! !モジュール読み込み ! use mpi_wrapper, only : myrank use dc_types, only : STRING use axesset, only : dx, &! x 方向の格子点間隔 & dy, &! y 方向の格子点間隔 & dz ! z 方向の格子点間隔 use dc_iounit, only : FileOpen use dc_message, only : MessageNotify use timeset, only : DelTimeLong use gridset, only : FlagCalc3D use namelist_util, only : namelist_filename !暗黙の型宣言禁止 ! implicit none ! 作業変数 ! integer :: unit character(STRING) :: FlagTurbulence = "" !乱流計算方法の選択 !------------------------------------------------------------------- ! NAMELIST から情報を取得 ! NAMELIST /turbulence_kw1978_nml/ & & Cm, KmMax, & & FlagDExnerDtTurb, FlagTurbulence call FileOpen(unit, file=namelist_filename, mode='r') read(unit, NML=turbulence_kw1978_nml) close(unit) !------------------------------------------------------------------- ! 混合距離 ! 2 次元計算の場合には DelY に依存しないようにするために if 文を利用. ! if ( FlagCalc3D ) then MixLen = (dx * dy * dz ) ** (1.0d0 / 3.0d0) else MixLen = sqrt( dx * dz ) end if ! まとめておく ! Cm_Cm_MixLen_MixLen = Cm * Cm * MixLen * MixLen !------------------------------------------------------------------- ! 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 !------------------------------------------------------------------- ! 計算方法の選択 ! std: 従来の微分平均演算を利用 ! 3d, 2d: 微分平均演算を手動でインライン展開した版 ! if ( FlagTurbulence == "std" ) then IDTurbulence = IDTurbulence_std else if ( FlagCalc3D ) then IDTurbulence = IDTurbulence_3d else IDTurbulence = IDTurbulence_2d end if end if !------------------------------------------------------------------- ! tendency の出力 ! call turbulence_kw1978_output !------------------------------------------------------------------- ! Output ! if (myrank == 0) then call MessageNotify( "M", module_name, "Cm = %f", d=(/Cm/)) call MessageNotify( "M", module_name, "KmMax = %f", d=(/KmMax/)) call MessageNotify( "M", module_name, "MixLen= %f", d=(/MixLen/)) call MessageNotify( "M", module_name, & & "FlagDExnerDtTurb = %b", l=(/ FlagDExnerDtTurb /)) call MessageNotify( "M", module_name, & & "IDTurbulence = %d", i=(/ IDTurbulence /)) end if 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, & & pyz_DVelXDt, xqz_DVelYDt, xyr_DVelZDt, & & xyz_DPTempDt,xyz_DExnerDt, xyzf_DQMixDt, & & xyz_DKmDt, & & xyz_KmAl, xyz_KhAl & ) use dc_types, only : DP use gtool_historyauto, only : & & HistoryAutoPut use timeset, only : DelTimeLong, & & TimeN use gridset, only : imin, &! x 方向の配列の下限 & imax, &! x 方向の配列の上限 & jmin, &! y 方向の配列の下限 & jmax, &! y 方向の配列の上限 & kmin, &! z 方向の配列の下限 & kmax, &! z 方向の配列の上限 & ncmax, & & nx,ny,nz use basicset, only : xyz_PTempBZ, &!基本場の温位 & xyz_ExnerBZ, &!基本場のエクスナー関数 & xyzf_QMixBZ !基本場の混合比 use setmargin, only : SetMargin_xyz use DExnerDt, only : xyz_DExnerDt_xyz use composition,only: SpcWetSymbol ! 暗黙の型宣言禁止 ! implicit none ! 入出力変数 ! real(DP),intent(in) :: pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax) !水平速度 real(DP),intent(in) :: xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax) !水平速度 real(DP),intent(in) :: xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax) !鉛直速度 real(DP),intent(in) :: xyz_PTempBl(imin:imax,jmin:jmax,kmin:kmax) !温位 real(DP),intent(in) :: xyz_ExnerBl(imin:imax,jmin:jmax,kmin:kmax) !無次元圧力 real(DP),intent(in) :: xyzf_QMixBl(imin:imax,jmin:jmax,kmin:kmax, ncmax) !凝縮成分の混合比 real(DP),intent(in) :: xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax) !乱流拡散係数 real(DP),intent(in) :: xyz_KhBl(imin:imax,jmin:jmax,kmin:kmax) !乱流拡散係数 real(DP),intent(inout):: pyz_DVelXDt(imin:imax,jmin:jmax,kmin:kmax) !スカラー量の水平乱流拡散 real(DP),intent(inout):: xqz_DVelYDt(imin:imax,jmin:jmax,kmin:kmax) !スカラー量の水平乱流拡散 real(DP),intent(inout):: xyr_DVelZDt(imin:imax,jmin:jmax,kmin:kmax) !スカラー量の水平乱流拡散 real(DP),intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax) real(DP),intent(inout):: xyz_DExnerDt(imin:imax,jmin:jmax,kmin:kmax) real(DP),intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax) real(DP),intent(inout):: xyz_DKmDt(imin:imax,jmin:jmax,kmin:kmax) real(DP),intent(out):: xyz_KmAl(imin:imax,jmin:jmax,kmin:kmax) !乱流拡散係数 real(DP),intent(out):: xyz_KhAl(imin:imax,jmin:jmax,kmin:kmax) !乱流拡散係数 ! 作業変数 ! real(DP) :: xyz_Buoy(imin:imax,jmin:jmax,kmin:kmax) !渦粘性係数の real(DP) :: xyz_BuoyT(imin:imax,jmin:jmax,kmin:kmax) !渦粘性係数の real(DP) :: xyz_BuoyM(imin:imax,jmin:jmax,kmin:kmax) !渦粘性係数の real(DP) :: xyz_Shear(imin:imax,jmin:jmax,kmin:kmax) !渦粘性係数の real(DP) :: xyz_Diff(imin:imax,jmin:jmax,kmin:kmax) !渦粘性係数の real(DP) :: xyz_Disp(imin:imax,jmin:jmax,kmin:kmax) !乱流エネルギーの消散 real(DP) :: xyz_DispPI(imin:imax,jmin:jmax,kmin:kmax) !乱流エネルギーの消散 real(DP) :: xyz_DispHeat(imin:imax,jmin:jmax,kmin:kmax) !乱流エネルギーの消散 real(DP) :: xyz_Turb(imin:imax,jmin:jmax,kmin:kmax) ! real(DP) :: xyzf_Turb(imin:imax,jmin:jmax,kmin:kmax, ncmax) !スカラー量の水平乱流拡散 real(DP) :: pyz_Turb(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_Turb(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_Turb(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_TempBlAll(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_PTempBlAll(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyzf_QMixBlAll(imin:imax,jmin:jmax,kmin:kmax, ncmax) integer :: f !---------------------------------- ! 初期化 ! xyz_PTempBlAll = xyz_PTempBl + xyz_PTempBZ xyzf_QMixBlAll = xyzf_QMixBl + xyzf_QMixBZ xyz_TempBlAll = xyz_PTempBlAll * ( xyz_ExnerBl + xyz_ExnerBZ ) !---------------------------------- ! 各予報変数に対して乱流拡散による tendency を計算する. ! 内部サブルーチンを利用. ! select case ( IDTurbulence ) case ( IDTurbulence_std ) call turbulence_kw1978_std case ( IDTurbulence_3d ) call turbulence_kw1978_center2_3d case ( IDTurbulence_2d ) call turbulence_kw1978_center2_2d end select !----------------------------------------- ! Km の時間積分 ! ! tendency ! xyz_DKmDt = xyz_DKmDt + xyz_Buoy + xyz_Shear + xyz_Diff + xyz_Disp ! tendency の保管 ! call HistoryAutoPut(TimeN, 'KmBuoyT', xyz_BuoyT(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'KmBuoyM', xyz_BuoyM(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'KmShear', xyz_Shear(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'KmDiff', xyz_Diff(1:nx,1:ny,1:nz) ) call HistoryAutoPut(TimeN, 'KmDisp', xyz_Disp(1:nx,1:ny,1:nz) ) ! 時間積分 ! xyz_KmAl = xyz_KmBl + (2.0d0 * DelTimeLong) * xyz_DKmDt ! 上限値の設定 ! xyz_KmAl = max( 0.0d0, min( xyz_KmAl, KmMax ) ) ! Set Margin ! call SetMargin_xyz(xyz_KmAl) ! Kh の設定 ! xyz_KhAl = 3.0d0 * xyz_KmAl !-------------------------------- ! 温位の tendency ! xyz_DPTempDt = xyz_DPTempDt + 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 ! xyzf_DQMixDt = xyzf_DQMixDt + xyzf_Turb do f = 1, ncmax call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_Turb', xyzf_Turb(1:nx,1:ny,1:nz,f)) end do !-------------------------------- ! 各速度成分の tendency ! pyz_DVelXDt = pyz_DVelXDt + pyz_Turb xqz_DVelYDt = xqz_DVelYDt + xqz_Turb xyr_DVelZDt = xyr_DVelZDt + xyr_Turb call HistoryAutoPut(TimeN, 'VelXTurb', pyz_Turb(1:nx, 1:ny, 1:nz)) call HistoryAutoPut(TimeN, 'VelYTurb', xqz_Turb(1:nx, 1:ny, 1:nz)) call HistoryAutoPut(TimeN, 'VelZTurb', xyr_Turb(1:nx, 1:ny, 1:nz)) !-------------------- ! tendency of Exner function ! if ( FlagDExnerDtTurb ) then xyz_DispPI = xyz_DExnerDt_xyz( xyz_DispHeat ) else xyz_DispPi = 0.0d0 end if xyz_DExnerDt = xyz_DExnerDt + xyz_DispPI call HistoryAutoPut(TimeN, 'ExnerDisp', xyz_DispPI(1:nx, 1:ny, 1:nz)) contains subroutine turbulence_kw1978_std !モジュールの読み込み ! use axesset, only : xyz_pyz, xyz_xqz, xyz_xyr, & & xyz_pqz, xyz_pyr, xyz_xqr, & & pqz_xyz, pyz_xyz, & & xqz_xyz, xqr_xyz, & & xyr_xyz, pyr_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, & & xyr_dx_pyr, xyr_dy_xqr, xqr_dz_xqz, & & xqz_dx_pqz, xqr_dy_xyr, xqz_dz_xqr use basicset, only : xyz_ExnerBZ, &!基本場のエクスナー関数 & xyz_DensBZ, &!基本場の密度 & pyz_DensBZ, &!基本場の密度 & xqz_DensBZ, &!基本場の密度 & xyr_DensBZ !基本場の密度 use constants, only : CpDry ! 暗黙の型宣言禁止 ! implicit none ! 作業変数 ! real(DP) :: xyz_KmBlKmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DensBZKhBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DensBZKmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DensBZKmBlKmBl(imin:imax, jmin:jmax, kmin:kmax) !--------------------------------------------------------- ! 作業変数の定義 ! xyz_KmBlKmBl = xyz_KmBl * xyz_KmBl xyz_DensBZKhBl = xyz_DensBZ * xyz_KhBl xyz_DensBZKmBl = xyz_DensBZ * xyz_KmBl xyz_DensBZKmBlKmBl = xyz_DensBZ * xyz_KmBl * xyz_KmBl !--------------------------------------------------------- ! Km の時間発展方程式の各項 ! 散逸項 ! xyz_Disp = - xyz_KmBlKmBl * 5.0d-1 / (MixLen ** 2.0d0) ! 浮力項 ! call turbulence_kw1978_BuoyancyKm ! シアー項 ! xyz_Shear = & & Cm_Cm_MixLen_MixLen & & * ( & & ( 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_pyr( pyr_dz_pyz( pyz_VelXBl ) ) & & + xyz_pyr( pyr_dx_xyr( xyr_VelZBl ) ) & & ) ** 2.0d0 & & + ( & & xyz_xqr( xqr_dy_xyr( xyr_VelZBl ) ) & & + xyz_xqr( xqr_dz_xqz( xqz_VelYBl ) ) & & ) ** 2.0d0 & & + ( & & xyz_pqz( pqz_dx_xqz( xqz_VelYBl ) ) & & + xyz_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 ! 拡散項 ! xyz_Diff = & & 5.0d-1 & & * ( & & xyz_dx_pyz(pyz_dx_xyz(xyz_KmBlKmBl)) & & + xyz_dy_xqz(xqz_dy_xyz(xyz_KmBlKmBl)) & & + xyz_dz_xyr(xyr_dz_xyz(xyz_KmBlKmBl)) & & ) & & + ( & & (xyz_pyz(pyz_dx_xyz(xyz_KmBl))) ** 2.0d0 & & + (xyz_xqz(xqz_dy_xyz(xyz_KmBl))) ** 2.0d0 & & + (xyz_xyr(xyr_dz_xyz(xyz_KmBl))) ** 2.0d0 & & ) !-------------------------------- ! 温位の tendency ! xyz_Turb = & & xyz_dx_pyz( pyz_xyz( xyz_KhBl ) * pyz_dx_xyz( xyz_PTempBlAll ) ) & & + xyz_dy_xqz( xqz_xyz( xyz_KhBl ) * xqz_dy_xyz( xyz_PTempBlAll ) ) & & + xyz_dz_xyr( xyr_xyz( xyz_DensBZKhBl ) & & * xyr_dz_xyz( xyz_PTempBlAll ) ) / xyz_DensBZ xyz_DispHeat = (xyz_KmBl ** 3.0d0) & & / (xyz_ExnerBZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0)) !-------------------------------- ! 混合比の tendency ! do f = 1, ncmax xyzf_Turb(:,:,:,f) = & & xyz_dx_pyz( pyz_xyz( xyz_KhBl ) & & * pyz_dx_xyz( xyzf_QMixBlAll(:,:,:,f) ) ) & & + xyz_dy_xqz( xqz_xyz( xyz_KhBl ) & & * xqz_dy_xyz( xyzf_QMixBlAll(:,:,:,f) ) ) & & + xyz_dz_xyr( xyr_xyz( xyz_DensBZKhBl ) & & * xyr_dz_xyz( xyzf_QMixBlAll(:,:,:,f) ) ) & & / xyz_DensBZ end do !-------------------------------- ! Turb.u ! pyz_Turb = & & 2.0d0 * pyz_dx_xyz( xyz_KmBl * xyz_dx_pyz( pyz_VelXBl ) ) & & + pyz_dy_pqz( & & pqz_xyz( xyz_KmBl ) * pqz_dx_xqz( xqz_VelYBl ) & & + pqz_xyz( xyz_KmBl ) * pqz_dy_pyz( pyz_VelXBl ) & & ) & & + pyz_dz_pyr( & & pyr_xyz( xyz_DensBZKmBl ) * pyr_dx_xyr( xyr_VelZBl ) & & + pyr_xyz( xyz_DensBZKmBl ) * pyr_dz_pyz( pyz_VelXBl ) & & ) / pyz_DensBZ & & - 2.0d0 * pyz_dx_xyz( xyz_KmBlKmBl ) & & / ( 3.0d0 * Cm_Cm_MixLen_MixLen ) !----------------------------------------------- ! Turb.v ! xqz_Turb = & & 2.0d0 * xqz_dy_xyz( xyz_KmBl * xyz_dy_xqz( xqz_VelYBl ) ) & & + xqz_dx_pqz( & & pqz_xyz( xyz_KmBl ) * pqz_dy_pyz( pyz_VelXBl ) & & + pqz_xyz( xyz_KmBl ) * pqz_dx_xqz( xqz_VelYBl ) & & ) & & + xqz_dz_xqr( & & xqr_xyz( xyz_DensBZKmBl ) * xqr_dy_xyr( xyr_VelZBl ) & & + xqr_xyz( xyz_DensBZKmBl ) * xqr_dz_xqz( xqz_VelYBl ) & & ) / xqz_DensBZ & & - 2.0d0 * xqz_dy_xyz( xyz_KmBlKmBl ) & & / ( 3.0d0 * Cm_Cm_MixLen_MixLen ) !----------------------------------------------- ! Turb.w ! xyr_Turb = & & + 2.0d0 * xyr_dz_xyz( xyz_DensBZKmBl * xyz_dz_xyr( xyr_VelZBl ) ) & & / xyr_DensBZ & & + xyr_dx_pyr( & & pyr_xyz( xyz_KmBl ) * pyr_dz_pyz( pyz_VelXBl ) & & + pyr_xyz( xyz_KmBl ) * pyr_dx_xyr( xyr_VelZBl ) & & ) & & + xyr_dy_xqr( & & xqr_xyz( xyz_KmBl ) * xqr_dz_xqz( xqz_VelYBl ) & & + xqr_xyz( xyz_KmBl ) * xqr_dy_xyr( xyr_VelZBl ) & & ) & & - 2.0d0 * xyr_dz_xyz( xyz_DensBZKmBlKmBl ) & & / ( 3.0d0 * Cm_Cm_MixLen_MixLen ) & & / xyr_DensBZ end subroutine Turbulence_kw1978_std !!---------------------------------------------------------------------------!! subroutine turbulence_kw1978_center2_3d ! モジュール ! use axesset, only : dx, dy, dz use gridset, only : imin, &! x 方向の配列の下限 & imax, &! x 方向の配列の上限 & jmin, &! y 方向の配列の下限 & jmax, &! y 方向の配列の上限 & kmin, &! z 方向の配列の下限 & kmax, &! z 方向の配列の上限 & ncmax, &! & nx,ny,nz use basicset, only : xyz_ExnerBZ, & & xyz_DensBZ, & & pyz_DensBZ, & & xqz_DensBZ, & & xyr_DensBZ use constants, only : CpDry ! 暗黙の型宣言禁止 ! implicit none ! 作業変数 ! real(DP) :: xyz_KmBlKmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: pqz_KmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: pyr_KmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xqr_KmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DensBZKhBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DensBZKmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DensBZKmBlKmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: pyr_DensBZKmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xqr_DensBZKmBl(imin:imax, jmin:jmax, kmin:kmax) integer :: i, j, k, f ! Km の自乗 ! xyz_KmBlKmBl = xyz_KmBl * xyz_KmBl ! Km の格子点位置の変更 ! do k = kmin, kmax do j = jmin, jmax-1 do i = imin, imax-1 pqz_KmBl(i,j,k) = & & ( & & xyz_KmBl(i, j, k) & & + xyz_KmBl(i+1, j, k) & & + xyz_KmBl(i, j+1, k) & & + xyz_KmBl(i+1, j+1, k) & & ) * 0.25d0 end do end do end do ! Km の格子点位置の変更 ! do k = kmin, kmax-1 do j = jmin, jmax do i = imin, imax-1 pyr_KmBl(i,j,k) = & & ( & & + xyz_KmBl(i, j, k ) & & + xyz_KmBl(i+1, j, k ) & & + xyz_KmBl(i, j, k+1) & & + xyz_KmBl(i+1, j, k+1) & & ) * 0.25d0 end do end do end do ! Km の格子点位置の変更 ! do k = kmin, kmax-1 do j = jmin, jmax-1 do i = imin, imax xqr_KmBl(i,j,k) = & & ( & & + xyz_KmBl(i, j, k ) & & + xyz_KmBl(i, j+1, k ) & & + xyz_KmBl(i, j, k+1) & & + xyz_KmBl(i, j+1, k+1) & & ) * 0.25d0 end do end do end do ! Dens * Kh ! xyz_DensBZKhBl = xyz_DensBZ * xyz_KhBl ! Dens * Km ! xyz_DensBZKmBl = xyz_DensBZ * xyz_KmBl xyz_DensBZKmBlKmBl = xyz_DensBZ * xyz_KmBl * xyz_KmBl ! Dens * Km (pyr) ! do k = kmin, kmax-1 do j = jmin, jmax do i = imin, imax-1 pyr_DensBZKmBl(i,j,k) = & & ( & & + xyz_DensBZKmBl(i, j, k ) & & + xyz_DensBZKmBl(i+1, j, k ) & & + xyz_DensBZKmBl(i, j, k+1) & & + xyz_DensBZKmBl(i+1, j, k+1) & & ) * 0.25d0 end do end do end do ! Dens * Km (xqr) ! do k = kmin, kmax-1 do j = jmin, jmax-1 do i = imin, imax xqr_DensBZKmBl(i,j,k) = & & ( & & + xyz_DensBZKmBl(i, j, k ) & & + xyz_DensBZKmBl(i, j+1, k ) & & + xyz_DensBZKmBl(i, j, k+1) & & + xyz_DensBZKmBl(i, j+1, k+1) & & ) * 0.25d0 end do end do end do !--------------------------------------- ! Km の時間発展方程式の各項 ! ! 散逸項 ! xyz_Disp = - xyz_KmBlKmBl * 5.0d-1 / ( MixLen ** 2.0d0 ) ! 浮力項 ! call turbulence_kw1978_BuoyancyKm ! シアー項 ! do k = 1, nz do j = 1, ny do i = 1, nx xyz_Shear(i,j,k) = & & Cm_Cm_MixLen_MixLen & & * ( & & + ( & & ( & & pyz_VelXBl( i, j, k ) & & - pyz_VelXBl( i-1,j, k ) & & ) / dx & & ) ** 2.0d0 & & + ( & & ( & & xqz_VelYBl( i, j, k ) & & - xqz_VelYBl( i, j-1, k ) & & ) / dy & & ) ** 2.0d0 & & + ( & & ( & & xyr_VelZBl( i, j, k ) & & - xyr_VelZBl( i, j, k-1 ) & & ) / dz & & ) ** 2.0d0 & & + 5.0d-1 & & * ( & & ( & & + ( & & + pyz_VelXBl( i, j, k+1 ) & & + pyz_VelXBl( i-1, j, k+1 ) & & - pyz_VelXBl( i, j, k-1 ) & & - pyz_VelXBl( i-1, j, k-1 ) & & ) * 0.25d0 / dz & & + ( & & + xyr_VelZBl( i+1, j, k ) & & + xyr_VelZBl( i+1, j, k-1 ) & & - xyr_VelZBl( i-1, j, k ) & & - xyr_VelZBl( i-1, j, k-1 ) & & ) * 0.25d0 / dx & & ) ** 2.0d0 & & + ( & & + ( & & + xyr_VelZBl( i, j+1, k ) & & + xyr_VelZBl( i, j+1, k-1 ) & & - xyr_VelZBl( i, j-1, k ) & & - xyr_VelZBl( i, j-1, k-1 ) & & ) * 0.25d0 / dy & & + ( & & + xqz_VelYBl( i, j, k+1 ) & & + xqz_VelYBl( i, j-1, k+1 ) & & - xqz_VelYBl( i, j, k-1 ) & & - xqz_VelYBl( i, j-1, k-1 ) & & ) * 0.25d0 / dz & & ) ** 2.0d0 & & + ( & & + ( & & + xqz_VelYBl( i+1, j, k ) & & + xqz_VelYBl( i+1, j-1, k ) & & - xqz_VelYBl( i-1, j, k ) & & - xqz_VelYBl( i-1, j-1, k ) & & ) * 0.25d0 / dx & & + ( & & + pyz_VelXBl( i, j+1, k ) & & + pyz_VelXBl( i-1, j+1, k ) & & - pyz_VelXBl( i, j-1, k ) & & - pyz_VelXBl( i-1, j-1, k ) & & ) * 0.25d0 / dy & & ) ** 2.0d0 & & ) & & ) & & - xyz_KmBl( i, j, k ) & & * ( & & ( & & pyz_VelXBl( i, j, k ) & & - pyz_VelXBl( i-1, j, k ) & & ) / dx & & + ( & & xqz_VelYBl( i, j, k ) & & - xqz_VelYBl( i, j-1, k ) & & ) / dy & & + ( & & xyr_VelZBl( i, j, k ) & & - xyr_VelZBl( i, j, k-1 ) & & ) / dz & & ) / 3.0d0 end do end do end do ! 拡散項 ! do k = 1, nz do j = 1, ny do i = 1, nx xyz_Diff(i, j, k ) = & & + 5.0d-1 & & * ( & & + ( & & + xyz_KmBlKmBl( i+1, j, k ) & & + xyz_KmBlKmBl( i-1, j, k ) & & - xyz_KmBlKmBl( i, j, k ) * 2.0d0 & & ) / ( dx * dx ) & & + ( & & + xyz_KmBlKmBl( i, j+1, k ) & & + xyz_KmBlKmBl( i, j-1, k ) & & - xyz_KmBlKmBl( i, j, k ) * 2.0d0 & & ) / ( dy * dy ) & & + ( & & + xyz_KmBlKmBl( i, j, k+1 ) & & + xyz_KmBlKmBl( i, j, k-1 ) & & - xyz_KmBlKmBl( i, j, k ) * 2.0d0 & & ) / ( dz * dz ) & & ) & & + ( & & + ( & & ( & & + xyz_KmBl( i+1, j, k ) & & - xyz_KmBl( i-1, j, k ) & & ) * 0.5d0 / dx & & ) ** 2.0d0 & & + ( & & ( & & + xyz_KmBl( i, j+1, k ) & & - xyz_KmBl( i, j-1, k ) & & ) * 0.5d0 / dy & & ) ** 2.0d0 & & + ( & & ( & & + xyz_KmBl( i, j, k+1 ) & & - xyz_KmBl( i, j, k-1 ) & & ) * 0.5d0 / dz & & ) ** 2.0d0 & & ) end do end do end do !----------------------------------------------- ! Turb.\theta ! do k = 1, nz do j = 1, ny do i = 1, nx xyz_Turb(i, j, k ) = & & + ( & & + ( & & xyz_KhBl(i+1,j, k ) & & + xyz_KhBl(i, j, k ) & & ) & & * ( & & xyz_PTempBlAll(i+1,j, k ) & & - xyz_PTempBlAll(i, j, k ) & & ) & & - ( & & xyz_KhBl(i, j, k ) & & + xyz_KhBl(i-1,j, k ) & & ) & & * ( & & xyz_PTempBlAll(i, j, k ) & & - xyz_PTempBlAll(i-1,j, k ) & & ) & & ) * 0.5d0 / ( dx * dx ) & & + ( & & + ( & & xyz_KhBl(i, j+1,k ) & & + xyz_KhBl(i, j, k ) & & ) & & * ( & & xyz_PTempBlAll(i, j+1,k ) & & - xyz_PTempBlAll(i, j, k ) & & ) & & - ( & & xyz_KhBl(i, j, k ) & & + xyz_KhBl(i, j-1,k ) & & ) & & * ( & & xyz_PTempBlAll(i, j, k ) & & - xyz_PTempBlAll(i, j-1,k ) & & ) & & ) * 0.5d0 / ( dy * dy ) & & + ( & & + ( & & xyz_DensBZKhBl(i, j, k+1) & & + xyz_DensBZKhBl(i, j, k ) & & ) & & * ( & & xyz_PTempBlAll(i, j, k+1) & & - xyz_PTempBlAll(i, j, k ) & & ) & & - ( & & xyz_DensBZKhBl(i, j, k ) & & + xyz_DensBZKhBl(i, j, k-1) & & ) & & * ( & & xyz_PTempBlAll(i, j, k ) & & - xyz_PTempBlAll(i, j, k-1) & & ) & & ) * 0.5d0 / ( dz * dz ) & & / xyz_DensBZ(i, j, k ) end do end do end do xyz_DispHeat = (xyz_KmBl ** 3.0d0) & & / (xyz_ExnerBZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0)) !----------------------------------------------- ! Turb.Q ! do f = 1, ncmax do k = 1, nz do j = 1, ny do i = 1, nx xyzf_Turb(i, j, k, f ) = & & + ( & & + ( & & xyz_KhBl(i+1,j, k ) & & + xyz_KhBl(i, j, k ) & & ) & & * ( & & xyzf_QMixBlAll(i+1,j, k, f ) & & - xyzf_QMixBlAll(i, j, k, f ) & & ) & & - ( & & xyz_KhBl(i, j, k ) & & + xyz_KhBl(i-1,j, k ) & & ) & & * ( & & xyzf_QMixBlAll(i, j, k, f ) & & - xyzf_QMixBlAll(i-1,j, k, f ) & & ) & & ) * 0.5d0 / ( dx * dx ) & & + ( & & + ( & & xyz_KhBl(i, j+1,k ) & & + xyz_KhBl(i, j, k ) & & ) & & * ( & & xyzf_QMixBlAll(i, j+1,k, f ) & & - xyzf_QMixBlAll(i, j, k, f ) & & ) & & - ( & & xyz_KhBl(i, j, k ) & & + xyz_KhBl(i, j-1,k ) & & ) & & * ( & & xyzf_QMixBlAll(i, j, k, f ) & & - xyzf_QMixBlAll(i, j-1,k, f ) & & ) & & ) * 0.5d0 / ( dy * dy ) & & + ( & & + ( & & xyz_DensBZKhBl(i, j, k+1) & & + xyz_DensBZKhBl(i, j, k ) & & ) & & * ( & & xyzf_QMixBlAll(i, j, k+1,f ) & & - xyzf_QMixBlAll(i, j, k, f ) & & ) & & - ( & & xyz_DensBZKhBl(i, j, k ) & & + xyz_DensBZKhBl(i, j, k-1) & & ) & & * ( & & xyzf_QMixBlAll(i, j, k, f ) & & - xyzf_QMixBlAll(i, j, k-1,f ) & & ) & & ) * 0.5d0 / ( dz * dz ) & & / xyz_DensBZ(i, j, k ) end do end do end do end do !----------------------------------------------- ! Turb.u ! do k = 1, nz do j = 1, ny do i = 1, nx pyz_Turb(i,j,k) = & & + 2.0d0 & & * ( & & + xyz_KmBl(i+1,j, k ) & & * ( pyz_VelXBl(i+1,j, k ) - pyz_VelXBl(i, j, k ) ) & & - xyz_KmBl(i, j, k ) & & * ( pyz_VelXBl(i, j, k ) - pyz_VelXBl(i-1,j, k ) ) & & ) / ( dx * dx ) & & + ( & & + pqz_KmBl( i, j, k ) & & * ( & & ( xqz_VelYBl(i+1,j, k ) - xqz_VelYBl(i, j, k ) ) / dx & & + ( pyz_VelXBl(i, j+1,k ) - pyz_VelXBl(i, j, k ) ) / dy & & ) & & - pqz_KmBl( i, j-1,k ) & & * ( & & ( xqz_VelYBl(i+1,j-1,k ) - xqz_VelYBl(i, j-1,k ) ) / dx & & + ( pyz_VelXBl(i, j, k ) - pyz_VelXBl(i, j-1,k ) ) / dy & & ) & & ) / dy & & + ( & & + pyr_DensBZKmBl(i, j, k ) & & * ( & & ( xyr_VelZBl(i+1,j, k ) - xyr_VelZBl(i, j, k ) ) / dx & & + ( pyz_VelXBl(i, j, k+1) - pyz_VelXBl(i, j, k ) ) / dz & & ) & & - pyr_DensBZKmBl(i, j, k-1) & & * ( & & ( xyr_VelZBl(i+1,j, k-1) - xyr_VelZBl(i, j, k-1) ) / dx & & + ( pyz_VelXBl(i, j, k ) - pyz_VelXBl(i, j, k-1) ) / dz & & ) & & ) / dz / pyz_DensBZ(i, j, k) & & - 2.0d0 & & * ( xyz_KmBlKmBl(i+1,j, k ) - xyz_KmBlKmBl(i, j, k ) ) / dx & & / ( 3.0d0 * Cm_Cm_MixLen_MixLen ) end do end do end do !----------------------------------------------- ! Turb.v ! do k = 1, nz do j = 1, ny do i = 1, nx xqz_Turb(i,j,k) = & & + 2.0d0 & & * ( & & + xyz_KmBl(i, j+1,k ) & & * ( xqz_VelYBl(i, j+1,k ) - xqz_VelYBl(i, j, k ) ) & & - xyz_KmBl(i, j, k ) & & * ( xqz_VelYBl(i, j, k ) - xqz_VelYBl(i, j-1,k ) ) & & ) / ( dy * dy ) & & + ( & & + pqz_KmBl( i, j, k ) & & * ( & & ( pyz_VelXBl(i, j+1,k ) - pyz_VelXBl(i, j, k ) ) / dy & & + ( xqz_VelYBl(i+1,j, k ) - xqz_VelYBl(i, j, k ) ) / dx & & ) & & - pqz_KmBl( i-1,j, k ) & & * ( & & ( pyz_VelXBl(i-1,j+1,k ) - pyz_VelXBl(i-1,j, k ) ) / dy & & + ( xqz_VelYBl(i, j, k ) - xqz_VelYBl(i-1,j, k ) ) / dx & & ) & & ) / dx & & + ( & & + xqr_DensBZKmBl(i, j, k ) & & * ( & & ( xyr_VelZBl(i, j+1,k ) - xyr_VelZBl(i, j, k ) ) / dy & & + ( xqz_VelYBl(i, j, k+1) - xqz_VelYBl(i, j, k ) ) / dz & & ) & & - xqr_DensBZKmBl(i, j, k-1) & & * ( & & ( xyr_VelZBl(i, j+1,k-1) - xyr_VelZBl(i, j, k-1) ) / dy & & + ( xqz_VelYBl(i, j ,k ) - xqz_VelYBl(i, j, k-1) ) / dz & & ) & & ) / dz / xyr_DensBZ(i, j, k) & & - 2.0d0 & & * ( xyz_KmBlKmBl(i, j+1,k ) - xyz_KmBlKmBl(i, j, k ) ) / dy & & / ( 3.0d0 * Cm_Cm_MixLen_MixLen ) end do end do end do !----------------------------------------------- ! Turb.w ! do k = 1, nz do j = 1, ny do i = 1, nx xyr_Turb(i, j, k ) = & & + 2.0d0 & & * ( & & + xyz_DensBZKmBl(i, j, k+1) & & * ( xyr_VelZBl(i, j, k+1) - xyr_VelZBl(i, j, k ) ) & & - xyz_DensBZKmBl(i, j, k ) & & * ( xyr_VelZBl(i, j, k ) - xyr_VelZBl(i, j, k-1) ) & & ) / ( dz * dz ) & & / xyr_DensBZ(i, j, k ) & & + ( & & + pyr_KmBl( i, j, k ) & & * ( & & ( pyz_VelXBl(i, j, k+1) - pyz_VelXBl(i, j, k ) ) / dz & & + ( xyr_VelZBl(i+1,j, k ) - xyr_VelZBl(i, j, k ) ) / dx & & ) & & - pyr_KmBl( i-1,j, k ) & & * ( & & ( pyz_VelXBl(i-1,j, k+1) - pyz_VelXBl(i-1,j, k ) ) / dz & & + ( xyr_VelZBl(i, j, k ) - xyr_VelZBl(i-1,j, k ) ) / dx & & ) & & ) / dx & & + ( & & + xqr_KmBl(i, j, k ) & & * ( & & ( xqz_VelYBl(i, j, k+1) - xqz_VelYBl(i, j, k ) ) / dz & & + ( xyr_VelZBl(i, j+1,k ) - xyr_VelZBl(i, j, k ) ) / dy & & ) & & - xqr_KmBl(i, j-1,k ) & & * ( & & ( xqz_VelYBl(i, j-1,k+1) - xqz_VelYBl(i, j-1,k ) ) / dz & & + ( xyr_VelZBl(i, j, k ) - xyr_VelZBl(i, j-1,k ) ) / dy & & ) & & ) / dy & & - 2.0d0 & & * ( & & xyz_DensBZKmBlKmBl(i, j, k+1) & & - xyz_DensBZKmBlKmBl(i, j, k ) & & ) / dz & & / ( 3.0d0 * Cm_Cm_MixLen_MixLen ) & & / xyr_DensBZ(i, j, k ) end do end do end do end subroutine Turbulence_kw1978_center2_3d !!---------------------------------------------------------------------------!! subroutine turbulence_kw1978_center2_2d ! モジュール ! use axesset, only : dx, dz use gridset, only : imin, &! x 方向の配列の下限 & imax, &! x 方向の配列の上限 & jmin, &! y 方向の配列の下限 & jmax, &! y 方向の配列の上限 & kmin, &! z 方向の配列の下限 & kmax, &! z 方向の配列の上限 & ncmax, &! & nx, nz use basicset, only : xyz_ExnerBZ, & & xyz_DensBZ, & & pyz_DensBZ, & & xyr_DensBZ use constants, only : CpDry ! 暗黙の型宣言禁止 ! implicit none ! 作業変数 ! real(DP) :: xyz_KmBlKmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: pyr_KmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DensBZKhBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DensBZKmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: xyz_DensBZKmBlKmBl(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: pyr_DensBZKmBl(imin:imax, jmin:jmax, kmin:kmax) integer :: i, k, f integer, parameter :: j = 1 ! Km の自乗 ! xyz_KmBlKmBl = xyz_KmBl * xyz_KmBl ! Km の格子点位置の変更 ! do k = kmin, kmax-1 do i = imin, imax-1 pyr_KmBl(i,j,k) = & & ( & & + xyz_KmBl(i, j, k ) & & + xyz_KmBl(i+1, j, k ) & & + xyz_KmBl(i, j, k+1) & & + xyz_KmBl(i+1, j, k+1) & & ) * 0.25d0 end do end do ! Dens * Kh ! xyz_DensBZKhBl = xyz_DensBZ * xyz_KhBl ! Dens * Km ! xyz_DensBZKmBl = xyz_DensBZ * xyz_KmBl xyz_DensBZKmBlKmBl = xyz_DensBZ * xyz_KmBl * xyz_KmBl ! Dens * Km (pyr) ! do k = kmin, kmax-1 do i = imin, imax-1 pyr_DensBZKmBl(i,j,k) = & & ( & & + xyz_DensBZKmBl(i, j, k ) & & + xyz_DensBZKmBl(i+1, j, k ) & & + xyz_DensBZKmBl(i, j, k+1) & & + xyz_DensBZKmBl(i+1, j, k+1) & & ) * 0.25d0 end do end do !--------------------------------------- ! Km の時間発展方程式の各項 ! ! 散逸項 ! xyz_Disp = - xyz_KmBlKmBl * 5.0d-1 / ( MixLen ** 2.0d0 ) ! 浮力項 ! call turbulence_kw1978_BuoyancyKm ! シアー項 ! do k = 1, nz do i = 1, nx xyz_Shear(i,j,k) = & & Cm_Cm_MixLen_MixLen & & * ( & & + ( & & ( & & pyz_VelXBl( i, j, k ) & & - pyz_VelXBl( i-1,j, k ) & & ) / dx & & ) ** 2.0d0 & & + ( & & ( & & xyr_VelZBl( i, j, k ) & & - xyr_VelZBl( i, j, k-1 ) & & ) / dz & & ) ** 2.0d0 & & + 5.0d-1 & & * ( & & ( & & + ( & & + pyz_VelXBl( i, j, k+1 ) & & + pyz_VelXBl( i-1, j, k+1 ) & & - pyz_VelXBl( i, j, k-1 ) & & - pyz_VelXBl( i-1, j, k-1 ) & & ) * 0.25d0 / dz & & + ( & & + xyr_VelZBl( i+1, j, k ) & & + xyr_VelZBl( i+1, j, k-1 ) & & - xyr_VelZBl( i-1, j, k ) & & - xyr_VelZBl( i-1, j, k-1 ) & & ) * 0.25d0 / dx & & ) ** 2.0d0 & & ) & & ) & & - xyz_KmBl( i, j, k ) & & * ( & & ( & & pyz_VelXBl( i, j, k ) & & - pyz_VelXBl( i-1, j, k ) & & ) / dx & & + ( & & xyr_VelZBl( i, j, k ) & & - xyr_VelZBl( i, j, k-1 ) & & ) / dz & & ) / 3.0d0 end do end do ! 拡散項 ! do k = 1, nz do i = 1, nx xyz_Diff(i, j, k ) = & & + 5.0d-1 & & * ( & & + ( & & + xyz_KmBlKmBl( i+1, j, k ) & & + xyz_KmBlKmBl( i-1, j, k ) & & - xyz_KmBlKmBl( i, j, k ) * 2.0d0 & & ) / ( dx * dx ) & & + ( & & + xyz_KmBlKmBl( i, j, k+1 ) & & + xyz_KmBlKmBl( i, j, k-1 ) & & - xyz_KmBlKmBl( i, j, k ) * 2.0d0 & & ) / ( dz * dz ) & & ) & & + ( & & + ( & & ( & & + xyz_KmBl( i+1, j, k ) & & - xyz_KmBl( i-1, j, k ) & & ) * 0.5d0 / dx & & ) ** 2.0d0 & & + ( & & ( & & + xyz_KmBl( i, j, k+1 ) & & - xyz_KmBl( i, j, k-1 ) & & ) * 0.5d0 / dz & & ) ** 2.0d0 & & ) end do end do !----------------------------------------------- ! Turb.\theta ! do k = 1, nz do i = 1, nx xyz_Turb(i, j, k ) = & & + ( & & + ( & & xyz_KhBl(i+1,j, k ) & & + xyz_KhBl(i, j, k ) & & ) & & * ( & & xyz_PTempBlAll(i+1,j, k ) & & - xyz_PTempBlAll(i, j, k ) & & ) & & - ( & & xyz_KhBl(i, j, k ) & & + xyz_KhBl(i-1,j, k ) & & ) & & * ( & & xyz_PTempBlAll(i, j, k ) & & - xyz_PTempBlAll(i-1,j, k ) & & ) & & ) * 0.5d0 / ( dx * dx ) & & + ( & & + ( & & xyz_DensBZKhBl(i, j, k+1) & & + xyz_DensBZKhBl(i, j, k ) & & ) & & * ( & & xyz_PTempBlAll(i, j, k+1) & & - xyz_PTempBlAll(i, j, k ) & & ) & & - ( & & xyz_DensBZKhBl(i, j, k ) & & + xyz_DensBZKhBl(i, j, k-1) & & ) & & * ( & & xyz_PTempBlAll(i, j, k ) & & - xyz_PTempBlAll(i, j, k-1) & & ) & & ) * 0.5d0 / ( dz * dz ) & & / xyz_DensBZ(i, j, k ) end do end do xyz_DispHeat = (xyz_KmBl ** 3.0d0) & & / (xyz_ExnerBZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0)) !----------------------------------------------- ! Turb.Q ! do f = 1, ncmax do k = 1, nz do i = 1, nx xyzf_Turb(i, j, k, f ) = & & + ( & & + ( & & xyz_KhBl(i+1,j, k ) & & + xyz_KhBl(i, j, k ) & & ) & & * ( & & xyzf_QMixBlAll(i+1,j, k, f ) & & - xyzf_QMixBlAll(i, j, k, f ) & & ) & & - ( & & xyz_KhBl(i, j, k ) & & + xyz_KhBl(i-1,j, k ) & & ) & & * ( & & xyzf_QMixBlAll(i, j, k, f ) & & - xyzf_QMixBlAll(i-1,j, k, f ) & & ) & & ) * 0.5d0 / ( dx * dx ) & & + ( & & + ( & & xyz_DensBZKhBl(i, j, k+1) & & + xyz_DensBZKhBl(i, j, k ) & & ) & & * ( & & xyzf_QMixBlAll(i, j, k+1,f ) & & - xyzf_QMixBlAll(i, j, k, f ) & & ) & & - ( & & xyz_DensBZKhBl(i, j, k ) & & + xyz_DensBZKhBl(i, j, k-1) & & ) & & * ( & & xyzf_QMixBlAll(i, j, k, f ) & & - xyzf_QMixBlAll(i, j, k-1,f ) & & ) & & ) * 0.5d0 / ( dz * dz ) & & / xyz_DensBZ(i, j, k ) end do end do end do !----------------------------------------------- ! Turb.u ! do k = 1, nz do i = 1, nx pyz_Turb(i,j,k) = & & + 2.0d0 & & * ( & & + xyz_KmBl(i+1,j, k ) & & * ( pyz_VelXBl(i+1,j, k ) - pyz_VelXBl(i, j, k ) ) & & - xyz_KmBl(i, j, k ) & & * ( pyz_VelXBl(i, j, k ) - pyz_VelXBl(i-1,j, k ) ) & & ) / ( dx * dx ) & & + ( & & + pyr_DensBZKmBl(i, j, k ) & & * ( & & ( xyr_VelZBl(i+1,j, k ) - xyr_VelZBl(i, j, k ) ) / dx & & + ( pyz_VelXBl(i, j, k+1) - pyz_VelXBl(i, j, k ) ) / dz & & ) & & - pyr_DensBZKmBl(i, j, k-1) & & * ( & & ( xyr_VelZBl(i+1,j, k-1) - xyr_VelZBl(i, j, k-1) ) / dx & & + ( pyz_VelXBl(i, j, k ) - pyz_VelXBl(i, j, k-1) ) / dz & & ) & & ) / dz / pyz_DensBZ(i, j, k) & & - 2.0d0 & & * ( xyz_KmBlKmBl(i+1,j, k ) - xyz_KmBlKmBl(i, j, k ) ) / dx & & / ( 3.0d0 * Cm_Cm_MixLen_MixLen ) end do end do !----------------------------------------------- ! Turb.v ! xqz_Turb = 0.0d0 !----------------------------------------------- ! Turb.w ! do k = 1, nz do i = 1, nx xyr_Turb(i, j, k ) = & & + 2.0d0 & & * ( & & + xyz_DensBZKmBl(i, j, k+1) & & * ( xyr_VelZBl(i, j, k+1) - xyr_VelZBl(i, j, k ) ) & & - xyz_DensBZKmBl(i, j, k ) & & * ( xyr_VelZBl(i, j, k ) - xyr_VelZBl(i, j, k-1) ) & & ) / ( dz * dz ) & & / xyr_DensBZ(i, j, k ) & & + ( & & + pyr_KmBl( i, j, k ) & & * ( & & ( pyz_VelXBl(i, j, k+1) - pyz_VelXBl(i, j, k ) ) / dz & & + ( xyr_VelZBl(i+1,j, k ) - xyr_VelZBl(i, j, k ) ) / dx & & ) & & - pyr_KmBl( i-1,j, k ) & & * ( & & ( pyz_VelXBl(i-1,j, k+1) - pyz_VelXBl(i-1,j, k ) ) / dz & & + ( xyr_VelZBl(i, j, k ) - xyr_VelZBl(i-1,j, k ) ) / dx & & ) & & ) / dx & & - 2.0d0 & & * ( & & xyz_DensBZKmBlKmBl(i, j, k+1) & & - xyz_DensBZKmBlKmBl(i, j, k ) & & ) / dz & & / ( 3.0d0 * Cm_Cm_MixLen_MixLen ) & & / xyr_DensBZ(i, j, k ) end do end do end subroutine Turbulence_kw1978_center2_2d !!!------------------------------------------------------------------------------!!! subroutine turbulence_kw1978_BuoyancyKm ! use ChemCalc, only: xyz_LatentHeat !潜熱 use basicset, only: xyz_PTempBZ, &!基本場の温位 & xyz_VPTempBZ, &! 仮温位 & xyz_QMixBZ, &!基本場の混合比 & xyz_QMixBZPerMolWt, &!基本場の混合比 & xyz_EffMolWtBZ, &!基本場の混合比 & xyz_ExnerBZ !基本場のエクスナー関数 use constants,only: Grav, & & MolWtDry, & & CpDry use composition, only: MolWtWet, & & SpcWetID, & & CondNum, &!凝結過程の数 & GasNum, &!気体の数 & IdxCG, &!凝結過程(蒸気)の配列添え字 & IdxCC, &!凝結過程(雲)の配列添え字 & IdxG !気体の配列添え字 use xyz_deriv_module, only: & & xyr_dz_xyz use axesset, only : xyz_xyr !暗黙の型宣言禁止 implicit none !変数定義 ! real(DP) :: xyzf_LatentHeat(imin:imax,jmin:jmax,kmin:kmax, GasNum) !潜熱 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 !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 ! xyzf_LatentHeat = 0.0d0 !作業配列の初期化. 気体のみ利用 ! do s = 1, GasNum xyzf_QMixPerMolWt(:,:,:,s) = xyzf_QMixBl(:,:,:,IdxG(s)) / MolWtWet(IdxG(s)) end do !温度の効果 ! xyz_EffPTemp = xyz_PTempBl / xyz_PTempBZ !分子量効果 + 引きづりの効果 ! xyz_EffMolWt = & & + sum( xyzf_QMixPerMolWt, 4 ) & & / ( 1.0d0 / MolWtDry + xyz_QMixBZPerMolWt ) & & - sum(xyzf_QMixBl, 4) / ( 1.0d0 + xyz_QMixBZ ) !蒸気が蒸発する場合の潜熱を計算 ! 分子量の部分はいつでも効くが潜熱は飽和していないと効かないので, ! 雲の混合比がゼロの時には, 潜熱の寄与はゼロとなるように調節している ! do s = 1, CondNum xyzf_LatentHeat(:,:,:,s) = & & xyz_LatentHeat( SpcWetID(IdxCC(s)), xyz_TempBlAll ) & & * xyzf_QMixBlAll(:,:,:,IdxCG(s)) & & * ( 5.0d-1 + sign( 5.0d-1, (xyzf_QMixBlAll(:,:,:,IdxCC(s)) - 1.0d-4) ) ) end do xyz_EffHeat = ( sum( xyzf_LatentHeat, 4 ) * xyz_EffMolWtBZ ) & & / ( CpDry * xyz_ExnerBZ ) !乱流拡散係数の時間発展式の浮力項を決める ! xyz_Buoy = & & - 3.0d0 * Grav * Cm_Cm_MixLen_MixLen & & * xyz_xyr( & & xyr_dz_xyz( & & xyz_EffHeat & & + xyz_VPTempBZ * ( 1.0d0 + xyz_EffPTemp + xyz_EffMolWt ) & & ) & & ) & & / ( 2.0d0 * xyz_VPTempBZ ) xyz_BuoyT = & & - 3.0d0 * Grav * Cm_Cm_MixLen_MixLen & & * xyz_xyr( & & xyr_dz_xyz( & & xyz_PTempBlAll & & ) & & ) & & / ( 2.0d0 * xyz_PTempBZ ) xyz_BuoyM = xyz_Buoy - xyz_BuoyT end subroutine Turbulence_kw1978_BuoyancyKm end subroutine Turbulence_KW1978_forcing !!!------------------------------------------------------------------------!!! subroutine turbulence_kw1978_output ! モジュール ! use gtool_historyauto, only : HistoryAutoAddVariable use gridset, only : ncmax use composition, only : SpcWetSymbol ! 暗黙の型宣言禁止 ! implicit none ! 作業変数 ! integer :: l call HistoryAutoAddVariable( & & varname='VelXTurb',& & dims=(/'x','y','z','t'/), & & longname='Turbulence term of velocity (x)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelYTurb',& & dims=(/'x','y','z','t'/), & & longname='Turbulence term of velocity (y)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelZTurb',& & dims=(/'x','y','z','t'/), & & longname='Turbulence term of velocity (z)', & & units='m.s-2', & & 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='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 call HistoryAutoAddVariable( & & varname='KmDiff',& & dims=(/'x','y','z','t'/), & & longname='Diffusion term of Km', & & units='s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='KmBuoyT',& & dims=(/'x','y','z','t'/), & & longname='Buoyancy term of Km (temperature)', & & units='s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='KmBuoyM',& & dims=(/'x','y','z','t'/), & & longname='Buoyancy term of Km (composition)', & & units='s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='KmShear',& & dims=(/'x','y','z','t'/), & & longname='Shear term of Km', & & units='s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='KmDisp',& & dims=(/'x','y','z','t'/), & & longname='Dissipation term of Km', & & units='s-1', & & xtype='float') end subroutine turbulence_kw1978_output end module Turbulence_kw1978