!= Module DynamicsHEVI ! ! Authors:: 杉山耕一朗(SUGIYAMA Ko-ichiro) ! Version:: $Id: dynamics_hevi_v3.f90,v 1.4 2014/07/08 00:55:26 sugiyama Exp $ ! Tag Name:: $Name: $ ! Copyright:: Copyright (C) GFD Dennou Club, 2014. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] module DynamicsHEVI ! ! 力学コア. ! タイムスプリット法を利用. 音波モードとそれ以外を別々の時間刻みで解く. ! 短い時間ステップの計算には HE-VI 法を利用. ! !モジュール読み込み ! use dc_types, only : DP !暗黙の型宣言禁止 ! implicit none !属性の指定 ! private ! 変数の定義 ! character(*), parameter:: module_name = 'dynamics_hevi' ! モジュールの名称. ! Module name real(DP), save :: FactorBuoyTemp = 1.0d0 !浮力 (温度の寄与) の有無 !考慮しない場合は値をゼロにする. real(DP), save :: FactorBuoyMolWt = 1.0d0 !浮力 (分子量効果) の有無 !考慮しない場合は値をゼロにする. real(DP), save :: FactorBuoyLoading = 1.0d0 !浮力 (荷重効果) の有無 !考慮しない場合は値をゼロにする. real(DP), save :: FactorExnerAdv = 1.0d0 !エクスナー関数の移流の有無 !考慮しない場合は値をゼロにする. real(DP), save :: FactorExnerExpnd = 1.0d0 !エクスナー関数の膨張項の有無 !考慮しない場合は値をゼロにする. integer, save :: IDAdvection = 0 ! 移流計算方法 integer, parameter :: IDAdvectionCenter4_std = 1 ! 1: 4 次中央差分 (微分平均モジュール利用) integer, parameter :: IDAdvectionCenter4_2D_dry = 2 ! 2: 4 次中央差分 (2D 版計算ルーチン) integer, parameter :: IDAdvectionCenter4_2D = 3 ! 2: 4 次中央差分 (2D 版計算ルーチン) integer, parameter :: IDAdvectionCenter4_3D_dry = 4 ! 3: 4 次中央差分 (3D 版計算ルーチン) integer, parameter :: IDAdvectionCenter4_3D = 5 ! 3: 4 次中央差分 (3D 版計算ルーチン) integer, save :: IDAcousticmode = 0 ! 短い時間ステップの計算方法 integer, parameter :: IDAcousticmode_std = 1 ! 1: 微分平均モジュールの利用. integer, parameter :: IDAcousticmode_2D = 2 ! 2: 2D 版計算ルーチン integer, parameter :: IDAcousticmode_3D = 3 ! 3: 3D 版計算ルーチン ! public ! public Dynamics_Init public Dynamics_Long_integrate public Dynamics_Long_forcing public Dynamics_Short_forcing contains subroutine Dynamics_Init ! ! 力学コア 3D 版の初期化ルーチン ! ! モジュール読み込み ! use dc_types, only : DP, STRING use dc_iounit, only : FileOpen use dc_message, only : MessageNotify use mpi_wrapper,only : myrank use gridset, only : FlagCalc3D, FlagCalcMoist use namelist_util, & & only : namelist_filename use advection_center4_std, & & only : advection_center4_std_init use advection_center4_2D, & & only : advection_center4_2D_init use advection_center4_3D, & & only : advection_center4_3D_init use acousticmode_std, & & only : acousticmode_std_init use acousticmode_2D, & & only : acousticmode_2D_init use acousticmode_3D, & & only : acousticmode_3D_init !暗黙の型宣言禁止 ! implicit none ! 変数の定義 ! real(DP) :: AlphaSound = 5.0d-2 !音波減衰項の係数 (気象庁数値予報課報告・別冊49 より) real(DP) :: AlphaNDiff = 1.0d-3 !4次の数値拡散の係数. CReSS マニュアルより real(DP) :: NDiffRatio = 1.0d0 !速度に対する粘性を上げる場合は数字を 1 以上にする. integer :: unit !装置番号 character(STRING) :: FlagAdvection = "" !移流計算方法の選択 character(STRING) :: FlagAcousticmode = "" !音波モードの計算方法の選択 !------------------------------------------------------------------- ! Namelist から情報を取得する ! NAMELIST /Dynamics_nml/ & & AlphaSound, AlphaNDiff, NDiffRatio, & & FactorBuoyTemp, FactorBuoyMolWt, FactorBuoyLoading, & & FactorExnerAdv, FactorExnerExpnd, & & FlagAdvection, FlagAcousticmode call FileOpen(unit, file=namelist_filename, mode='r') read(unit, NML=dynamics_nml) close(unit) !------------------------------------------------------------------- ! フラグ処理 ! if ( FlagAdvection == "Center4" .OR. FlagAdvection == "" ) then ! ! 4 次精度中心差分. 微分平均モジュールを使わない場合. (デフォルト) ! if ( FlagCalc3D ) then ! if ( FlagCalcTracer .OR. FlagCalcMoist ) then if ( FlagCalcMoist ) then IDAdvection = IDAdvectionCenter4_3D else IDAdvection = IDAdvectionCenter4_3D_dry end if else ! if ( FlagCalcTracer .OR. FlagCalcMoist ) then if ( FlagCalcMoist ) then IDAdvection = IDAdvectionCenter4_2D else IDAdvection = IDAdvectionCenter4_2D_dry end if end if else ! ! 4 次精度中心差分 with 微分平均モジュール ! IDAdvection = IDAdvectionCenter4_std end if if ( FlagAcousticmode == "Center2" .OR. FlagAcousticmode == "" ) then ! ! 微分平均モジュールを使わない場合. (デフォルト) ! if ( FlagCalc3D ) then IDAcousticmode = IDAcousticmode_3D else IDAcousticmode = IDAcousticmode_2D end if else ! ! with 微分平均モジュール ! IDAcousticmode = IDAcousticmode_std end if !------------------------------------------------------------------- ! 移流計算用の計算モジュールの初期化: ! select case ( IDAdvection ) case ( IDAdvectionCenter4_std ) call advection_center4_std_init( AlphaNDiff, NDiffRatio ) case ( IDAdvectionCenter4_3D_dry, IDAdvectionCenter4_3D ) call advection_center4_3D_init( AlphaNDiff, NDiffRatio ) case ( IDAdvectionCenter4_2D_dry, IDAdvectionCenter4_2D ) call advection_center4_2D_init( AlphaNDiff, NDiffRatio ) end select !------------------------------------------------------------------- ! 音波モードの計算モジュールの初期化: ! select case ( IDAcousticmode ) case ( IDAcousticmode_std ) call acousticmode_std_init( AlphaSound ) case ( IDAcousticmode_3D ) call acousticmode_3D_init( AlphaSound ) case ( IDAcousticmode_2D ) call acousticmode_2D_init( AlphaSound ) end select !------------------------------------------------------------------- ! tendency の出力 ! call Dynamics_Tendency_output !------------------------------------------------------------------- ! 出力 ! if (myrank == 0) then call MessageNotify( "M", module_name, "AlphaSound = %f", d=(/AlphaSound/) ) call MessageNotify( "M", module_name, "AlphaNDiff = %f", d=(/AlphaNDiff/) ) call MessageNotify( "M", module_name, "FactorBuoyTemp = %f", d=(/FactorBuoyTemp/) ) call MessageNotify( "M", module_name, "FactorBuoyMolWt = %f", d=(/FactorBuoyMolWt/) ) call MessageNotify( "M", module_name, "FactorBuoyLoading= %f", d=(/FactorBuoyLoading/) ) call MessageNotify( "M", module_name, "FactorExnerAdv = %f", d=(/FactorExnerAdv/) ) call MessageNotify( "M", module_name, "FactorExnerExpnd = %f", d=(/FactorExnerExpnd/) ) call MessageNotify( "M", module_name, "IDAdvection = %d", i=(/IDAdvection/) ) call MessageNotify( "M", module_name, "IDAcousticmode = %d", i=(/IDAcousticmode/) ) end if end subroutine Dynamics_Init !!!-------------------------------------------------------------------!!! subroutine Dynamics_Long_forcing( & & pyz_VelXBl, pyz_VelXNl, & ! (in) & xqz_VelYBl, xqz_VelYNl, & ! (in) & xyr_VelZBl, xyr_VelZNl, & ! (in) & xyz_PTempBl, xyz_PTempNl, & ! (in) & xyz_ExnerBl, xyz_ExnerNl, & ! (in) & xyzf_QMixBl, xyzf_QMixNl, & ! (in) & xyz_KmBl, xyz_KmNl, & ! (in) & pyz_DVelXDtNl, & ! (inout) & xqz_DVelYDtNl, & ! (inout) & xyr_DVelZDtNl, & ! (inout) & xyz_DPTempDtNl, & ! (inout) & xyz_DExnerDtNl, & ! (inout) & xyzf_DQMixDtNl, & ! (inout) & xyz_DKmDtNl & ! (inout) & ) ! ! 長い時間ステップで評価する移流項, 浮力項の計算. ! ! モジュール読み込み ! use dc_types, only : DP use gridset, only : FlagCalcMoist, & & imin, &! x 方向の配列の下限 & imax, &! x 方向の配列の上限 & jmin, &! y 方向の配列の下限 & jmax, &! y 方向の配列の上限 & kmin, &! z 方向の配列の下限 & kmax, &! z 方向の配列の上限 & nx, ny, nz, &! 物理領域の配列の上限 & ncmax ! 物質数 use composition, & & only : SpcWetSymbol use timeset, only : TimeN use gtool_historyauto, only : HistoryAutoPut use advection_center4_std, & & only : advection_center4_std_main use advection_center4_2D, & & only : advection_center4_2D_dry, advection_center4_2D_tracer use advection_center4_3D, & & only : advection_center4_3D_dry, advection_center4_3D_tracer ! 暗黙の型宣言禁止 ! implicit none ! 変数の定義 ! real(DP), intent(in) :: pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_PTempBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_PTempNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_ExnerBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_ExnerNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyzf_QMixBl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP), intent(in) :: xyzf_QMixNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP), intent(in) :: xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_KmNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: pyz_DVelXDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: xqz_DVelYDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: xyz_DPTempDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: xyz_DExnerDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: xyzf_DQMixDtNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP), intent(inout) :: xyz_DKmDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_VelXAdv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_VelYAdv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_VelZAdv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_KmAdv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_ExnerAdv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_PTempAdv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyzf_QMixAdv(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) real(DP) :: pyz_VelXnDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_VelYnDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_VelZnDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_KmNDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_ExnerNDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_PTempNDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyzf_QMixNDiff(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) real(DP) :: xyr_BuoyT(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_BuoyM(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_BuoyD(imin:imax,jmin:jmax,kmin:kmax) integer :: f !-------------------------------------------------------------------- ! 移流項の計算 ! select case ( IDAdvection ) case ( IDAdvectionCenter4_std ) call advection_center4_std_main( & & pyz_VelXBl, pyz_VelXNl, & ! (in) & xqz_VelYBl, xqz_VelYNl, & ! (in) & xyr_VelZBl, xyr_VelZNl, & ! (in) & xyz_PTempBl, xyz_PTempNl, & ! (in) & xyz_ExnerBl, xyz_ExnerNl, & ! (in) & xyzf_QMixBl, xyzf_QMixNl, & ! (in) & xyz_KmBl, xyz_KmNl, & ! (in) & pyz_VelXAdv, pyz_VelXnDiff, & !(out) & xqz_VelYAdv, xqz_VelYnDiff, & !(out) & xyr_VelZAdv, xyr_VelZnDiff, & !(out) & xyz_PTempAdv, xyz_PTempNDiff, & !(out) & xyz_ExnerAdv, xyz_ExnerNDiff, & !(out) & xyzf_QMixAdv, xyzf_QMixNDiff, & !(out) & xyz_KmAdv, xyz_KmNDiff & !(out) & ) case ( IDAdvectionCenter4_3D_dry ) call advection_center4_3D_dry( & & pyz_VelXBl, pyz_VelXNl, & ! (in) & xqz_VelYBl, xqz_VelYNl, & ! (in) & xyr_VelZBl, xyr_VelZNl, & ! (in) & xyz_PTempBl, xyz_PTempNl, & ! (in) & xyz_ExnerBl, xyz_ExnerNl, & ! (in) & xyz_KmBl, xyz_KmNl, & ! (in) & pyz_VelXAdv, pyz_VelXnDiff, & !(out) & xqz_VelYAdv, xqz_VelYnDiff, & !(out) & xyr_VelZAdv, xyr_VelZnDiff, & !(out) & xyz_PTempAdv, xyz_PTempNDiff, & !(out) & xyz_ExnerAdv, xyz_ExnerNDiff, & !(out) & xyz_KmAdv, xyz_KmNDiff & !(out) & ) case ( IDAdvectionCenter4_3D ) call advection_center4_3D_dry( & & pyz_VelXBl, pyz_VelXNl, & ! (in) & xqz_VelYBl, xqz_VelYNl, & ! (in) & xyr_VelZBl, xyr_VelZNl, & ! (in) & xyz_PTempBl, xyz_PTempNl, & ! (in) & xyz_ExnerBl, xyz_ExnerNl, & ! (in) & xyz_KmBl, xyz_KmNl, & ! (in) & pyz_VelXAdv, pyz_VelXnDiff, & !(out) & xqz_VelYAdv, xqz_VelYnDiff, & !(out) & xyr_VelZAdv, xyr_VelZnDiff, & !(out) & xyz_PTempAdv, xyz_PTempNDiff, & !(out) & xyz_ExnerAdv, xyz_ExnerNDiff, & !(out) & xyz_KmAdv, xyz_KmNDiff & !(out) & ) call advection_center4_3D_tracer( & & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & ! (in) & xyzf_QMixBl, xyzf_QMixNl, & ! (in) & xyzf_QMixAdv, xyzf_QMixNDiff & !(out) & ) case ( IDAdvectionCenter4_2D_dry ) ! Y 方向は tendency は零. ! xqz_VelYAdv = 0.0d0 xqz_VelYnDiff = 0.0d0 call advection_center4_2D_dry( & & pyz_VelXBl, pyz_VelXNl, & ! (in) & xyr_VelZBl, xyr_VelZNl, & ! (in) & xyz_PTempBl, xyz_PTempNl, & ! (in) & xyz_ExnerBl, xyz_ExnerNl, & ! (in) & xyz_KmBl, xyz_KmNl, & ! (in) & pyz_VelXAdv, pyz_VelXnDiff, & !(out) & xyr_VelZAdv, xyr_VelZnDiff, & !(out) & xyz_PTempAdv, xyz_PTempNDiff, & !(out) & xyz_ExnerAdv, xyz_ExnerNDiff, & !(out) & xyz_KmAdv, xyz_KmNDiff & !(out) & ) case ( IDAdvectionCenter4_2D ) ! Y 方向は tendency は零. ! xqz_VelYAdv = 0.0d0 xqz_VelYnDiff = 0.0d0 call advection_center4_2D_dry( & & pyz_VelXBl, pyz_VelXNl, & ! (in) & xyr_VelZBl, xyr_VelZNl, & ! (in) & xyz_PTempBl, xyz_PTempNl, & ! (in) & xyz_ExnerBl, xyz_ExnerNl, & ! (in) & xyz_KmBl, xyz_KmNl, & ! (in) & pyz_VelXAdv, pyz_VelXnDiff, & !(out) & xyr_VelZAdv, xyr_VelZnDiff, & !(out) & xyz_PTempAdv, xyz_PTempNDiff, & !(out) & xyz_ExnerAdv, xyz_ExnerNDiff, & !(out) & xyz_KmAdv, xyz_KmNDiff & !(out) & ) call advection_center4_2D_tracer( & & pyz_VelXNl, xyr_VelZNl, & ! (in) & xyzf_QMixBl, xyzf_QMixNl, & ! (in) & xyzf_QMixAdv, xyzf_QMixNDiff & !(out) & ) end select !-------------------------------------------------------------------- ! tendency の更新 ! ! 拡散係数 ! xyz_DKmDtNl = xyz_DKmDtNl + xyz_KmAdv + xyz_KmNDiff ! output ! call HistoryAutoPut(TimeN, 'KmAdv', xyz_KmAdv(1:nx,1:ny,1:nz) ) call HistoryAutoPut(TimeN, 'KmNDiff', xyz_KmNDiff(1:nx,1:ny,1:nz) ) ! エクスナー関数 ! 移流項を考慮しない場合は, FactorExnerAdv = 0 とする. ! xyz_DExnerDtNl = & & xyz_DExnerDtNl & & + ( xyz_ExnerAdv + xyz_ExnerNDiff ) * FactorExnerAdv ! output ! call HistoryAutoPut(TimeN, 'ExnerAdv', xyz_ExnerAdv(1:nx,1:ny,1:nz) ) call HistoryAutoPut(TimeN, 'ExnerNDiff', xyz_ExnerNDiff(1:nx,1:ny,1:nz) ) ! 温位 ! xyz_DPTempDtNl = xyz_DPTempDtNl + xyz_PTempAdv + xyz_PTempNDiff ! output ! call HistoryAutoPut(TimeN, 'PTempAdv', xyz_PTempAdv(1:nx,1:ny,1:nz) ) call HistoryAutoPut(TimeN, 'PTempNDiff', xyz_PTempNDiff(1:nx,1:ny,1:nz) ) ! 混合比 ! xyzf_DQMixDtNl = xyzf_DQMixDtNl + xyzf_QMixAdv + xyzf_QMixNDiff ! output ! do f = 1, ncmax call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_Adv', & & xyzf_QMixAdv(1:nx,1:ny,1:nz,f)) call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_NDiff', & & xyzf_QMixNDiff(1:nx,1:ny,1:nz,f)) end do ! VelX ! pyz_DVelXDtNl = pyz_DVelXDtNl + pyz_VelXAdv + pyz_VelXNDiff call HistoryAutoPut(TimeN, 'VelXAdv', pyz_VelXAdv(1:nx,1:ny,1:nz) ) call HistoryAutoPut(TimeN, 'VelXNDiff', pyz_VelXnDiff(1:nx,1:ny,1:nz) ) ! VelY ! xqz_DVelYDtNl = xqz_DVelYDtNl + xqz_VelYAdv + xqz_VelYNDiff call HistoryAutoPut(TimeN, 'VelYAdv', xqz_VelYAdv(1:nx,1:ny,1:nz) ) call HistoryAutoPut(TimeN, 'VelYNDiff', xqz_VelYnDiff(1:nx,1:ny,1:nz) ) ! VelZ ! call BuoyancyLong_xyr ! 浮力項 xyr_DVelZDtNl = & & xyr_DVelZDtNl & & + xyr_VelZAdv & & + xyr_VelZnDiff & & + xyr_BuoyT * FactorBuoyTemp & & + xyr_BuoyM * FactorBuoyMolWt & & + xyr_BuoyD * FactorBuoyLoading call HistoryAutoPut(TimeN, 'VelZAdv', xyr_VelZAdv(1:nx,1:ny,1:nz) ) call HistoryAutoPut(TimeN, 'VelZNDiff', xyr_VelZnDiff(1:nx,1:ny,1:nz) ) call HistoryAutoPut(TimeN, 'VelZBuoyT', xyr_BuoyT(1:nx,1:ny,1:nz) ) call HistoryAutoPut(TimeN, 'VelZBuoyM', xyr_BuoyM(1:nx,1:ny,1:nz) ) call HistoryAutoPut(TimeN, 'VelZBuoyD', xyr_BuoyD(1:nx,1:ny,1:nz) ) contains subroutine BuoyancyLong_xyr ! ! 浮力項の計算 ! use composition, only: GasNum, &! & IdxG, &! & MolWtWet ! 湿潤成分の分子量 use constants, only: MolWtDry, &! 乾燥成分の分子量 & Grav ! 重力加速度 use basicset, only : xyr_PTempBZ, &! 基本場の温位 & xyr_QMixBZ, &! 基本場の混合比 & xyr_QMixBZPerMolWt ! 暗黙の型宣言禁止 ! implicit none ! 作業変数 ! real(DP) :: xyzf_QMixPerMolWt(imin:imax,jmin:jmax,kmin:kmax, 1:GasNum) real(DP) :: xyz_QMixPerMolWtSum(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_QMixNlSumGas(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_QMixNlSumCnd(imin:imax,jmin:jmax,kmin:kmax) integer :: i, j, k, f, n ! Buoyancy due to temperature disturbunce ! do k = kmin, kmax-1 do j = jmin, jmax do i = imin, imax xyr_BuoyT(i,j,k) = & & Grav & & * ( & & xyz_PTempNl(i,j,k+1) & & + xyz_PTempNl(i,j,k) & & ) * 5.0d-1 & & / xyr_PTempBZ(i,j,k) end do end do end do ! 穴埋め ! xyr_BuoyT(:,:,kmax) = 0.0d0 if ( .NOT. FlagCalcMoist ) then ! 乾燥の場合は BuoyD, BuoyM は零. ! xyr_BuoyM = 0.0d0 xyr_BuoyD = 00d0 ! サブルーチンを抜ける ! return end if ! Buoyancy due to molecular weight ! do f = 1, GasNum n = IdxG(f) xyzf_QMixPerMolWt(:,:,:,f) = xyzf_QMixNl(:,:,:,n) / MolWtWet(n) end do xyz_QMixPerMolWtSum = sum(xyzf_QMixPerMolWt, 4) xyz_QMixNlSumGas = sum(xyzf_QMixNl(:,:,:,1:GasNum), 4) do k = kmin, kmax-1 do j = jmin, jmax do i = imin, imax xyr_BuoyM(i,j,k) = & & + Grav & & * ( & & xyz_QMixPerMolWtSum(i,j,k+1) & & + xyz_QMixPerMolWtSum(i,j,k) & & ) * 5.0d-1 & & / ( 1.0d0 / MolWtDry + xyr_QMixBZPerMolWt(i,j,k) ) & & - Grav & & * ( & & xyz_QMixNlSumGas(i,j,k+1) & & + xyz_QMixNlSumGas(i,j,k) & & ) * 5.0d-1 & & / ( 1.0d0 + xyr_QmixBZ(i,j,k) ) end do end do end do ! 穴埋め ! xyr_BuoyM(:,:,kmax) = 0.0d0 ! Buoyancy due to loading ! xyz_QMixNlSumCnd = sum( xyzf_QMixNl(:,:,:,GasNum+1:ncmax), 4 ) do k = kmin, kmax-1 do j = jmin, jmax do i = imin, imax xyr_BuoyD(i,j,k) = & & - Grav & & * ( & & xyz_QMixNlSumCnd(i,j,k+1) & & + xyz_QMixNlSumCnd(i,j,k) & & ) * 5.0d-1 & & / ( 1.0d0 + xyr_QMixBZ(i,j,k) ) end do end do end do ! 穴埋め ! xyr_BuoyD(:,:,kmax) = 0.0d0 end subroutine BuoyancyLong_xyr end subroutine Dynamics_Long_forcing !!!-------------------------------------------------------------!!! subroutine Dynamics_Long_integrate( & & xyz_VarBl, xyz_DVarDtNl, & !(IN) & xyzf_VarBl, xyzf_DVarDtNl, & !(IN) & xyz_VarAl, & !(OUT) & xyzf_VarAl & !(OUT) ) ! モジュール読み込み ! use dc_types, only : DP use gridset, only : imin, &! x 方向の配列の下限 & imax, &! x 方向の配列の上限 & jmin, &! y 方向の配列の下限 & jmax, &! y 方向の配列の上限 & kmin, &! z 方向の配列の下限 & kmax, &! z 方向の配列の上限 & ncmax ! 物質数 use timeset, only : DelTimeLong use setmargin,only : SetMargin_xyz, & & SetMargin_xyzf use fillnegative, & & only : FillNegativeQMix ! 暗黙の型宣言禁止 ! implicit none ! 変数の定義 ! real(DP), intent(in) :: xyz_VarBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_DVarDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyzf_VarBl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) real(DP), intent(in) :: xyzf_DVarDtNl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) real(DP), intent(out) :: xyz_VarAl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyzf_VarAl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) !------------------------------------------------------------------ ! スカラー量の積分 ! Integration ! xyz_VarAl = xyz_VarBl + (2.0d0 * DelTimeLong) * xyz_DVarDtNl ! Set Margin ! call SetMargin_xyz(xyz_VarAl) !------------------------------------------------------------------ ! 混合比の積分 ! Integration ! xyzf_VarAl = xyzf_VarBl + (2.0d0 * DelTimeLong) * xyzf_DVarDtNl call SetMargin_xyzf(xyzf_VarAl) ! QMix については, 周囲を削って負の値を埋める ! 穴埋めするためには, 先に setmargin しておく必要がある. ! call FillNegativeQMix(xyzf_VarAl) ! のり代の値の設定. Set Margin ! call SetMargin_xyzf(xyzf_VarAl) end subroutine Dynamics_Long_integrate !!!---------------------------------------------------------------!!! subroutine Dynamics_Short_forcing( & & pyz_VelXNs, & ! (in) & xqz_VelYNs, & ! (in) & xyr_VelZNs, & ! (in) & xyz_ExnerNs, & ! (in) & pyz_DVelXDtNl, & ! (in) & xqz_DVelYDtNl, & ! (in) & xyr_DVelZDtNl, & ! (in) & xyz_DExnerDtNl, & ! (in) & xyz_DExnerDtNs, & ! (in) & pyz_VelXAs, & ! (out) & xqz_VelYAs, & ! (out) & xyr_VelZAs, & ! (out) & xyz_ExnerAs & ! (out) & ) ! モジュール読み込み ! use dc_types, only : DP use gridset, only : imin, &! x 方向の配列の下限 & imax, &! x 方向の配列の上限 & jmin, &! y 方向の配列の下限 & jmax, &! y 方向の配列の上限 & kmin, &! z 方向の配列の下限 & kmax, &! z 方向の配列の上限 & nx, ny, nz ! 物理領域の大きさ use timeset, only : DelTimeShort, TimeN use setmargin,only : SetMargin_xyzf, SetMargin_xyz, & & SetMargin_pyz, SetMargin_xqz, SetMargin_xyr use fillnegative, & & only : FillNegativeQMix use gtool_historyauto, & & only : HistoryAutoPut use acousticmode_std, & & only : acousticmode_std_exp, acousticmode_std_imp use acousticmode_2D, & & only : acousticmode_2D_exp, acousticmode_2D_imp use acousticmode_3D, & & only : acousticmode_3D_exp, acousticmode_3D_imp use constants,only : CvDry, &! 乾燥成分の比熱 & GasRDry ! 暗黙の型宣言禁止 ! implicit none ! 変数の定義 ! real(DP), intent(in) :: pyz_VelXNs(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xqz_VelYNs(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_VelZNs(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_ExnerNs(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: pyz_DVelXDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xqz_DVelYDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_DExnerDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: xyz_DExnerDtNs(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: pyz_VelXAs(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xqz_VelYAs(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyr_VelZAs(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyz_ExnerAs(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_DVelXDtNs(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_DVelYDtNs(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_DVelZDtNs(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_VelDivNs(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_Expnd(imin:imax,jmin:jmax,kmin:kmax) ! subroutine の引数となる配列は, 物理領域の大きさに. ! real(DP) :: pyz_PGrad(1:nx,1:ny,1:nz) real(DP) :: xqz_PGrad(1:nx,1:ny,1:nz) real(DP) :: xyr_PGrad(1:nx,1:ny,1:nz) real(DP) :: pyz_SWF(1:nx,1:ny,1:nz) real(DP) :: xqz_SWF(1:nx,1:ny,1:nz) real(DP) :: xyr_SWF(1:nx,1:ny,1:nz) !------------------------------------------------------------ ! 初期化 ! pyz_DVelXDtNs = 0.0d0 xqz_DVelYDtNs = 0.0d0 xyr_DVelZDtNs = 0.0d0 !------------------------------------------------------------ ! 水平方向: 陽解法 ! select case ( IDAcousticmode ) case ( IDAcousticmode_std ) call acousticmode_std_exp( & & pyz_VelXNs, xqz_VelYNs, xyr_VelZNs, & !(IN) & xyz_ExnerNs, & !(IN) & xyz_VelDivNs, & !(OUT) & pyz_PGrad, pyz_SWF, & !(OUT) & xqz_PGrad, xqz_SWF & !(OUT) & ) case ( IDAcousticmode_3D ) call acousticmode_3d_exp( & & pyz_VelXNs, xqz_VelYNs, xyr_VelZNs, & !(IN) & xyz_ExnerNs, & !(IN) & xyz_VelDivNs, & !(OUT) & pyz_PGrad, pyz_SWF, & !(OUT) & xqz_PGrad, xqz_SWF & !(OUT) & ) case ( IDAcousticmode_2D ) xqz_PGrad = 0.0d0 xqz_SWF = 0.0d0 call acousticmode_2d_exp( & & pyz_VelXNs, xyr_VelZNs, & !(IN) & xyz_ExnerNs, & !(IN) & xyz_VelDivNs, & !(OUT) & pyz_PGrad, pyz_SWF & !(OUT) & ) end select ! tendency ! pyz_DVelXDtNs(1:nx,1:ny,1:nz) = pyz_PGrad + pyz_SWF xqz_DVelYDtNs(1:nx,1:ny,1:nz) = xqz_PGrad + xqz_SWF ! 値の保管 ! call HistoryAutoPut(TimeN, 'VelXPGrad', pyz_PGrad ) call HistoryAutoPut(TimeN, 'VelXSWF', pyz_SWF ) call HistoryAutoPut(TimeN, 'VelYPGrad', xqz_PGrad ) call HistoryAutoPut(TimeN, 'VelYSWF', xqz_SWF ) ! Time integration ! pyz_VelXAs = pyz_VelXNs + DelTimeShort * ( pyz_DVelXDtNl + pyz_DVelXDtNs ) xqz_VelYAs = xqz_VelYNs + DelTimeShort * ( xqz_DVelYDtNl + xqz_DVelYDtNs ) ! Set Margin ! call SetMargin_pyz( pyz_VelXAs ) ! (inout) call SetMargin_xqz( xqz_VelYAs ) ! (inout) !------------------------------------------------------------ ! 鉛直方向: 陰解法 ! 凝結項を短い時間ステップで評価する場合を考えて, 長い時間ステップと ! 短い時間ステップで評価された tendency を足し合わせたものを引数に. ! ! 短い時間ステップで評価する圧力の式の tendency ! xyz_Expnd = GasRDry * xyz_ExnerNs * xyz_VelDivNs / CvDry & & * FactorExnerExpnd ! tendency の合計 ! xyz_DExnerDtNs = xyz_DExnerDtNs + xyz_Expnd ! 値の保管 ! call HistoryAutoPut(TimeN, 'ExnerExpnd', xyz_Expnd(1:nx,1:ny,1:nz)) select case ( IDAcousticmode ) case ( IDAcousticmode_std ) call acousticmode_std_imp( & & pyz_VelXAs, xqz_VelYAs, xyr_VelZNs, xyz_VelDivNs, & !(IN) & xyz_ExnerNs, & !(IN) & xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs, & !(IN) & xyz_ExnerAs, & !(OUT) & xyr_PGrad, xyr_SWF & !(OUT) & ) case ( IDAcousticmode_3D ) call acousticmode_3D_imp( & & pyz_VelXAs, xqz_VelYAs, xyr_VelZNs, xyz_VelDivNs, & !(IN) & xyz_ExnerNs, & !(IN) & xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs, & !(IN) & xyz_ExnerAs, & !(OUT) & xyr_PGrad, xyr_SWF & !(OUT) & ) case ( IDAcousticmode_2D ) call acousticmode_2D_imp( & & pyz_VelXAs, xyr_VelZNs, xyz_VelDivNs, & !(IN) & xyz_ExnerNs, & !(IN) & xyr_DVelZDtNl, xyz_DExnerDtNl, xyz_DExnerDtNs, & !(IN) & xyz_ExnerAs, & !(OUT) & xyr_PGrad, xyr_SWF & !(OUT) & ) end select ! tendency ! xyr_DVelZDtNs(1:nx,1:ny,1:nz) = xyr_PGrad + xyr_SWF ! Time integration ! xyr_VelZAs = xyr_VelZNs + DelTimeShort * ( xyr_DVelZDtNl + xyr_DVelZDtNs ) ! Set Margin ! call SetMargin_xyz( xyz_ExnerAs ) ! (inout) call SetMargin_xyr( xyr_VelZAs ) ! (inout) ! 値の保管 ! call HistoryAutoPut(TimeN, 'VelZPGrad', xyr_PGrad ) call HistoryAutoPut(TimeN, 'VelZSWF', xyr_SWF ) end subroutine Dynamics_Short_forcing !!!------------------------------------------------------------------------!!! subroutine Dynamics_Tendency_Output use gtool_historyauto, only : HistoryAutoAddVariable use gridset, only : ncmax ! 物質数 use composition, only : SpcWetSymbol implicit none integer :: f call HistoryAutoAddVariable( & & varname='PTempAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of potential temperature', & & units='K.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='PTempNDiff', & & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of potential temperature', & & units='K.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='ExnerAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of exner function', & & units='s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='ExnerNDiff', & & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of exner function', & & units='s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='ExnerExpnd', & & dims=(/'x','y','z','t'/), & & longname='Expanding term of exner function', & & units='s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='CDensAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of cloud density', & & units='kg.m-3.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='CDensNDiff', & & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of cloud density', & & units='kg.m-3.s-1', & & xtype='float') do f = 1, ncmax call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(f))//'_Adv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of ' & & //trim(SpcWetSymbol(f))//' mixing ratio', & & units='kg.kg-1.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(f))//'_NDiff', & & dims=(/'x','y','z','t'/), & & longname='Diffusion term of ' & & //trim(SpcWetSymbol(f))//' mixing ratio', & & units='kg.kg-1.s-1', & & xtype='float') end do call HistoryAutoAddVariable( & & varname='VelXAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of velocity (x)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelXNDiff', & & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of velocity (x)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelXPGrad', & & dims=(/'x','y','z','t'/), & & longname='Pressure gradient term of velocity (x)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelXSWF', & & dims=(/'x','y','z','t'/), & & longname='Filter for acoustic mode (x)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelYAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of velocity (y)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelYNDiff', & & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of velocity (y)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelYPGrad', & & dims=(/'x','y','z','t'/), & & longname='Pressure gradient term of velocity (y)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelYSWF', & & dims=(/'x','y','z','t'/), & & longname='Filter for acoustic mode (y)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelZAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of velocity (z)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelZNDiff', & & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of Velocity (z)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelZBuoyT', & & dims=(/'x','y','z','t'/), & & longname='Buoyancy (Temperature)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelZBuoyM', & & dims=(/'x','y','z','t'/), & & longname='Buoyancy (MolWt)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelZBuoyD', & & dims=(/'x','y','z','t'/), & & longname='Buoyancy (Drag)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelZPGrad', & & dims=(/'x','y','z','t'/), & & longname='Pressure gradient term of velocity (z)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='VelZSWF', & & dims=(/'x','y','z','t'/), & & longname='Filter for acoustic mode (z)', & & units='m.s-2', & & xtype='float') call HistoryAutoAddVariable( & & varname='KmAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection of Km', & & units='s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname='KmNDiff', & & dims=(/'x','y','z','t'/), & & longname='Diffusion term of Km', & & units='s-1', & & xtype='float') end subroutine Dynamics_Tendency_Output end module DynamicsHEVI