!= Module DynamicsHEVI ! ! Authors:: 杉山耕一朗(SUGIYAMA Ko-ichiro), ODAKA Masatsugu ! Version:: $Id: dynamics_hevi_v2.test.f90,v 1.1 2013-06-15 15:31:41 sugiyama Exp $ ! Tag Name:: $Name: $ ! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! !モデルの力学過程を計算するために必要となる関数群を束ねたモジュール !具体的には以下の項を計算するための関数を格納する. ! * 移流項 ! * 浮力項 ! * 気圧傾度力項 ! !== Error Handling ! !== Known Bugs ! !== Note ! ! * エクスナー関数の空間方向の離散化において, 2 次精度の離散化を陽に利 ! 用しているため, 気圧傾度力項の計算プログラムにおいて ! differentiate_center4 モジュールを指定することはできないので注意. ! !== Future Plans ! module DynamicsHEVI ! !陽解法を用いた力学過程の各項の計算モジュール. !具体的には以下の項を計算するための関数を格納する. ! * 移流項 ! * 浮力項 ! * 気圧傾度力項 ! !モジュール読み込み use dc_types, only : DP, STRING use dc_iounit, only : FileOpen use dc_message, only : MessageNotify use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut use axesset use xyz_deriv_module use xyz_deriv_c4_module use mpi_wrapper,only: myrank use gridset, only: & & imin, &! x 方向の配列の下限 & imax, &! x 方向の配列の上限 & jmin, &! y 方向の配列の下限 & jmax, &! y 方向の配列の上限 & kmin, &! z 方向の配列の下限 & kmax, &! z 方向の配列の上限 & nx, &! x 方向の物理領域の上限 & ny, &! x 方向の物理領域の上限 & nz, &! y 方向の物理領域の上限 & ncmax, &! 物質数 & FlagCalc3D use constants,only: CpDry ! 乾燥成分の比熱 use composition, only: SpcWetSymbol use timeset, only: DelTimeShort, DelTimeLong, TimeN ! use axesset, only: dx, dy, dz ! 格子間隔 use basicset, only: xyz_VelSW, &!基本場の音速 & xyz_VelSoundBZ, &!基本場の音速 & xyz_DensBZ, &!基本場の密度 & xyz_PTempBZ, &!基本場の温位 & xyz_VPTempBZ, &!基本場の温位 & pyz_VPTempBZ, &!基本場の温位 & xqz_VPTempBZ, &!基本場の温位 & xyr_VPTempBZ, &!基本場の温位 & xyz_ExnerBZ, & & xyzf_QMixBZ, & & xyr_QMixBZPerMolWt, & & xyr_QMixBZ, xyz_EffMolWtBZ use setmargin,only: SetMargin_xyzf, SetMargin_xyz, & & SetMargin_pyz, SetMargin_xqz, SetMargin_xyr use fillnegative, only: FillNegativeQMix use namelist_util, only: namelist_filename !暗黙の型宣言禁止 implicit none !属性の指定 private real(DP), save :: beta = 1.0d0 !クランクニコルソン法なら 0.5 !完全陰解法なら 1 integer, save :: N = 10 !係数行列/改行列の次数, 整合寸法 integer, save :: M = 10 !方程式の組数 integer, save :: NUD = 1 !係数行列の上三角部分の帯幅 integer, save :: NLD = 1 !係数行列の下三角部分の帯幅 integer, save :: NAL = 1 !LU 分解の結果 L の整合寸法 integer, save :: NA = 3 !NUD + NLD + 1 real(DP), allocatable, save :: xyz_F3BZ(:,:,:) !係数行列の計算に用いる配列 real(DP), allocatable, save :: xyz_F1BZ(:,:,:) !係数行列の計算に用いる配列 real(DP), allocatable, save :: xyr_F2BZ(:,:,:) !係数行列の計算に用いる配列 real(DP), allocatable, save :: xyr_CpVPTempBZ(:,:,:) !係数行列の計算に用いる配列 real(DP), allocatable, save :: xyr_CpDensVPTemp2BZ(:,:,:) !係数行列の計算に用いる配列 real(DP), allocatable, save :: xyr_DensVPTempBZ(:,:,:) !係数行列の計算に用いる配列 real(DP), allocatable, save :: A(:) !係数行列の対角成分 real(DP), allocatable, save :: B(:) !係数行列の上三角部分 real(DP), allocatable, save :: C(:) !係数行列の下三角部分 real(DP), allocatable, save :: AL1(:) !LU 分解の結果 L (1 次元配列) integer, allocatable, save :: IP(:) !部分ピボット交換の情報を格納 real(DP), save :: Alpha = 0.0d0 !音波減衰項の減衰係数 real(DP), save :: NuHh = 0.0d0 !数値粘性の係数 (水平方向) real(DP), save :: NuVh = 0.0d0 !数値粘性の係数 (鉛直方向) real(DP), save :: NuHm = 0.0d0 !数値粘性の係数 (水平方向) real(DP), save :: NuVm = 0.0d0 !数値粘性の係数 (鉛直方向) character(*), parameter:: module_name = 'DynamicHEVI' ! モジュールの名称. ! Module name real(DP), save :: FactorBuoyTemp = 1.0d0 !浮力 (温度の寄与) の有無 !考慮しない場合は値をゼロにする. real(DP), save :: FactorBuoyMolWt = 1.0d0 !浮力 (分子量効果) の有無 !考慮しない場合は値をゼロにする. real(DP), save :: FactorBuoyLoading = 1.0d0 !浮力 (荷重効果) の有無 !考慮しない場合は値をゼロにする. !public public Dynamics_Init public Dynamics_Long_forcing public Dynamics_Short_forcing public Dynamics_Km_forcing contains subroutine Dynamics_Init !暗黙の型宣言禁止 implicit none real(DP) :: DelXMin, DelYMin, DelZMin real(DP) :: AlphaSound = 5.0d-2 !音波減衰項の係数 (気象庁数値予報課報告・別冊49 より) real(DP) :: AlphaNDiff = 1.0d-3 !4次の数値拡散の係数. CReSS マニュアルより real(DP) :: NDiffRatio = 1.0d0 !速度に対する粘性を上げる場合は数字を 1 以上にする. integer :: unit !装置番号 NAMELIST /Dynamics_nml/ & & AlphaSound, AlphaNDiff, NDiffRatio, beta, & & FactorBuoyTemp, FactorBuoyMolWt, FactorBuoyLoading !ファイルオープン. 情報取得. call FileOpen(unit, file=namelist_filename, mode='r') read(unit, NML=dynamics_nml) close(unit) !------------------------------------------------------------------- ! 音波減衰項の減衰率 Min(DelX, DelZ) ** 2.0 に比例 ! 2 次元計算の場合には DelY に依存しないようにするために if 文を利用. ! DelXMin = dx DelYMin = dy DelZMin = dz if ( FlagCalc3D ) then Alpha = AlphaSound * ( Min(DelXMin, DelYMin, DelZMin) ** 2.0d0 ) / DelTimeShort else Alpha = AlphaSound * ( Min(DelXMin, DelZMin) ** 2.0d0 ) / DelTimeShort end if if ( FlagCalc3D ) then NuHh = AlphaNDiff * ( SQRT(DelXMin*DelYMin) ** 4.0d0 ) / (2.0d0 * DelTimeLong) else NuHh = AlphaNDiff * ( DelXMin ** 4.0d0 ) / (2.0d0 * DelTimeLong) end if NuVh = AlphaNDiff * ( DelZMin ** 4.0d0 ) / (2.0d0 * DelTimeLong) NuHm = NuHh * NDiffRatio NuVm = NuVh * NDiffRatio if (myrank == 0) then call MessageNotify( "M", module_name, "Alpha = %f", d=(/Alpha/) ) call MessageNotify( "M", module_name, "NuHh = %f", d=(/NuHh/) ) call MessageNotify( "M", module_name, "NuVh = %f", d=(/NuVh/) ) call MessageNotify( "M", module_name, "NuHm = %f", d=(/NuHm/) ) call MessageNotify( "M", module_name, "NuVm = %f", d=(/NuVm/) ) 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/) ) end if ! 鉛直陰解法を用いるために, 行列の準備を行う. ! call Dynamics_VI_init ! tendency の出力 ! call Dynamics_Tendency_output 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) & xyzf_QMixBl, xyzf_QMixNl, & ! (in) & pyz_DVelXDtNl, & ! (inout) & xqz_DVelYDtNl, & ! (inout) & xyr_DVelZDtNl, & ! (inout) & xyz_DPTempDtNl, & ! (inout) & xyzf_DQMixDtNl, & ! (inout) & xyz_PTempAl, & ! (out) & xyzf_QMixAl & ! (out) & ) 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) :: 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(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) :: xyzf_DQMixDtNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP), intent(out) :: xyz_PTempAl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyzf_QMixAl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP) :: pyz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_nDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_nDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_nDiff(imin:imax,jmin:jmax,kmin:kmax) 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) real(DP) :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyzf_Adv(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP) :: xyzf_nDiff(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP) :: xyz_PTempAll(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) integer :: f real(DP) :: xyz_Adv2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyzf_Adv2(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) real(DP) :: xyzf_nDiff2(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) real(DP) :: pyz_Adv2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_Adv2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_nDiff2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_Adv2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_nDiff2(imin:imax,jmin:jmax,kmin:kmax) character(15) :: info integer :: k !-------------------------------------------------------------------- ! 温位 ! 移流については, 基本場も考慮する. ! xyz_PTempAll = xyz_PTempNl + xyz_PTempBZ call AdvC4_nDiff_xyz( & & xyz_PTempBl, xyz_PTempAll, & !(IN) & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN) & xyz_Adv, xyz_nDiff ) !(OUT) xyz_Adv2 = & & - xyz_avr_pyz(pyz_VelXNl * pyz_c4dx_xyz(xyz_PTempAll)) & & - xyz_avr_xqz(xqz_VelYNl * xqz_c4dy_xyz(xyz_PTempAll)) & & - xyz_avr_xyr(xyr_VelZNl * xyr_c4dz_xyz(xyz_PTempAll)) xyz_nDiff2 = & & - NuHh * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyz_PTempBl ))))) & & - NuHh * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_PTempBl ))))) & & - NuVh * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyz_PTempBl ))))) info = "PTemp, Adv: " call CHK_Val( xyz_Adv(imin:imax,jmin:jmax,kmin:kmax), xyz_Adv2(imin:imax,jmin:jmax,kmin:kmax), info ) info = "PTemp, nDiff: " call CHK_Val( xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax), xyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax), info ) ! tendency の更新 ! xyz_DPTempDtNl = xyz_DPTempDtNl + xyz_Adv + xyz_nDiff ! output ! call HistoryAutoPut(TimeN, 'PTempAdv', xyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'PTempNDiff', xyz_nDiff(1:nx,1:ny,1:nz)) !-------------------------------------------------------------------- ! 混合比 ! ! 移流については, 基本場も考慮する. ! xyzf_QMixAll = xyzf_QMixNl + xyzf_QMixBZ call AdvC4_nDiff_xyzf( & & xyzf_QMixBl, xyzf_QMixAll, & !(IN) & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN) & xyzf_Adv, xyzf_nDiff ) !(OUT) do f = 1, ncmax xyzf_Adv2(:,:,:,f) = & & - xyz_avr_pyz(pyz_VelXNl * pyz_c4dx_xyz(xyzf_QMixAll(:,:,:,f))) & & - xyz_avr_xqz(xqz_VelYNl * xqz_c4dy_xyz(xyzf_QMixAll(:,:,:,f))) & & - xyz_avr_xyr(xyr_VelZNl * xyr_c4dz_xyz(xyzf_QMixAll(:,:,:,f))) xyzf_nDiff2(:,:,:,f) = & & - NuHh * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyzf_QMixBl(:,:,:,f) ))))) & & - NuHh * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyzf_QMixBl(:,:,:,f) ))))) & & - NuVh * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyzf_QMixBl(:,:,:,f) ))))) end do do f = 1, ncmax info = "QMix, Adv: " call CHK_Val( xyzf_Adv(imin:imax,jmin:jmax,kmin:kmax,f), xyzf_Adv2(imin:imax,jmin:jmax,kmin:kmax,f), info ) info = "QMix, nDiff: " call CHK_Val( xyzf_nDiff(imin:imax,jmin:jmax,kmin:kmax,f), xyzf_nDiff2(imin:imax,jmin:jmax,kmin:kmax,f), info ) end do ! tendency の更新 ! xyzf_DQMixDtNl = xyzf_DQMixDtNl + xyzf_Adv + xyzf_nDiff ! output ! do f = 1, ncmax call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_Adv', xyzf_Adv(1:nx,1:ny,1:nz,f)) call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_NDiff', xyzf_NDiff(1:nx,1:ny,1:nz,f)) end do !------------------------------------------------------------------ ! VelX, VelY, VelZ ! ! 移流項・数値拡散項をまとめて計算 ! call AdvC4_nDiff_pyz_xqz_xyr( & & pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, & !(IN) & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN) & pyz_Adv, xqz_Adv, xyr_Adv, & !(OUT) & pyz_nDiff, xqz_nDiff, xyr_nDiff ) !(OUT) pyz_Adv2 = & & - pyz_VelXNl * pyz_avr_xyz( xyz_c4dx_pyz( pyz_VelXNl ) ) & & - pyz_avr_pqz( pqz_avr_xqz( xqz_VelYNl ) * pqz_c4dy_pyz( pyz_VelXNl ) ) & & - pyz_avr_pyr( pyr_avr_xyr( xyr_VelZNl ) * pyr_c4dz_pyz( pyz_VelXNl ) ) pyz_nDiff2 = & & - NuHm * ( pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz( pyz_VelXBl ))))) & & - NuHm * ( pyz_dy_pqz(pqz_dy_pyz(pyz_dy_pqz(pqz_dy_pyz( pyz_VelXBl ))))) & & - NuVm * ( pyz_dz_pyr(pyr_dz_pyz(pyz_dz_pyr(pyr_dz_pyz( pyz_VelXBl ))))) xqz_Adv2 = & & - xqz_avr_pqz( pqz_avr_pyz( pyz_VelXNl ) * pqz_c4dx_xqz( xqz_VelYNl ) ) & & - xqz_VelYNl * xqz_avr_xyz( xyz_c4dy_xqz( xqz_VelYNl ) ) & & - xqz_avr_xqr( xqr_avr_xyr( xyr_VelZNl ) * xqr_c4dz_xqz( xqz_VelYNl ) ) xqz_nDiff2 = & & - NuHm * ( xqz_dx_pqz(pqz_dx_xqz(xqz_dx_pqz(pqz_dx_xqz( xqz_VelYBl ))))) & & - NuHm * ( xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz( xqz_VelYBl ))))) & & - NuVm * ( xqz_dz_xqr(xqr_dz_xqz(xqz_dz_xqr(xqr_dz_xqz( xqz_VelYBl ))))) xyr_Adv2 = & & - xyr_avr_pyr( pyr_avr_pyz( pyz_VelXNl ) * pyr_c4dx_xyr( xyr_VelZNl ) ) & & - xyr_avr_xqr( xqr_avr_xqz( xqz_VelYNl ) * xqr_c4dy_xyr( xyr_VelZNl ) ) & & - xyr_VelZNl * xyr_avr_xyz( xyz_c4dz_xyr( xyr_VelZNl ) ) xyr_nDiff2 = & & - NuHm * ( xyr_dx_pyr(pyr_dx_xyr(xyr_dx_pyr(pyr_dx_xyr( xyr_VelZBl ))))) & & - NuHm * ( xyr_dy_xqr(xqr_dy_xyr(xyr_dy_xqr(xqr_dy_xyr( xyr_VelZBl ))))) & & - NuVm * ( xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr( xyr_VelZBl ))))) info = "VelX, Adv: " call CHK_Val( pyz_Adv(imin:imax,jmin:jmax,kmin:kmax), pyz_Adv2(imin:imax,jmin:jmax,kmin:kmax), info ) info = "VelX, nDiff: " call CHK_Val( pyz_nDiff(imin:imax,jmin:jmax,kmin:kmax), pyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax), info ) info = "VelY, Adv: " call CHK_Val( xqz_Adv(imin:imax,jmin:jmax,kmin:kmax), xqz_Adv2(imin:imax,jmin:jmax,kmin:kmax), info ) info = "VelY, nDiff: " call CHK_Val( xqz_nDiff(imin:imax,jmin:jmax,kmin:kmax), xqz_nDiff2(imin:imax,jmin:jmax,kmin:kmax), info ) info = "VelZ, Adv: " call CHK_Val( xyr_Adv(imin:imax,jmin:jmax,kmin:kmax), xyr_Adv2(imin:imax,jmin:jmax,kmin:kmax), info ) info = "VelZ, nDiff: " call CHK_Val( xyr_nDiff(imin:imax,jmin:jmax,kmin:kmax), xyr_nDiff2(imin:imax,jmin:jmax,kmin:kmax), info ) ! tendency of VelX ! pyz_DVelXDtNl = pyz_DVelXDtNl + pyz_Adv + pyz_NDiff call HistoryAutoPut(TimeN, 'VelXAdv', pyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelXNDiff', pyz_nDiff(1:nx,1:ny,1:nz)) ! tendency of VelY ! xqz_DVelYDtNl = xqz_DVelYDtNl + xqz_Adv + xqz_NDiff call HistoryAutoPut(TimeN, 'VelYAdv', xqz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelYNDiff', xqz_nDiff(1:nx,1:ny,1:nz)) ! Buoyancy ! call BuoyancyLong_xyr( & & xyz_PTempNl, xyzf_QMixNl, & !(IN) & xyr_BuoyT, xyr_BuoyM, xyr_BuoyD ) !(OUT) ! tendency of VelZ ! xyr_DVelZDtNl = xyr_DVelZDtNl + xyr_Adv + xyr_NDiff & & + xyr_BuoyT * FactorBuoyTemp & & + xyr_BuoyM * FactorBuoyMolWt & & + xyr_BuoyD * FactorBuoyLoading call HistoryAutoPut(TimeN, 'VelZAdv', xyr_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelZNDiff', xyr_NDiff(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)) !------------------------------------------------------------------ ! 長い時間ステップでの積分 ! Integration ! xyz_PTempAl = xyz_PTempBl + (2.0d0 * DelTimeLong) * xyz_DPTempDtNl xyzf_QMixAl = xyzf_QMixBl + (2.0d0 * DelTimeLong) * xyzf_DQMixDtNl ! Set Margin ! call SetMargin_xyz(xyz_PTempAl) call SetMargin_xyzf(xyzf_QMixAl) ! QMix については, 周囲を削って負の値を埋める ! 穴埋めするためには, 先に setmargin しておく必要がある. ! call FillNegativeQMix(xyzf_QMixAl) ! のり代の値の設定. Set Margin ! call SetMargin_xyzf(xyzf_QMixAl) end subroutine Dynamics_Long_forcing subroutine Dynamics_Long_Mars_forcing( & & pyz_VelXBl, pyz_VelXNl, & ! (in) & xqz_VelYBl, xqz_VelYNl, & ! (in) & xyr_VelZBl, xyr_VelZNl, & ! (in) & xyz_PTempBl, xyz_PTempNl, & ! (in) & xyz_CDensBl, xyz_CDensNl, & ! (in) & pyz_DVelXDtNl, & ! (inout) & xqz_DVelYDtNl, & ! (inout) & xyr_DVelZDtNl, & ! (inout) & xyz_DPTempDtNl, & ! (inout) & xyz_DCDensDtNl, & ! (inout) & xyz_PTempAl & ! (out) & ) 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_CDensBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_CDensNl(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_DCDensDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out), optional :: xyz_PTempAl(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_nDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_nDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_nDiff(imin:imax,jmin:jmax,kmin:kmax) 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) real(DP) :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_PTempAll(imin:imax,jmin:jmax,kmin:kmax) !-------------------------------------------------------------------- ! 温位 ! 移流については, 基本場も考慮する. ! xyz_PTempAll = xyz_PTempNl + xyz_PTempBZ call AdvC4_nDiff_xyz( & & xyz_PTempBl, xyz_PTempAll, & !(IN) & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN) & xyz_Adv, xyz_nDiff ) !(OUT) ! tendency の更新 ! 基本場の移流の分を加えている. ! xyz_DPTempDtNl = xyz_DPTempDtNl + xyz_Adv + xyz_nDiff ! output ! call HistoryAutoPut(TimeN, 'PTempAdv', xyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'PTempNDiff', xyz_nDiff(1:nx,1:ny,1:nz)) !-------------------------------------------------------------------- ! 雲密度 ! call AdvC4_nDiff_xyz( & & xyz_CDensBl, xyz_CDensNl, & !(IN) & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN) & xyz_Adv, xyz_nDiff ) !(OUT) ! tendency の更新 ! xyz_DCDensDtNl = xyz_DCDensDtNl + xyz_Adv + xyz_nDiff ! 値の保管 ! call HistoryAutoPut(TimeN, 'CDensAdv', xyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'CDensNDiff', xyz_nDiff(1:nx,1:ny,1:nz)) !------------------------------------------------------------------ ! VelX, VelY, VelZ ! ! 移流項・数値拡散項をまとめて計算 ! call AdvC4_nDiff_pyz_xqz_xyr( & & pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, & !(IN) & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN) & pyz_Adv, xqz_Adv, xyr_Adv, & !(OUT) & pyz_nDiff, xqz_nDiff, xyr_nDiff ) !(OUT) ! tendency of VelX ! pyz_DVelXDtNl = pyz_DVelXDtNl + pyz_Adv + pyz_NDiff call HistoryAutoPut(TimeN, 'VelXAdv', pyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelXNDiff', pyz_nDiff(1:nx,1:ny,1:nz)) ! tendency of VelY ! xqz_DVelYDtNl = xqz_DVelYDtNl + xqz_Adv + xqz_NDiff call HistoryAutoPut(TimeN, 'VelYAdv', xqz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelYNDiff', xqz_nDiff(1:nx,1:ny,1:nz)) ! Buoyancy ! 火星版の浮力は? 要チェック. CHECK!! ! ! call BuoyancyLong_xyr( & ! & xyz_PTempNl, xyzf_QMixNl, & !(IN) ! & xyr_BuoyT, xyr_BuoyM, xyr_BuoyD ) !(OUT) xyr_BuoyT = 0.0d0; xyr_BuoyM = 0.0d0; xyr_BuoyD = 0.0d0 ! tendency of VelZ ! xyr_DVelZDtNl = xyr_DVelZDtNl + xyr_Adv + xyr_NDiff & & + xyr_BuoyT * FactorBuoyTemp & & + xyr_BuoyM * FactorBuoyMolWt & & + xyr_BuoyD * FactorBuoyLoading call HistoryAutoPut(TimeN, 'VelZAdv', xyr_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelZNDiff', xyr_NDiff(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)) !------------------------------------------------------------------ ! 長い時間ステップでの積分 ! Integration ! if (present(xyz_PTempAl)) then ! 積分 ! xyz_PTempAl = xyz_PTempBl + (2.0d0 * DelTimeLong) * xyz_DPTempDtNl ! Set Margin ! call SetMargin_xyz(xyz_PTempAl) end if end subroutine Dynamics_Long_Mars_forcing 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) & ) 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(inout) :: xyz_DExnerDtNl(imin:imax,jmin:jmax,kmin:kmax) !test real(DP), intent(inout) :: xyz_DExnerDtNs(imin:imax,jmin:jmax,kmin:kmax) !test 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) :: pyz_PGrad(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_PGrad(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_PGrad(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_SWF(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_SWF(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_SWF(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_PGrad2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_PGrad2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_PGrad2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_SWF2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_SWF2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_SWF2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_VelDivNs2(imin:imax,jmin:jmax,kmin:kmax) integer :: k character(15) :: info ! initialize: Divergence of velocity ! call VelDivC2( & & pyz_VelXNs, xqz_VelYNs, xyr_VelZNs, & !(IN) & xyz_VelDivNs ) !(OUT) xyz_VelDivNs2 = & & xyz_dx_pyz( pyz_VelXNs ) & & + xyz_dy_xqz( xqz_VelYNs ) & & + xyz_dz_xyr( xyr_VelZNs ) info = "Vel, Div: " call CHK_Val( xyz_VelDivNs(imin:imax,jmin:jmax,kmin:kmax), xyz_VelDivNs2(imin:imax,jmin:jmax,kmin:kmax), info ) !------------------------------------------------------------ ! VelX, VelY ! 水平方向は陽解法で解く. ! call PGrad_HE( & & xyz_VelDivNs, xyz_ExnerNs, & !(IN) & pyz_PGrad, pyz_SWF, & !(OUT) & xqz_PGrad, xqz_SWF ) !(OUT) pyz_PGrad2 = - pyz_avr_xyz( CpDry * xyz_VPTempBZ ) * pyz_dx_xyz( xyz_ExnerNs ) pyz_SWF2 = Alpha * pyz_dx_xyz( xyz_VelDivNs ) xqz_PGrad2 = - xqz_avr_xyz( CpDry * xyz_VPTempBZ ) * xqz_dy_xyz( xyz_ExnerNs ) xqz_SWF2 = Alpha * xqz_dy_xyz( xyz_VelDivNs ) info = "VelX, PGrad: " call CHK_Val( pyz_PGrad(imin:imax,jmin:jmax,kmin:kmax), pyz_PGrad2(imin:imax,jmin:jmax,kmin:kmax), info ) info = "VelX, SWF: " call CHK_Val( pyz_SWF(imin:imax,jmin:jmax,kmin:kmax), pyz_SWF2(imin:imax,jmin:jmax,kmin:kmax), info ) info = "VelY, PGrad: " call CHK_Val( xqz_PGrad(imin:imax,jmin:jmax,kmin:kmax), xqz_PGrad2(imin:imax,jmin:jmax,kmin:kmax), info ) info = "VelY, SWF: " call CHK_Val( xqz_SWF(imin:imax,jmin:jmax,kmin:kmax), xqz_SWF2(imin:imax,jmin:jmax,kmin:kmax), info ) ! tendency ! pyz_DVelXDtNs = pyz_PGrad + pyz_SWF xqz_DVelYDtNs = xqz_PGrad + xqz_SWF ! 値の保管 ! call HistoryAutoPut(TimeN, 'VelXPGrad', pyz_PGrad(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelXSWF', pyz_SWF(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelYPGrad', xqz_PGrad(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelYSWF', xqz_SWF(1:nx,1:ny,1:nz)) ! 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) !------------------------------------------------------------ ! Exner function ! 積分値を返すことに注意. ! call Exner_HEVI( & & pyz_VelXAs, xqz_VelYAs, xyr_VelZNs, & !(IN) & xyz_VelDivNs, xyz_ExnerNs, & !(IN) & xyr_DVelZDtNl, xyz_DExnerDtNl, & !(IN) & xyz_ExnerAs & !(OUT) & ) !------------------------------------------------------------ ! VelZ ! call PGrad_VI( & & Alpha, beta, dz, CpDry, xyr_VPTempBZ, & & imin, imax, jmin, jmax, kmin, kmax, & & xyz_VelDivNs, xyz_ExnerNs, xyz_ExnerAs, &!(IN) & xyr_PGrad, xyr_SWF ) !(OUT) xyr_SWF2 = Alpha * xyr_dz_xyz( xyz_VelDivNs ) xyr_PGrad2 = & & - xyr_avr_xyz(CpDry * xyz_VPTempBZ ) & & * ( & & beta * xyr_dz_xyz( xyz_ExnerAs ) & & + (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerNs ) & & ) info = "VelZ, PGrad: " call CHK_Val( xyr_PGrad(imin:imax,jmin:jmax,kmin:kmax), xyr_PGrad2(imin:imax,jmin:jmax,kmin:kmax), info ) info = "VelZ, SWF: " call CHK_Val( xyr_SWF(imin:imax,jmin:jmax,kmin:kmax), xyr_SWF2(imin:imax,jmin:jmax,kmin:kmax), info ) ! tendency ! xyr_DVelZDtNs = xyr_PGrad + xyr_SWF ! 値の保管 ! call HistoryAutoPut(TimeN, 'VelZPGrad', xyr_PGrad(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelZSWF', xyr_SWF(1:nx,1:ny,1:nz)) ! Time integration ! xyr_VelZAs = xyr_VelZNs + DelTimeShort * (xyr_DVelZDtNl + xyr_DVelZDtNs) ! Set Margin ! call SetMargin_xyr( xyr_VelZAs ) ! (inout) end subroutine Dynamics_Short_forcing !!!--------------------------------------------------------------------!!! subroutine Dynamics_VI_init() ! !エクスナー関数を陰解法で解く際に必要となる, 係数行列の要素を決め, !LU 分解を行う. ! !暗黙の型宣言禁止 implicit none real(DP) :: dt ! 短い時間格子 integer :: INFO !解のコンディションチェック integer :: i, j, k real(8) :: A2(1:nz) real(8) :: B2(2:nz) real(8) :: C2(1:nz-1) real(8) :: xyz_F1BZ2(imin:imax, jmin:jmax, kmin:kmax) real(8) :: xyr_F2BZ2(imin:imax, jmin:jmax, kmin:kmax) real(DP) :: dts ! 短い時間格子 character(15) :: info2 !---------------------------------------------------------------- ! 初期化 ! 変数名が長すぎたので, 名前を置き換える ! dt = DelTimeShort dts = DelTimeShort ! 配列の割り付け ! allocate( A(1:nz) ) allocate( B(2:nz) ) allocate( C(1:nz-1) ) allocate( xyz_F3BZ(imin:imax, jmin:jmax, kmin:kmax) ) allocate( xyz_F1BZ(1:nx,1:ny,1:nz) ) allocate( xyr_F2BZ(1:nx,1:ny,kmin:kmax) ) allocate( xyr_CpDensVPTemp2BZ(1:nx,1:ny,kmin:kmax) ) allocate( xyr_DensVPTempBZ(1:nx,1:ny,kmin:kmax) ) allocate( xyr_CpVPTempBZ(1:nx,1:ny,kmin:kmax) ) !---------------------------------------------------------------- ! 係数行列および共通して利用される配列の値を決める ! A, B, C を求める際, 基本場 (BZ) の量は X 方向に一様なので. ! nx, ny の値で代表させることとした. ! do k = 1, nz do j = 1, ny do i = 1, nx xyz_F1BZ(i,j,k) = & & ( xyz_VelSW(i,j,k) ** 2.0d0 ) & & / ( CpDry * xyz_DensBZ(i,j,k) * (xyz_VPTempBZ(i,j,k) ** 2.0d0) ) end do end do end do xyz_F1BZ2 = & & ( xyz_VelSoundBZ ** 2.0d0 ) & & / (CpDry * xyz_DensBZ * (xyz_VPTempBZ ** 2.0d0)) xyz_F3BZ = xyz_F1BZ2 info2 = "VI_init, F1BZ: " call CHK2_Val( xyz_F1BZ(1:nx,1:ny,1:nz), xyz_F1BZ2(1:nx,1:ny,1:nz), info2 ) do k = kmin, kmax-1 do j = 1, ny do i = 1, nx xyr_CpDensVPTemp2BZ(i,j,k)= & & CpDry & & * ( & & + xyz_DensBZ(i,j,k+1) & & * ( xyz_VPTempBZ(i,j,k+1) ** 2.0d0 ) & & + xyz_DensBZ(i,j,k) & & * ( xyz_VPTempBZ(i,j,k) ** 2.0d0 ) & & ) * 5.0d-1 end do end do end do xyr_CpDensVPTemp2BZ(:,:,kmax) = 0.0d0 !穴埋め xyr_F2BZ2 = & & xyr_avr_xyz( & & CpDry * xyz_DensBZ * ( xyz_VPTempBZ ** 2.0d0 ) & & ) info2 = "VI_init, xyr_CpDensVPTemp2BZ: " call CHK2_Val( xyr_CpDensVPTemp2BZ(1:nx,1:ny,1:nz) , xyr_F2BZ2(1:nx,1:ny,1:nz) , info2 ) xyr_F2BZ = xyr_CpDensVPTemp2BZ do k = kmin, kmax-1 do j = 1, ny do i = 1, nx xyr_DensVPTempBZ(i,j,k) = & & + ( & & + xyz_DensBZ(i,j,k+1) * xyz_VPTempBZ(i,j,k+1) & & + xyz_DensBZ(i,j,k) * xyz_VPTempBZ(i,j,k) & & ) * 5.0d-1 end do end do end do xyr_DensVPTempBZ(:,:,kmax) = 0.0d0 !穴埋め do k = kmin, kmax-1 do j = 1, ny do i = 1, nx xyr_CpVPTempBZ(i,j,k) = & & CpDry & & * ( & & + xyz_VPTempBZ(i,j,k+1) & & + xyz_VPTempBZ(i,j,k) & & ) * 5.0d-1 end do end do end do xyr_CpVPTempBZ(:,:,kmax) = 0.0d0 do k = 2, nz-1 A(k) = & & + 1.0d0 & & + ( beta ** 2.0d0 ) & & * xyz_F1BZ(nx,ny,k) * ( dt * dt ) & & * ( & & xyr_CpDensVPTemp2BZ(nx,ny,k) & & + xyr_CpDensVPTemp2BZ(nx,ny,k-1) & & ) & & / ( dz * dz ) end do A(1) = & & + 1.0d0 & & + ( beta ** 2.0d0 ) & & * xyz_F1BZ(nx,ny,1) * ( dt * dt ) & & * xyr_CpDensVPTemp2BZ(nx,ny,1) & & / ( dz * dz ) A(nz) = & & + 1.0d0 & & + ( beta ** 2.0d0 ) & & * xyz_F1BZ(nx,ny,nz) * ( dt * dt ) & & * xyr_CpDensVPTemp2BZ(nx,ny,nz-1) & & / ( dz * dz ) do k = 2, nz B(k) = & & - ( beta ** 2.0d0 ) & & * xyz_F1BZ(nx,ny,k-1) * ( dt * dt ) & & * xyr_CpDensVPTemp2BZ(nx,ny,k-1) & & / ( dz * dz ) end do do k = 1, nz-1 C(k) = & & - ( beta ** 2.0d0 ) & & * xyz_F1BZ(nx,ny,k+1) * (dt * dt ) & & * xyr_CpDensVPTemp2BZ(nx,ny,k) & & / ( dz * dz ) end do A2(1+1: nz-1) = & & (beta ** 2.0d0) & & * xyz_F1BZ2(nx,ny,1+1: nz-1) & & * (DTS ** 2.0d0) & & * ( & & xyr_F2BZ2(nx,ny,1+1: nz-1) & & / r_dz(1+1: nz-1) & & + xyr_F2BZ2(nx,ny,1 : nz-2) & & / r_dz(1: nz-2) & & ) & & / z_dz(1+1: nz-1) & & + 1.0d0 A2(1) = & & (beta ** 2.0d0) & & * xyz_F1BZ2(nx,ny,1) & & * xyr_F2BZ2(nx,ny,1) & & / r_dz(1) & & * (DTS ** 2.0d0) & & / z_dz(1) & & + 1.0d0 A2(nz) = & & (beta ** 2.0d0) & & * xyz_F1BZ2(nx,ny,nz) & & * xyr_F2BZ2(nx,ny,nz-1) & & / r_dz(nz-1) & & * (DTS ** 2.0d0) & & / z_dz(nz) & & + 1.0d0 B2(1+1:nz) = & & - (beta ** 2.0d0) & & * xyz_F1BZ2(nx,ny,1:nz-1) & & * xyr_F2BZ2(nx,ny,1:nz-1) & & * (DTS ** 2.0d0) & & / ( r_dz(1:nz-1) * z_dz(1:nz-1) ) C2(1: nz-1) = & & - ( beta ** 2.0d0 ) & & * xyz_F1BZ2(nx,ny,1+1:nz) & & * xyr_F2BZ2(nx,ny,1 :nz-1) & & * (DTS ** 2.0d0) & & / ( r_dz(1:nz-1) * z_dz(1+1:nz) ) write(*,*) "Dynamics_VI_init, A: ", & & minval( A(1:nz) - A2(1:nz) ), & & maxval( A(1:nz) - A2(1:nz) ), & & minval( (A(1:nz) - A2(1:nz)) / A2(1:nz)), & & maxval( (A(1:nz) - A2(1:nz)) / A2(1:nz)) write(*,*) "Dynamics_VI_init, B: ", & & minval( B(2:nz) - B2(2:nz) ), & & maxval( B(2:nz) - B2(2:nz) ), & & minval( (B(2:nz) - B2(2:nz)) / B2(2:nz)), & & maxval( (B(2:nz) - B2(2:nz)) / B2(2:nz)) write(*,*) "Dynamics_VI_init, C: ", & & minval( C(1:nz-1) - C2(1:nz-1) ), & & maxval( C(1:nz-1) - C2(1:nz-1) ), & & minval( (C(1:nz-1) - C2(1:nz-1)) / C2(1:nz-1)), & & maxval( (C(1:nz-1) - C2(1:nz-1)) / C2(1:nz-1)) ! write(*,*) A(1), A2(1), C(1), C2(1) ! do k = 1, nz-1 ! write(*,*) A(k), A2(k), B(k), B2(k), C(k), C2(k) ! end do ! write(*,*) A(nz), A2(nz), B(nz), B2(nz) !---------------------------------------------------------------- ! 係数行列を LU 分解 ! !配列の大きさを保管 N = nz !係数行列/改行列の次数, 整合寸法 M = nx * ny !方程式の組数 NUD = 1 !係数行列の上三角部分の帯幅 NLD = 1 !係数行列の下三角部分の帯幅 NAL = NLD !LU 分解の結果 L の整合寸法 NA = NUD + NLD + 1 INFO = 0 ! 配列の割り当て ! allocate( AL1(N), IP(N) ) ! 解行列の計算. LAPACK を使用. ! call DGTTRF(N, C, A, B, AL1, IP, INFO) ! 解のコンディションをチェック. ! if (INFO /= 0) then call MessageNotify("Error", "lapack_linear", "INFO is not 0") stop end if end subroutine Dynamics_VI_init !!!--------------------------------------------------------------------!!! subroutine Exner_HEVI( & & pyz_VelXAs, xqz_VelYAs, xyr_VelZNs, & & xyz_VelDivNs, xyz_ExnerNs, & & xyr_DVelZDtNl, xyz_DExnerDtNl, & & xyz_ExnerAs & & ) ! !陰解法を用いたエクスナー関数の計算. ! !暗黙の型宣言禁止 implicit none !入出力変数 real(DP), intent(in) :: pyz_VelXAs(imin:imax,jmin:jmax,kmin:kmax) !速度 u [τ+Δτ] real(DP), intent(in) :: xqz_VelYAs(imin:imax,jmin:jmax,kmin:kmax) !速度 v [τ+Δτ] real(DP), intent(in) :: xyr_VelZNs(imin:imax,jmin:jmax,kmin:kmax) !速度 w [τ] real(DP), intent(in) :: xyz_VelDivNs(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax) !Z 方向の外力項[t] real(DP), intent(in) :: xyz_DExnerDtNl(imin:imax,jmin:jmax,kmin:kmax) !Z 方向の外力項[t] real(DP), intent(in) :: xyz_ExnerNs(imin:imax,jmin:jmax,kmin:kmax) !無次元圧力 real(DP), intent(out) :: xyz_ExnerAs(imin:imax,jmin:jmax,kmin:kmax) !無次元圧力[τ+Δτ] !作業変数定義 real(DP) :: D1(1:nx,1:ny,1:nz) real(DP) :: D(nx*ny,nz) real(DP) :: E(1:nx,1:ny,0:nz) ! real(DP) :: E(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: F(1:nx,1:ny,1:nz) real(DP) :: F0(1:nx,1:ny,kmin:kmax-1) ! real(DP) :: xyz_DivVelNs(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: dt ! 短い時間格子間隔 integer :: i, j, k real(DP) :: X(M, N) !定数/解行列 real(DP) :: TX(N, M) !解行列を転置したもの integer :: NRHS integer :: INFO integer :: LDB character(1),parameter :: TRANS = 'N' real(DP) :: dts real(DP) :: nD1(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: nD2(1:nx,1:ny,1:nz) ! real(DP) :: nD(nx*ny,nz) real(DP) :: nE(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: nF(imin:imax,jmin:jmax,kmin:kmax) character(15) :: info2 real(DP) :: tmpA(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: tmpB(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: tmpC(imin:imax,jmin:jmax,kmin:kmax) ! Initialize ! dt = DelTimeShort dts = DelTimeShort !--------------------------------------------------------------- !行列計算のための係数 ! 添字の範囲は, 1:nx, 1:ny, 0:nz ! D(:,:,1) を求める時に D(:,:,0) の値が必要になるため. ! tmpA = xyr_dz_xyz( xyz_ExnerNs ) tmpB = xyr_dz_xyz( xyz_VelDivNs ) tmpC = xyr_avr_xyz( CpDry * xyz_VPTempBZ ) do k = kmin, kmax-1 do j = jmin, jmax-1 do i = imin, imax-1 tmpC(i,j,k) = & & ( & & + xyz_VelDivNs(i,j,k+1) & & - xyz_VelDivNs(i,j,k) & & ) / dz end do end do end do tmpC(imax,:,:) = 1.0d10 tmpC(:,jmax,:) = 1.0d10 tmpC(:,:,kmax) = 1.0d10 ! i = 25 ! j = 23 ! write(*,*) "***", kmin ! do k = kmin, kmax-1 ! write(*,*) k, xyz_VelDivNs(i,j,k+1), xyz_VelDivNs(i,j,k), tmpC(i,j,k), tmpB(i,j,k), tmpC(i,j,k) - tmpB(i,j,k) ! end do do k = 0, nz do j = 1, ny do i = 1, nx E(i,j,k) = & & - ( 1.0d0 - beta ) & & * ( & & + xyz_ExnerNs(i,j,k+1) & ! & - xyz_ExnerNs(i,j,k) & ! xyz => xyr & ) / dz & & + ( & & + Alpha & & * ( & & + xyz_VelDivNs(i,j,k+1) & ! & - xyz_VelDivNs(i,j,k) & ! xyz => xyr & ) / dz & & + xyr_DVelZDtNl(i,j,k) & & ) / xyr_CpVPTempBZ(i,j,k) end do end do end do nE = & & - ( 1.0d0 - beta ) * xyr_dz_xyz( xyz_ExnerNs ) & & + ( Alpha * xyr_dz_xyz( xyz_VelDivNs ) + xyr_DVelZDtNl ) & & / xyr_avr_xyz( CpDry * xyz_VPTempBZ ) info2 = "VI, E: " call CHK2_Val( E(1:nx,1:ny,1:nz), nE(1:nx,1:ny,1:nz), info2 ) ! 被微分関数 ! 配列 F0 の添字の範囲は, 1:nx, 1:ny, kmin:kmax ! 配列 F を求める際に F0 を z 方向に微分するため. ! do k = kmin, kmax-1 do j = 1, ny do i = 1, nx F0(i,j,k) = & & + xyr_DensVPTempBZ(i,j,k) & & * ( & & + xyr_VelZNs(i,j,k) & & - xyr_CpVPTempBZ(i,j,k) & & * (1.0d0 - beta) & & * ( & & xyz_ExnerNs(i,j,k+1) & ! & - xyz_ExnerNs(i,j,k) & ! xyz => xyr & ) / dz * dt & & + Alpha & & * ( & & xyz_VelDivNs(i,j,k+1) & & - xyz_VelDivNs(i,j,k) & & ) / dz * dt & & + xyr_DVelZDtNl(i,j,k) * dt & & ) end do end do end do !行列計算のための係数 ! 添字の範囲は, 1:nx, 1:ny, 1:nz ! do k = 1, nz do j = 1, ny do i = 1, nx F(i,j,k) = & & - beta * xyz_F1BZ(i,j,k) * dt & & * ( & & F0(i,j,k) & & - F0(i,j,k-1) & & ) / dz & & + xyz_DExnerDtNl(i,j,k) * dt end do end do end do !行列計算のための係数 ! 添字の範囲は, 1:nx, 1:ny, 1:nz ! do k = 1, nz do j = 1, ny do i = 1, nx D1(i,j,k) = & & + xyz_ExnerNs(i,j,k) & & - (1.0d0 - beta) & & * xyz_F1BZ(i,j,k) * dt & & * ( & & xyr_DensVPTempBZ(i,j,k) & & * xyr_VelZNs(i,j,k) & & - xyr_DensVPTempBZ(i,j,k-1) & & * xyr_VelZNs(i,j,k-1) & & ) / dz & & - (xyz_VelSW(i,j,k) ** 2.0d0) * dt & & / (CpDry * xyz_VPTempBZ(i,j,k)) & & * ( & & + ( & & pyz_VelXAs(i,j,k) & & - pyz_VelXAs(i-1,j,k) & & ) / dx & & + ( & & xqz_VelYAs(i,j,k) & & - xqz_VelYAs(i,j-1,k) & & ) / dy & & ) & & + F(i,j,k) end do end do end do ! 行列計算のための係数 ! do j = 1, ny do i = 1, nx D1(i,j,1) = & & + D1(i,j,1) & & - beta * xyz_F1BZ(i,j,1) * (dt ** 2.0d0) & & * xyr_CpDensVPTemp2BZ(i,j,0) & & * E(i,j,0) & & / dz D1(i,j,nz) = & & + D1(i,j,nz) & & + beta * xyz_F1BZ(i,j,nz) * (dt ** 2.0d0) & & * xyr_CpDensVPTemp2BZ(i,j,nz) & & * E(i,j,nz) & & / dz end do end do nF = - beta * xyz_F3BZ * DTS & & * xyz_dz_xyr( & & xyr_avr_xyz( xyz_DensBZ * xyz_VPTempBZ) & & * ( & & xyr_VelZNs & & - xyr_avr_xyz(CpDry * xyz_VPTempBZ) * DTS & & * (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerNs ) & & + Alpha * xyr_dz_xyz( xyz_VelDivNs ) * DTS & & + xyr_DVelZDtNl * DTS & & ) & & ) & & + ( xyz_DExnerDtNl ) * DTS nD1 = xyz_ExnerNs & & - (1.0d0 - beta) & & * xyz_F3BZ * DTS & & * xyz_dz_xyr( & & xyr_avr_xyz(xyz_DensBZ * xyz_VPTempBZ) * xyr_VelZNs & & ) & & - (xyz_VelSoundBZ ** 2.0d0) * DTS & & / (CpDry * xyz_VPTempBZ) * xyz_dx_pyz( pyz_VelXAs ) & & - (xyz_VelSoundBZ ** 2.0d0) * DTS & & / (CpDry * xyz_VPTempBZ) * xyz_dy_xqz( xqz_VelYAs ) & & + nF nD1(1:nx,1:ny,1) = nD1(1:nx,1:ny,1) & & - beta * xyz_F3BZ(1:nx,1:ny,1) * (DTS ** 2.0d0) & & * xyr_F2BZ(1:nx,1:ny, 1-1) * nE(1:nx,1:ny,1-1) & & / z_dz(1) nD1(1:nx,1:ny,nz) = nD1(1:nx,1:ny,nz) & & + beta * xyz_F3BZ(1:nx,1:ny,nz) * (DTS ** 2.0d0) & & * xyr_F2BZ(1:nx,1:ny,nz) * nE(1:nx,1:ny,nz) & & / z_dz(nz) nD2 = nD1(1:nx,1:ny,1:nz) info2 = "VI, F: " call CHK2_Val( F(1:nx,1:ny,1:nz), nF(1:nx,1:ny,1:nz), info2 ) info2 = "VI, D1: " call CHK2_Val( D1(1:nx,1:ny,1:nz), nD2(1:nx,1:ny,1:nz), info2 ) !----------------------------------------------------------- !連立一次方程式の解を求める !変数の初期化 ! NRHS = M INFO = 0 LDB = N ! LAPACK の仕様に合わせて変形 ! do k = 1, nz do j = 1, ny do i = 1, nx D(i + nx * (j - 1), k) = D1(i,j,k) end do end do end do TX = transpose( D ) !解行列の計算. LAPACK を使用. ! call DGTTRS(TRANS, N, NRHS, C, A, B, AL1, IP, TX, LDB, INFO) !解のコンディションをチェック. ! if (INFO /= 0) then call MessageNotify("Error", "lapack_linear", "INFO is not 0") stop end if !戻り値を出力 ! X = transpose( TX ) do k = 1, nz do j = 1, ny do i = 1, nx xyz_ExnerAs(i,j,k) = X(i + nx * (j - 1 ), k) end do end do end do xyz_ExnerAs(imin:0,:,:) = 0.0d0 xyz_ExnerAs(:,jmin:0,:) = 0.0d0 xyz_ExnerAs(:,:,kmin:0) = 0.0d0 xyz_ExnerAs(nx+1:imax,:,:) = 0.0d0 xyz_ExnerAs(:,ny+1:jmax,:) = 0.0d0 xyz_ExnerAs(:,:,nz+1:kmax) = 0.0d0 ! のり代に値を入れる ! call SetMargin_xyz( xyz_ExnerAs ) ! (inout) end subroutine Exner_HEVI subroutine AdvC4_nDiff_xyz( & & xyz_VarBl, xyz_VarNl, & !(IN) & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN) & xyz_Adv, xyz_nDiff ) !(OUT) implicit none real(DP), intent(in) :: xyz_VarBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_VarNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: fct1, fct2 integer :: i, j, k ! 微分に用いる係数を予め計算 ! fct1 = 9.0d0 / 8.0d0 fct2 = 1.0d0 / 24.0d0 ! 移流項の計算. 移流項: 4 次精度中心差分 ! do k = kmin + 2, kmax - 2 do j = jmin + 2, jmax - 2 do i = imin + 2, imax - 2 xyz_Adv(i,j,k) = & & - ( & & pyz_VelXNl(i,j,k) & & * ( & & fct1 * ( xyz_VarNl(i+1,j,k) - xyz_VarNl(i,j,k) ) & & - fct2 * ( xyz_VarNl(i+2,j,k) - xyz_VarNl(i-1,j,k) ) & & ) & & + pyz_VelXNl(i-1,j,k) & & * ( & & fct1 * ( xyz_VarNl(i,j,k) - xyz_VarNl(i-1,j,k) ) & & - fct2 * ( xyz_VarNl(i+1,j,k) - xyz_VarNl(i-2,j,k) ) & & ) & & ) * 5.0d-1 / dx & & - ( & & xqz_VelYNl(i,j,k) & & * ( & & fct1 * ( xyz_VarNl(i,j+1,k) - xyz_VarNl(i,j,k) ) & & - fct2 * ( xyz_VarNl(i,j+2,k) - xyz_VarNl(i,j-1,k) ) & & ) & & + xqz_VelYNl(i,j-1,k) & & * ( & & fct1 * ( xyz_VarNl(i,j,k) - xyz_VarNl(i,j-1,k) ) & & - fct2 * ( xyz_VarNl(i,j+1,k) - xyz_VarNl(i,j-2,k) ) & & ) & & ) * 5.0d-1 / dy & & - ( & & xyr_VelZNl(i,j,k) & & * ( & & fct1 * ( xyz_VarNl(i,j,k+1) - xyz_VarNl(i,j,k) ) & & - fct2 * ( xyz_VarNl(i,j,k+2) - xyz_VarNl(i,j,k-1) ) & & ) & & + xyr_VelZNl(i,j,k-1) & & * ( & & fct1 * ( xyz_VarNl(i,j,k) - xyz_VarNl(i,j,k-1) ) & & - fct2 * ( xyz_VarNl(i,j,k+1) - xyz_VarNl(i,j,k-2) ) & & ) & & ) * 5.0d-1 / dz end do end do end do xyz_Adv(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0 xyz_Adv(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0 xyz_Adv(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0 xyz_Adv(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0 xyz_Adv(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0 xyz_Adv(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0 ! 4 次の数値拡散: 2 次精度中心差分 ! do k = kmin + 2, kmax - 2 do j = jmin + 2, jmax - 2 do i = imin + 2, imax - 2 xyz_nDiff(i,j,k) = & & - ( & & + xyz_VarBl(i+2,j,k) & & + xyz_VarBl(i-2,j,k) & & - xyz_VarBl(i+1,j,k) * 4.0d0 & & - xyz_VarBl(i-1,j,k) * 4.0d0 & & + xyz_VarBl(i ,j,k) * 6.0d0 & & ) * NuHh / ( dx ** 4.0d0 ) & & - ( & & xyz_VarBl(i,j+2,k) & & + xyz_VarBl(i,j-2,k) & & - xyz_VarBl(i,j+1,k) * 4.0d0 & & - xyz_VarBl(i,j-1,k) * 4.0d0 & & + xyz_VarBl(i,j ,k) * 6.0d0 & & ) * NuHh / ( dy ** 4.0d0 ) & & - ( & & xyz_VarBl(i,j,k+2) & & + xyz_VarBl(i,j,k-2) & & - xyz_VarBl(i,j,k+1) * 4.0d0 & & - xyz_VarBl(i,j,k-1) * 4.0d0 & & + xyz_VarBl(i,j,k ) * 6.0d0 & & ) * NuVh / ( dz ** 4.0d0 ) end do end do end do xyz_nDiff(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0 xyz_nDiff(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0 xyz_nDiff(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0 xyz_nDiff(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0 xyz_nDiff(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0 xyz_nDiff(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0 end subroutine AdvC4_nDiff_xyz subroutine AdvC4_nDiff_xyzf( & & xyzf_VarBl, xyzf_VarNl, & !(IN) & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN) & xyzf_Adv, xyzf_nDiff ) !(OUT) implicit none real(DP), intent(in) :: xyzf_VarBl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) real(DP), intent(in) :: xyzf_VarNl(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) real(DP), intent(in) :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyzf_Adv(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) real(DP), intent(out) :: xyzf_nDiff(imin:imax,jmin:jmax,kmin:kmax,1:ncmax) real(DP) :: fct1, fct2 integer :: s, i, j, k ! 微分に用いる係数を予め計算 ! fct1 = 9.0d0 / 8.0d0 fct2 = 1.0d0 / 24.0d0 ! 移流項の計算. 移流項: 4 次精度中心差分 ! do s = 1, ncmax do k = kmin + 2, kmax - 2 do j = jmin + 2, jmax - 2 do i = imin + 2, imax - 2 xyzf_Adv(i,j,k,s) = & & - ( & & pyz_VelXNl(i,j,k) & & * ( & & fct1 * ( xyzf_VarNl(i+1,j,k,s) - xyzf_VarNl(i,j,k,s) ) & & - fct2 * ( xyzf_VarNl(i+2,j,k,s) - xyzf_VarNl(i-1,j,k,s) ) & & ) & & + pyz_VelXNl(i-1,j,k) & & * ( & & fct1 * ( xyzf_VarNl(i,j,k,s) - xyzf_VarNl(i-1,j,k,s) ) & & - fct2 * ( xyzf_VarNl(i+1,j,k,s) - xyzf_VarNl(i-2,j,k,s) ) & & ) & & ) * 5.0d-1 / dx & & - ( & & xqz_VelYNl(i,j,k) & & * ( & & fct1 * ( xyzf_VarNl(i,j+1,k,s) - xyzf_VarNl(i,j,k,s) ) & & - fct2 * ( xyzf_VarNl(i,j+2,k,s) - xyzf_VarNl(i,j-1,k,s) ) & & ) & & + xqz_VelYNl(i,j-1,k) & & * ( & & fct1 * ( xyzf_VarNl(i,j,k,s) - xyzf_VarNl(i,j-1,k,s) ) & & - fct2 * ( xyzf_VarNl(i,j+1,k,s) - xyzf_VarNl(i,j-2,k,s) ) & & ) & & ) * 5.0d-1 / dy & & - ( & & xyr_VelZNl(i,j,k) & & * ( & & fct1 * ( xyzf_VarNl(i,j,k+1,s) - xyzf_VarNl(i,j,k,s) ) & & - fct2 * ( xyzf_VarNl(i,j,k+2,s) - xyzf_VarNl(i,j,k-1,s) ) & & ) & & + xyr_VelZNl(i,j,k-1) & & * ( & & fct1 * ( xyzf_VarNl(i,j,k,s) - xyzf_VarNl(i,j,k-1,s) ) & & - fct2 * ( xyzf_VarNl(i,j,k+1,s) - xyzf_VarNl(i,j,k-2,s) ) & & ) & & ) * 5.0d-1 / dz end do end do end do end do xyzf_Adv(imin:imin+1,jmin:jmax,kmin:kmax,1:ncmax) = 0.0d0 xyzf_Adv(imax-1:imax,jmin:jmax,kmin:kmax,1:ncmax) = 0.0d0 xyzf_Adv(imin:imax,jmin:jmin+1,kmin:kmax,1:ncmax) = 0.0d0 xyzf_Adv(imin:imax,jmax-1:jmax,kmin:kmax,1:ncmax) = 0.0d0 xyzf_Adv(imin:imax,jmin:jmax,kmin:kmin+1,1:ncmax) = 0.0d0 xyzf_Adv(imin:imax,jmin:jmax,kmax-1:kmax,1:ncmax) = 0.0d0 ! 数値拡散: 2 次精度中心差分 ! do s = 1, ncmax do k = kmin + 2, kmax - 2 do j = jmin + 2, jmax - 2 do i = imin + 2, imax - 2 xyzf_nDiff(i,j,k,s) = & & - ( & & xyzf_VarBl(i+2,j,k,s) & & - xyzf_VarBl(i+1,j,k,s) * 4.0d0 & & + xyzf_VarBl(i ,j,k,s) * 6.0d0 & & - xyzf_VarBl(i-1,j,k,s) * 4.0d0 & & + xyzf_VarBl(i-2,j,k,s) & & ) * NuHh / ( dx ** 4.0d0 ) & & - ( & & xyzf_VarBl(i,j+2,k,s) & & - xyzf_VarBl(i,j+1,k,s) * 4.0d0 & & + xyzf_VarBl(i,j ,k,s) * 6.0d0 & & - xyzf_VarBl(i,j-1,k,s) * 4.0d0 & & + xyzf_VarBl(i,j-2,k,s) & & ) * NuHh / ( dy ** 4.0d0 ) & & - ( & & xyzf_VarBl(i,j,k+2,s) & & - xyzf_VarBl(i,j,k+1,s) * 4.0d0 & & + xyzf_VarBl(i,j,k ,s) * 6.0d0 & & - xyzf_VarBl(i,j,k-1,s) * 4.0d0 & & + xyzf_VarBl(i,j,k-2,s) & & ) * NuVh / ( dz ** 4.0d0 ) end do end do end do end do xyzf_nDiff(imin:imin+1,jmin:jmax,kmin:kmax,1:ncmax) = 0.0d0 xyzf_nDiff(imax-1:imax,jmin:jmax,kmin:kmax,1:ncmax) = 0.0d0 xyzf_nDiff(imin:imax,jmin:jmin+1,kmin:kmax,1:ncmax) = 0.0d0 xyzf_nDiff(imin:imax,jmax-1:jmax,kmin:kmax,1:ncmax) = 0.0d0 xyzf_nDiff(imin:imax,jmin:jmax,kmin:kmin+1,1:ncmax) = 0.0d0 xyzf_nDiff(imin:imax,jmin:jmax,kmax-1:kmax,1:ncmax) = 0.0d0 end subroutine AdvC4_nDiff_xyzf subroutine AdvC4_nDiff_pyz_xqz_xyr( & & pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, & !(IN) & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN) & pyz_Adv, xqz_Adv, xyr_Adv, & !(OUT) & pyz_nDiff, xqz_nDiff, xyr_nDiff ) !(OUT) 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) :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: pyz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xqz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyr_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: pyz_nDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xqz_nDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyr_nDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: fct1, fct2 integer :: i, j, k ! 微分に用いる係数を予め計算 ! fct1 = 9.0d0 / 8.0d0 fct2 = 1.0d0 / 24.0d0 ! 移流項の計算. 移流項: 4 次精度中心差分 ! do k = kmin + 2, kmax - 2 do j = jmin + 2, jmax - 2 do i = imin + 2, imax - 2 pyz_Adv(i,j,k) = & & - pyz_VelXNl(i,j,k) & & * ( & & fct1 * ( pyz_VelXNl(i+1,j,k) - pyz_VelXNl(i-1,j,k) ) & & - fct2 * ( pyz_VelXNl(i+2,j,k) + pyz_VelXNl(i+1,j,k) & & - pyz_VelXNl(i-1,j,k) - pyz_VelXNl(i-2,j,k) ) & & ) * 5.0d-1 / dx & & - ( & & + ( xqz_VelYNl(i+1,j,k) + xqz_VelYNl(i,j,k) ) & & * ( & & fct1 * ( pyz_VelXNl(i,j+1,k) - pyz_VelXNl(i,j,k) ) & & - fct2 * ( pyz_VelXNl(i,j+2,k) - pyz_VelXNl(i,j-1,k) ) & & ) & & + ( xqz_VelYNl(i+1,j-1,k) + xqz_VelYNl(i,j-1,k) ) & & * ( & & fct1 * ( pyz_VelXNl(i,j,k) - pyz_VelXNl(i,j-1,k) ) & & - fct2 * ( pyz_VelXNl(i,j+1,k) - pyz_VelXNl(i,j-2,k) ) & & ) & & ) * 2.5d-1 / dy & & - ( & & + ( xyr_VelZNl(i+1,j,k) + xyr_VelZNl(i,j,k) ) & & * ( & & fct1 * ( pyz_VelXNl(i,j,k+1) - pyz_VelXNl(i,j,k) ) & & - fct2 * ( pyz_VelXNl(i,j,k+2) - pyz_VelXNl(i,j,k-1) ) & & ) & & + ( xyr_VelZNl(i+1,j,k-1) + xyr_VelZNl(i,j,k-1) ) & & * ( & & fct1 * ( pyz_VelXNl(i,j,k) - pyz_VelXNl(i,j,k-1) ) & & - fct2 * ( pyz_VelXNl(i,j,k+1) - pyz_VelXNl(i,j,k-2) ) & & ) & & ) * 2.5d-1 / dz end do end do end do pyz_Adv(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0 pyz_Adv(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0 pyz_Adv(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0 pyz_Adv(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0 pyz_Adv(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0 pyz_Adv(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0 do k = kmin + 2, kmax - 2 do j = jmin + 2, jmax - 2 do i = imin + 2, imax - 2 xqz_Adv(i,j,k) = & & - ( & & + ( pyz_VelXNl(i,j+1,k) + pyz_VelXNl(i,j,k) ) & & * ( & & fct1 * ( xqz_VelYNl(i+1,j,k) - xqz_VelYNl(i,j,k) ) & & - fct2 * ( xqz_VelYNl(i+2,j,k) - xqz_VelYNl(i-1,j,k) ) & & ) & & + ( pyz_VelXNl(i-1,j+1,k) + pyz_VelXNl(i-1,j,k) ) & & * ( & & fct1 * ( xqz_VelYNl(i,j,k) - xqz_VelYNl(i-1,j,k) ) & & - fct2 * ( xqz_VelYNl(i+1,j,k) - xqz_VelYNl(i-2,j,k) ) & & ) & & ) * 2.5d-1 / dx & & - xqz_VelYNl(i,j,k) & & * ( & & fct1 * ( xqz_VelYNl(i,j+1,k) - xqz_VelYNl(i,j-1,k) ) & & - fct2 * ( xqz_VelYNl(i,j+2,k) + xqz_VelYNl(i,j+1,k) & & - xqz_VelYNl(i,j-1,k) - xqz_VelYNl(i,j-2,k) ) & & ) * 5.0d-1 / dy & & - ( & & + ( xyr_VelZNl(i,j+1,k) + xyr_VelZNl(i,j,k) ) & & * ( & & fct1 * ( xqz_VelYNl(i,j,k+1) - xqz_VelYNl(i,j,k) ) & & - fct2 * ( xqz_VelYNl(i,j,k+2) - xqz_VelYNl(i,j,k-1) ) & & ) & & + ( xyr_VelZNl(i,j+1,k-1) + xyr_VelZNl(i,j,k-1) ) & & * ( & & fct1 * ( xqz_VelYNl(i,j,k) - xqz_VelYNl(i,j,k-1) ) & & - fct2 * ( xqz_VelYNl(i,j,k+1) - xqz_VelYNl(i,j,k-2) ) & & ) & & ) * 2.5d-1 / dz end do end do end do xqz_Adv(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0 xqz_Adv(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0 xqz_Adv(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0 xqz_Adv(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0 xqz_Adv(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0 xqz_Adv(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0 do k = kmin + 2, kmax - 2 do j = jmin + 2, jmax - 2 do i = imin + 2, imax - 2 xyr_Adv(i,j,k) = & & - ( & & + ( pyz_VelXNl(i,j,k+1) + pyz_VelXNl(i,j,k) ) & & * ( & & fct1 * ( xyr_VelZNl(i+1,j,k) - xyr_VelZNl(i,j,k) ) & & - fct2 * ( xyr_VelZNl(i+2,j,k) - xyr_VelZNl(i-1,j,k) ) & & ) & & + ( pyz_VelXNl(i-1,j,k+1) + pyz_VelXNl(i-1,j,k) ) & & * ( & & fct1 * ( xyr_VelZNl(i,j,k) - xyr_VelZNl(i-1,j,k) ) & & - fct2 * ( xyr_VelZNl(i+1,j,k) - xyr_VelZNl(i-2,j,k) ) & & ) & & ) * 2.5d-1 / dx & & - ( & & + ( xqz_VelYNl(i,j,k+1) + xqz_VelYNl(i,j,k) ) & & * ( & & fct1 * ( xyr_VelZNl(i,j+1,k) - xyr_VelZNl(i,j,k) ) & & - fct2 * ( xyr_VelZNl(i,j+2,k) - xyr_VelZNl(i,j-1,k) ) & & ) & & + ( xqz_VelYNl(i,j-1,k+1) + xqz_VelYNl(i,j-1,k) ) & & * ( & & fct1 * ( xyr_VelZNl(i,j,k) - xyr_VelZNl(i,j-1,k) ) & & - fct2 * ( xyr_VelZNl(i,j+1,k) - xyr_VelZNl(i,j-2,k) ) & & ) & & ) * 2.5d-1 / dy & & - xyr_VelZNl(i,j,k) & & * ( & & fct1 * ( xyr_VelZNl(i,j,k+1) - xyr_VelZNl(i,j,k-1) ) & & - fct2 * ( xyr_VelZNl(i,j,k+2) + xyr_VelZNl(i,j,k+1) & & - xyr_VelZNl(i,j,k-1) - xyr_VelZNl(i,j,k-2) ) & & ) * 5.0d-1 / dz end do end do end do xyr_Adv(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0 xyr_Adv(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0 xyr_Adv(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0 xyr_Adv(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0 xyr_Adv(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0 xyr_Adv(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0 do k = kmin + 2, kmax - 2 do j = jmin + 2, jmax - 2 do i = imin + 2, imax - 2 pyz_nDiff(i,j,k) = & & - ( & & + pyz_VelXBl(i+2,j,k) & & + pyz_VelXBl(i-2,j,k) & & - pyz_VelXBl(i+1,j,k) * 4.0d0 & & - pyz_VelXBl(i-1,j,k) * 4.0d0 & & + pyz_VelXBl(i ,j,k) * 6.0d0 & & ) * NuHm / ( dx ** 4.0d0 ) & & - ( & & + pyz_VelXBl(i,j+2,k) & & + pyz_VelXBl(i,j-2,k) & & - pyz_VelXBl(i,j+1,k) * 4.0d0 & & - pyz_VelXBl(i,j-1,k) * 4.0d0 & & + pyz_VelXBl(i,j ,k) * 6.0d0 & & ) * NuHm / ( dy ** 4.0d0 ) & & - ( & & + pyz_VelXBl(i,j,k+2) & & + pyz_VelXBl(i,j,k-2) & & - pyz_VelXBl(i,j,k+1) * 4.0d0 & & - pyz_VelXBl(i,j,k-1) * 4.0d0 & & + pyz_VelXBl(i,j,k ) * 6.0d0 & & ) * NuVm / ( dz ** 4.0d0 ) end do end do end do pyz_nDiff(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0 pyz_nDiff(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0 pyz_nDiff(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0 pyz_nDiff(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0 pyz_nDiff(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0 pyz_nDiff(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0 do k = kmin + 2, kmax - 2 do j = jmin + 2, jmax - 2 do i = imin + 2, imax - 2 xqz_nDiff(i,j,k) = & & - ( & & + xqz_VelYBl(i+2,j,k) & & + xqz_VelYBl(i-2,j,k) & & - xqz_VelYBl(i+1,j,k) * 4.0d0 & & - xqz_VelYBl(i-1,j,k) * 4.0d0 & & + xqz_VelYBl(i ,j,k) * 6.0d0 & & ) * NuHm / ( dx ** 4.0d0 ) & & - ( & & + xqz_VelYBl(i,j+2,k) & & + xqz_VelYBl(i,j-2,k) & & - xqz_VelYBl(i,j+1,k) * 4.0d0 & & - xqz_VelYBl(i,j-1,k) * 4.0d0 & & + xqz_VelYBl(i,j ,k) * 6.0d0 & & ) * NuHm / ( dy ** 4.0d0 ) & & - ( & & + xqz_VelYBl(i,j,k+2) & & + xqz_VelYBl(i,j,k-2) & & - xqz_VelYBl(i,j,k+1) * 4.0d0 & & - xqz_VelYBl(i,j,k-1) * 4.0d0 & & + xqz_VelYBl(i,j,k ) * 6.0d0 & & ) * NuVm / ( dz ** 4.0d0 ) end do end do end do xqz_nDiff(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0 xqz_nDiff(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0 xqz_nDiff(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0 xqz_nDiff(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0 xqz_nDiff(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0 xqz_nDiff(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0 do k = kmin + 2, kmax - 2 do j = jmin + 2, jmax - 2 do i = imin + 2, imax - 2 xyr_nDiff(i,j,k) = & & - ( & & + xyr_VelZBl(i+2,j,k) & & + xyr_VelZBl(i-2,j,k) & & - xyr_VelZBl(i+1,j,k) * 4.0d0 & & - xyr_VelZBl(i-1,j,k) * 4.0d0 & & + xyr_VelZBl(i ,j,k) * 6.0d0 & & ) * NuHm / ( dx ** 4.0d0 ) & & - ( & & + xyr_VelZBl(i,j+2,k) & & + xyr_VelZBl(i,j-2,k) & & - xyr_VelZBl(i,j+1,k) * 4.0d0 & & - xyr_VelZBl(i,j-1,k) * 4.0d0 & & + xyr_VelZBl(i,j ,k) * 6.0d0 & & ) * NuHm / ( dy ** 4.0d0 ) & & - ( & & + xyr_VelZBl(i,j,k+2) & & + xyr_VelZBl(i,j,k-2) & & - xyr_VelZBl(i,j,k+1) * 4.0d0 & & - xyr_VelZBl(i,j,k-1) * 4.0d0 & & + xyr_VelZBl(i,j,k ) * 6.0d0 & & ) * NuVm / ( dz ** 4.0d0 ) end do end do end do xyr_nDiff(imin:imin+1,jmin:jmax,kmin:kmax) = 0.0d0 xyr_nDiff(imax-1:imax,jmin:jmax,kmin:kmax) = 0.0d0 xyr_nDiff(imin:imax,jmin:jmin+1,kmin:kmax) = 0.0d0 xyr_nDiff(imin:imax,jmax-1:jmax,kmin:kmax) = 0.0d0 xyr_nDiff(imin:imax,jmin:jmax,kmin:kmin+1) = 0.0d0 xyr_nDiff(imin:imax,jmin:jmax,kmax-1:kmax) = 0.0d0 end subroutine AdvC4_nDiff_pyz_xqz_xyr subroutine BuoyancyLong_xyr( & & xyz_PTempNl, xyzf_QMixNl, & !(IN) & xyr_BuoyT, xyr_BuoyM, xyr_BuoyD ) !(OUT) use composition, only: GasNum, &! & IdxG, &! & MolWtWet ! 湿潤成分の分子量 use constants,only: MolWtDry, &! 乾燥成分の分子量 & Grav ! 重力加速度 implicit none real(DP), intent(in) :: xyz_PTempNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyzf_QMixNl(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP), intent(out) :: xyr_BuoyT(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyr_BuoyM(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyr_BuoyD(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyzf_QMixPerMolWt(imin:imax,jmin:jmax,kmin:kmax, 1:GasNum) real(DP) :: tmp1(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: tmp2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: tmp3(imin:imax,jmin:jmax,kmin:kmax) integer :: i, j, k, f, n real(DP) :: xyr_BuoyT2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_BuoyM2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_BuoyD2(imin:imax,jmin:jmax,kmin:kmax) character(15) :: info2 do f = 1, GasNum n = IdxG(f) xyzf_QMixPerMolWt(:,:,:,f) = xyzf_QMixNl(:,:,:,n) / MolWtWet(n) end do ! 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_PTempBZ(i,j,k+1) & & + xyz_PTempNl(i,j,k) / xyz_PTempBZ(i,j,k) & & ) * 5.0d-1 end do end do end do xyr_BuoyT(:,:,kmax) = 0.0d0 ! Buoyancy due to molecular weight ! tmp1 = sum(xyzf_QMixPerMolWt, 4) tmp2 = 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 & & * ( tmp1(i,j,k+1) + tmp1(i,j,k) ) * 5.0d-1 & & / ( 1.0d0 / MolWtDry + xyr_QMixBZPerMolWt(i,j,k) ) & & - Grav & & * ( tmp2(i,j,k+1) + tmp2(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 ! tmp3 = 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 & & * ( tmp3(i,j,k+1) + tmp3(i,j,k) ) * 5.0d-1 & & / ( 1.0d0 + xyr_QMixBZ(i,j,k) ) end do end do end do xyr_BuoyD(:,:,kmax) = 0.0d0 xyr_BuoyT2 = Grav * xyr_avr_xyz( xyz_PTempNl / xyz_PTempBZ) xyr_BuoyM2 = & & + Grav * xyr_avr_xyz( sum(xyzf_QMixPerMolWt, 4) ) & & / ( 1.0d0 / MolWtDry + xyr_QMixBZPerMolWt ) & & - Grav * xyr_avr_xyz( sum(xyzf_QMixNl(:,:,:,1:GasNum), 4) ) & & / ( 1.0d0 + xyr_QmixBZ ) xyr_BuoyD2 = & & - Grav * xyr_avr_xyz( sum(xyzf_QMixNl(:,:,:,GasNum+1:ncmax), 4) ) & & / ( 1.0d0 + xyr_QMixBZ ) info2 = "VelZ, BuoyT: " call CHK_Val( xyr_BuoyT(imin:imax,jmin:jmax,kmin:kmax), xyr_BuoyT2(imin:imax,jmin:jmax,kmin:kmax), info2 ) info2 = "VelZ, BuoyM: " call CHK_Val( xyr_BuoyM(imin:imax,jmin:jmax,kmin:kmax), xyr_BuoyM2(imin:imax,jmin:jmax,kmin:kmax), info2 ) info2 = "VelZ, BuoyT: " call CHK_Val( xyr_BuoyD(imin:imax,jmin:jmax,kmin:kmax), xyr_BuoyD2(imin:imax,jmin:jmax,kmin:kmax), info2 ) end subroutine BuoyancyLong_xyr subroutine VelDivC2( & & pyz_VelX, xqz_VelY, xyr_VelZ, & !(IN) & xyz_VelDiv ) !(OUT) implicit none real(DP), intent(in) :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_VelZ(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyz_VelDiv(imin:imax,jmin:jmax,kmin:kmax) integer :: i, j, k do k = kmin + 1, kmax do j = jmin + 1, jmax do i = imin + 1, imax xyz_VelDiv(i,j,k) = & & + ( & & pyz_VelX(i,j,k) & & - pyz_VelX(i-1,j,k) & & ) / dx & & + ( & & xqz_VelY(i,j,k) & & - xqz_VelY(i,j-1,k) & & ) / dy & & + ( & & xyr_VelZ(i,j,k) & & - xyr_VelZ(i,j,k-1) & & ) / dz end do end do end do ! 不定を作らないために ! ! xyz_VelDiv(imin,:,:) = 0.0d0 ! xyz_VelDiv(:,jmin,:) = 0.0d0 ! xyz_VelDiv(:,:,kmin) = 0.0d0 xyz_VelDiv(imin,:,:) = 1.0d10 xyz_VelDiv(:,jmin,:) = 1.0d10 xyz_VelDiv(:,:,kmin) = 1.0d10 end subroutine VelDivC2 subroutine PGrad_HE( & & xyz_VelDiv, xyz_Exner, & !(IN) & pyz_PGrad, pyz_SWF, & !(OUT) & xqz_PGrad, xqz_SWF ) !(OUT) implicit none real(DP), intent(in) :: xyz_VelDiv(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: pyz_PGrad(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xqz_PGrad(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: pyz_SWF(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xqz_SWF(imin:imax,jmin:jmax,kmin:kmax) integer :: i, j, k !------------------------------------------------------------------ ! X 方向 ! 音波減衰項 ! do k = kmin, kmax do j = jmin, jmax do i = imin, imax-1 pyz_SWF(i,j,k) = & & Alpha & & * ( & & xyz_VelDiv(i+1,j,k) & & - xyz_VelDiv(i,j,k) & & ) / dx end do end do end do ! 圧力傾度力 ! do k = kmin, kmax do j = jmin, jmax do i = imin, imax-1 pyz_PGrad(i,j,k) = & & - CpDry * pyz_VPTempBZ(i,j,k) & & * ( & & xyz_Exner(i+1,j,k) & & - xyz_Exner(i,j,k) & & ) / dx end do end do end do ! write(*,*) CpDry ! write(*,*) pyz_VPTempBZ(1,1,:) ! 穴埋め ! pyz_SWF(imax,:,:) = 0.0d0 pyz_PGrad(imax,:,:) = 0.0d0 !------------------------------------------------------------------ ! Y 方向 ! 音波減衰項 ! do k = kmin, kmax do j = jmin, jmax-1 do i = imin, imax xqz_SWF(i,j,k) = & & Alpha & & * ( & & xyz_VelDiv(i,j+1,k) & & - xyz_VelDiv(i,j,k) & & ) / dy end do end do end do ! 圧力傾度力 ! do k = kmin, kmax do j = jmin, jmax-1 do i = imin, imax xqz_PGrad(i,j,k) = & & - CpDry * xqz_VPTempBZ(i,j,k) & & * ( & & xyz_Exner(i,j+1,k) & & - xyz_Exner(i,j,k) & & ) /dy end do end do end do ! 穴埋め ! xqz_SWF(:,jmax,:) = 0.0d0 xqz_PGrad(:,jmax,:) = 0.0d0 end subroutine PGrad_HE subroutine PGrad_VI( & & Alpha, beta, dz, CpDry, xyr_VPTempBZ, & & imin, imax, jmin, jmax, kmin, kmax, & & xyz_VelDiv, xyz_ExnerNs, xyz_ExnerAs, & !(IN) & xyr_PGrad, xyr_SWF ) !(OUT) implicit none real(DP), intent(in) :: Alpha, beta, dz, CpDry integer, intent(in) :: imin, imax, jmin, jmax, kmin, kmax real(DP), intent(in) :: xyr_VPTempBZ(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_VelDiv(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_ExnerNs(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_ExnerAs(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyr_PGrad(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(out) :: xyr_SWF(imin:imax,jmin:jmax,kmin:kmax) integer :: i, j, k ! 音波減衰項 ! do k = kmin, kmax-1 do j = jmin, jmax do i = imin, imax xyr_SWF(i,j,k) = & & + Alpha & & * ( & & xyz_VelDiv(i,j,k+1) & & - xyz_VelDiv(i,j,k) & & ) / dz end do end do end do ! 圧力傾度力 ! do k = kmin, kmax-1 do j = jmin, jmax do i = imin, imax xyr_PGrad(i,j,k) = & & - CpDry * xyr_VPTempBZ(i,j,k) & & * ( & & beta & & * ( & & xyz_ExnerAs(i,j,k+1) & & - xyz_ExnerAs(i,j,k) & & ) & & + (1.0d0 - beta) & & * ( & & xyz_ExnerNs(i,j,k+1) & & - xyz_ExnerNs(i,j,k) & & ) & & ) / dz end do end do end do ! 穴埋め ! xyr_PGrad(:,:,kmax) = 0.0d0 xyr_SWF(:,:,kmax) = 0.0d0 end subroutine PGrad_VI subroutine Dynamics_Tendency_Output 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='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 subroutine Dynamics_Km_forcing( & & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & ! (in) & xyz_KmBl, xyz_KmNl, & ! (in) & xyz_DKmDtNl & ! (inout) & ) implicit none real(DP), intent(in) :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax) 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) :: xyz_DKmDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_Adv2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax) character(15) :: info2 ! 移流および数値拡散 ! call AdvC4_nDiff_xyz( & & xyz_KmBl, xyz_KmNl, & !(IN) & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & !(IN) & xyz_Adv, xyz_nDiff ) !(OUT) xyz_Adv2 = & & - xyz_avr_pyz(pyz_VelXNl * pyz_c4dx_xyz(xyz_KmNl)) & & - xyz_avr_xqz(xqz_VelYNl * xqz_c4dy_xyz(xyz_KmNl)) & & - xyz_avr_xyr(xyr_VelZNl * xyr_c4dz_xyz(xyz_KmNl)) xyz_nDiff2 = & & - NuHm * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyz_KmBl ))))) & & - NuHm * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_KmBl ))))) & & - NuVm * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyz_KmBl ))))) info2 = "Km, Adv: " call CHK_Val( xyz_Adv(imin:imax,jmin:jmax,kmin:kmax), xyz_Adv2(imin:imax,jmin:jmax,kmin:kmax), info2 ) info2 = "Km, nDiff: " call CHK_Val( xyz_nDiff(imin:imax,jmin:jmax,kmin:kmax), xyz_nDiff2(imin:imax,jmin:jmax,kmin:kmax), info2 ) xyz_DKmDtNl = xyz_DKmDtNl + xyz_Adv + xyz_nDiff call HistoryAutoPut(TimeN, 'KmAdv', xyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'KmNDiff', xyz_nDiff(1:nx,1:ny,1:nz)) end subroutine Dynamics_Km_forcing subroutine CHK_Val(xyz_Var1, xyz_Var2, info) real(DP), intent(in) :: xyz_Var1(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_Var2(imin:imax,jmin:jmax,kmin:kmax) character(15), intent(in) :: info integer :: min1(3), min2(3), max1(3), max2(3) real(DP) :: tmp1(nx,ny,nz), tmp2(nx,ny,nz) real(DP), parameter :: EPS = 1.0d-12 tmp1 = xyz_Var1(1:nx,1:ny,1:nz) - xyz_Var2(1:nx,1:ny,1:nz) tmp2 = (xyz_Var1(1:nx,1:ny,1:nz) - xyz_Var2(1:nx,1:ny,1:nz)) / (xyz_Var2(1:nx,1:ny,1:nz) + 1.0d-60) min1 = minloc(tmp1); min2 = minloc(tmp2) max1 = maxloc(tmp1); max2 = maxloc(tmp2) if ( abs(minval( tmp1 )) > EPS & & .OR. abs(maxval( tmp1 )) > EPS & & .OR. abs(minval( tmp2 )) > EPS & & .OR. abs(maxval( tmp2 )) > EPS ) then write(*,*) "***", info write(*,*) minval( tmp1 ), maxval( tmp1 ), minval( tmp2 ), maxval( tmp2 ) write(*,*) min1, tmp1(min1(1),min1(2),min1(3)), xyz_Var1(min1(1),min1(2),min1(3)), xyz_Var2(min1(1),min1(2),min1(3)) write(*,*) max1, tmp1(max1(1),max1(2),max1(3)), xyz_Var1(max1(1),max1(2),max1(3)), xyz_Var2(max1(1),max1(2),max1(3)) write(*,*) min2, tmp2(min2(1),min2(2),min2(3)), xyz_Var1(min2(1),min2(2),min2(3)), xyz_Var2(min2(1),min2(2),min2(3)) write(*,*) max2, tmp2(max2(1),max2(2),max2(3)), xyz_Var1(max2(1),max2(2),max2(3)), xyz_Var2(max2(1),max2(2),max2(3)) end if end subroutine CHK_Val subroutine CHK2_Val(xyz_Var1, xyz_Var2, info) real(DP), intent(in) :: xyz_Var1(nx,ny,nz) real(DP), intent(in) :: xyz_Var2(nx,ny,nz) character(15), intent(in) :: info real(DP), parameter :: EPS = 1.0d-12 integer :: min1(3), min2(3), max1(3), max2(3) real(DP) :: tmp1(nx,ny,nz), tmp2(nx,ny,nz) tmp1 = xyz_Var1(1:nx,1:ny,1:nz) - xyz_Var2(1:nx,1:ny,1:nz) tmp2 = (xyz_Var1(1:nx,1:ny,1:nz) - xyz_Var2(1:nx,1:ny,1:nz)) / (xyz_Var2(1:nx,1:ny,1:nz) + 1.0d-60) min1 = minloc(tmp1); min2 = minloc(tmp2) max1 = maxloc(tmp1); max2 = maxloc(tmp2) if ( abs(minval( tmp1 )) > EPS & & .OR. abs(maxval( tmp1 )) > EPS & & .OR. abs(minval( tmp2 )) > EPS & & .OR. abs(maxval( tmp2 )) > EPS ) then write(*,*) "***", info write(*,*) minval( tmp1 ), maxval( tmp1 ), minval( tmp2 ), maxval( tmp2 ) write(*,*) min1, tmp1(min1(1),min1(2),min1(3)), xyz_Var1(min1(1),min1(2),min1(3)), xyz_Var2(min1(1),min1(2),min1(3)) write(*,*) max1, tmp1(max1(1),max1(2),max1(3)), xyz_Var1(max1(1),max1(2),max1(3)), xyz_Var2(max1(1),max1(2),max1(3)) write(*,*) min2, tmp2(min2(1),min2(2),min2(3)), xyz_Var1(min2(1),min2(2),min2(3)), xyz_Var2(min2(1),min2(2),min2(3)) write(*,*) max2, tmp2(max2(1),max2(2),max2(3)), xyz_Var1(max2(1),max2(2),max2(3)), xyz_Var2(max2(1),max2(2),max2(3)) end if end subroutine CHK2_Val end module DynamicsHEVI