!= Module DynamicsHEVI ! ! Authors:: 杉山耕一朗(SUGIYAMA Ko-ichiro), ODAKA Masatsugu ! Version:: $Id: dynamicshevi.f90,v 1.22 2015/02/19 02:17:23 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 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 ! 物質数 use constants,only: CpDry, &! 乾燥成分の比熱 & CvDry, &! 乾燥成分の定積比熱 & GasRDry, &! 乾燥成分の気体定数 & MolWtDry, &! 乾燥成分の分子量 & Grav ! 重力加速度 use composition, only: SpcWetSymbol, &! & GasNum, &! & IdxG, &! & MolWtWet ! 湿潤成分の分子量 use timeset, only: DelTimeShort, DelTimeLong, TimeN use axesset, only: x_dx, y_dy, z_dz, r_dz, &! 格子間隔 & xyz_avr_pyz, xyr_avr_pyr, xqz_avr_pqz, & & pyz_avr_xyz, pyr_avr_xyr, pqz_avr_xqz, & & xyz_avr_xqz, pyz_avr_pqz, xyr_avr_xqr, & & xqz_avr_xyz, pqz_avr_pyz, xqr_avr_xyr, & & xyz_avr_xyr, pyz_avr_pyr, xqz_avr_xqr, & & xyr_avr_xyz, pyr_avr_pyz, xqr_avr_xqz use basicset, only: xyz_VelSoundBZ, &!基本場の音速 & xyz_DensBZ, &!基本場の密度 & xyz_PTempBZ, &!基本場の温位 & xyz_VPTempBZ, &!基本場の温位 & xyz_ExnerBZ, & & xyzf_QMixBZ, & & xyr_QMixBZPerMolWt, & & xyr_QMixBZ, xyz_EffMolWtBZ use xyz_deriv_module use xyz_deriv_c4_module, only: pyz_c4dx_xyz, xqz_c4dy_xyz, xyr_c4dz_xyz, & & xyz_c4dx_pyz, pqz_c4dy_pyz, pyr_c4dz_pyz, & & pqz_c4dx_xqz, xyz_c4dy_xqz, xqr_c4dz_xqz, & & pyr_c4dx_xyr, xqr_c4dy_xyr, xyz_c4dz_xyr 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 real(DP), allocatable, save :: xyz_F1BZ(:,:,:) !係数行列の計算に用いる配列 real(DP), allocatable, save :: xyr_F2BZ(:,:,:) !係数行列の計算に用いる配列 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 :: A(:) !係数行列の対角成分 real(DP), allocatable, save :: B(:) !係数行列の上三角部分 real(DP), allocatable, save :: C(:) !係数行列の下三角部分 ! real(DP), allocatable, save :: AU2(:,:) !LU 分解の結果 U (2 次元配列) real(DP), allocatable, save :: AL1(:) !LU 分解の結果 L (1 次元配列) ! real(DP), allocatable, save :: AL2(:,:) !LU 分解の結果 L (2 次元配列) integer, allocatable, save :: IP(:) !部分ピボット交換の情報を格納 real(DP), save :: Alpha = 0.0d0 !音波減衰項の減衰係数 real(DP), save :: NuHh4 = 0.0d0 !数値粘性の係数 (水平方向) real(DP), save :: NuVh4 = 0.0d0 !数値粘性の係数 (鉛直方向) real(DP), save :: NuHm4 = 0.0d0 !数値粘性の係数 (水平方向) real(DP), save :: NuVm4 = 0.0d0 !数値粘性の係数 (鉛直方向) real(DP), save :: NuHh2 = 0.0d0 !数値粘性の係数 (水平方向) real(DP), save :: NuVh2 = 0.0d0 !数値粘性の係数 (鉛直方向) real(DP), save :: NuHm2 = 0.0d0 !数値粘性の係数 (水平方向) real(DP), save :: NuVm2 = 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_Km_forcing public Dynamics_Long_forcing public Dynamics_Short_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) :: AlphaNDiff2 = 0.0d0 !2次の数値拡散の係数. CReSS マニュアルより real(DP) :: NDiffRatio = 1.0d0 !速度に対する粘性を上げる場合は数字を 1 以上にする. integer :: unit !装置番号 integer :: l NAMELIST /Dynamics_nml/ & & AlphaSound, AlphaNDiff, AlphaNDiff2, 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 = minval(x_dx) DelYMin = minval(y_dy) DelZMin = minval(z_dz) if (jmin == jmax) then Alpha = AlphaSound * ( Min(DelXMin, DelZMin) ** 2.0d0 ) / DelTimeShort else Alpha = AlphaSound * ( Min(DelXMin, DelYMin, DelZMin) ** 2.0d0 ) / DelTimeShort end if if (jmin == jmax) then NuHh4 = AlphaNDiff * ( DelXMin ** 4.0d0 ) / (2.0d0 * DelTimeLong) else NuHh4 = AlphaNDiff * ( SQRT(DelXMin*DelYMin) ** 4.0d0 ) / (2.0d0 * DelTimeLong) end if NuVh4 = AlphaNDiff * ( DelZMin ** 4.0d0 ) / (2.0d0 * DelTimeLong) NuHm4 = NuHh4 * NDiffRatio NuVm4 = NuVh4 * NDiffRatio if (jmin == jmax) then NuHh2 = AlphaNDiff2 * ( DelXMin ** 2.0d0 ) / (2.0d0 * DelTimeLong) else NuHh2 = AlphaNDiff2 * ( SQRT(DelXMin*DelYMin) ** 2.0d0 ) / (2.0d0 * DelTimeLong) end if NuVh2 = AlphaNDiff2 * ( DelZMin ** 2.0d0 ) / (2.0d0 * DelTimeLong) NuHm2 = NuHh2 * NDiffRatio NuVm2 = NuVh2 * NDiffRatio if (myrank == 0) then if ( AlphaNDiff > 0.0d0 .AND. AlphaNDiff2 > 0.0d0 ) then write(*,*) "CHECK! AlphaNDiff == 0.0d0 .OR. AlphaNDiff2 == 0.0d0" stop end if call MessageNotify( "M", module_name, "Alpha = %f", d=(/Alpha/) ) call MessageNotify( "M", module_name, "NuHh = %f", d=(/NuHh4/) ) call MessageNotify( "M", module_name, "NuVh = %f", d=(/NuVh4/) ) call MessageNotify( "M", module_name, "NuHm = %f", d=(/NuHm4/) ) call MessageNotify( "M", module_name, "NuVm = %f", d=(/NuVm4/) ) call MessageNotify( "M", module_name, "NuHh = %f", d=(/NuHh2/) ) call MessageNotify( "M", module_name, "NuVh = %f", d=(/NuVh2/) ) call MessageNotify( "M", module_name, "NuHm = %f", d=(/NuHm2/) ) call MessageNotify( "M", module_name, "NuVm = %f", d=(/NuVm2/) ) 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 DynamicsVI_init() 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='PTempNDiff2',& & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of potential temperature (2 order)',& & 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='ExnerNDiff2',& & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of exner function (2 order)',& & 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') call HistoryAutoAddVariable( & & varname='CDensNDiff2',& & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of cloud density (2 order)',& & units='kg.m-3.s-1', & & xtype='float') do l = 1, ncmax call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(l))//'_Adv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of ' & & //trim(SpcWetSymbol(l))//' mixing ratio', & & units='kg.kg-1.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(l))//'_NDiff', & & dims=(/'x','y','z','t'/), & & longname='Diffusion term of ' & & //trim(SpcWetSymbol(l))//' mixing ratio', & & units='kg.kg-1.s-1', & & xtype='float') call HistoryAutoAddVariable( & & varname=trim(SpcWetSymbol(l))//'_NDiff2', & & dims=(/'x','y','z','t'/), & & longname='Diffusion term of ' & & //trim(SpcWetSymbol(l))//' mixing ratio (2 order)', & & 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='VelXNDiff2',& & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of velocity (x) (2 order)',& & 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='VelXTndNs', & ! & dims=(/'x','y','z','t'/), & ! & longname='Velocity Tendency (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='VelYNDiff2',& & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of velocity (y) (2 order)',& & 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='VelYTndNs', & ! & dims=(/'x','y','z','t'/), & ! & longname='Velocity Tendency (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='VelZNDiff2',& & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of Velocity (z) (2 order)',& & 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') call HistoryAutoAddVariable( & & varname='KmNDiff2', & & dims=(/'x','y','z','t'/), & & longname='Diffusion term of Km', & & units='s-1', & & xtype='float') end subroutine Dynamics_Init 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_Orig(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_NDiff4(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_NDiff2(imin:imax,jmin:jmax,kmin:kmax) !---------------------------------- ! 拡散係数 ! ! Initialize ! xyz_Orig = xyz_DKmDtNl ! Advection term ! xyz_Adv = & & - 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)) ! Numerical diffusion term ! xyz_NDiff4 = & & - NuHm4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyz_KmBl ))))) & & - NuHm4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_KmBl ))))) & & - NuVm4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyz_KmBl ))))) xyz_NDiff2 = & & + NuHm2 * (xyz_dx_pyz(pyz_dx_xyz( xyz_KmBl ))) & & + NuHm2 * (xyz_dy_xqz(xqz_dy_xyz( xyz_KmBl ))) & & + NuVm2 * (xyz_dz_xyr(xyr_dz_xyz( xyz_KmBl ))) xyz_DKmDtNl = xyz_Orig + xyz_Adv + xyz_NDiff4 + xyz_NDiff2 call HistoryAutoPut(TimeN, 'KmAdv', xyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'KmNDiff', xyz_NDiff4(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'KmNDiff2', xyz_NDiff2(1:nx,1:ny,1:nz)) ! Set Margin ! ! call SetMargin_xyz( xyz_DKmDtNl ) end subroutine Dynamics_Km_forcing 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) ! & xyz_CDensBl, xyz_CDensNl, & ! (in) & pyz_DVelXDtNl, & ! (inout) & xqz_DVelYDtNl, & ! (inout) & xyr_DVelZDtNl, & ! (inout) & xyz_DPTempDtNl, & ! (inout) & xyzf_DQMixDtNl, & ! (inout) ! & xyz_DCDensDtNl, & ! (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), 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) :: xyz_DCDensDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_Orig(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_NDiff4(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_NDiff2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_Orig(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_NDiff4(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_NDiff2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_Orig(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_NDiff4(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_NDiff2(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_Orig(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_NDiff4(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_NDiff2(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_PTempAll(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyzf_Orig(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP) :: xyzf_Adv(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP) :: xyzf_NDiff4(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP) :: xyzf_NDiff2(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, 1:ncmax) real(DP) :: xyzf_QMixPerMolWt(imin:imax,jmin:jmax,kmin:kmax, 1:GasNum) integer :: f !------------------------------ ! tendency of cloud density ! ! initialize ! xyz_Orig = xyz_DCDensDtNl ! フラックス項の計算. 4 次精度中心差分を用いて計算 ! ! xyz_Adv = & ! & - xyz_c4dx_pyz(pyz_VelXNl * pyz_avr_xyz(xyz_CDensNl)) & ! & - xyz_c4dy_xqz(xqz_VelYNl * xqz_avr_xyz(xyz_CDensNl)) & ! & - xyz_c4dz_xyr(xyr_VelZNl * xyr_avr_xyz(xyz_CDensNl)) ! 数値粘性項の計算 ! xyz_NDiff4 = & ! & - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_CDensBl))))) & ! & - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_CDensBl))))) & ! & - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_CDensBl))))) ! xyz_NDiff2 = & ! & + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz(xyz_CDensBl))) & ! & + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz(xyz_CDensBl))) & ! & + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz(xyz_CDensBl))) ! xyz_DCDensDtNl = xyz_Orig + xyz_Adv + xyz_NDiff4 + xyz_NDiff2 ! call HistoryAutoPut(TimeN, 'CDensAdv', xyz_Adv(1:nx,1:ny,1:nz)) ! call HistoryAutoPut(TimeN, 'CDensNDiff', xyz_NDiff4(1:nx,1:ny,1:nz)) ! call HistoryAutoPut(TimeN, 'CDensNDiff2', xyz_NDiff2(1:nx,1:ny,1:nz)) ! tendency of potential temperature ! ! initialize xyz_Orig = xyz_DPTempDtNl xyz_PTempAll = xyz_PTempNl + xyz_PTempBZ ! Advection term ! xyz_AdvScalar( xyz_PTempNl + xyz_PTempBZ, pyz_VelXNl, pyz_VelXNl, xyr_VelZNl) ! xyz_Adv = & & - 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)) ! numerical diffusion term ! xyz_Num = xyz_NumDiffScalar( xyz_PTempBl) ! xyz_NDiff4 = & & - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyz_PTempBl ))))) & & - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyz_PTempBl ))))) & & - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyz_PTempBl ))))) xyz_NDiff2 = & & + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz( xyz_PTempBl ))) & & + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz( xyz_PTempBl ))) & & + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz( xyz_PTempBl ))) ! sum ! xyz_DPTempDtNl = xyz_Orig + xyz_Adv + xyz_NDiff4 + xyz_NDiff2 ! output ! call HistoryAutoPut(TimeN, 'PTempAdv', xyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'PTempNDiff', xyz_NDiff4(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'PTempNDiff2', xyz_NDiff2(1:nx,1:ny,1:nz)) xyz_PTempAl = xyz_PTempBl + (2.0d0 * DelTimeLong) * xyz_DPTempDtNl ! Set Margin ! call SetMargin_xyz(xyz_PTempAl) !------------------------------ ! tendency of mixing ratio ! ! initialize xyzf_Orig = xyzf_DQMixDtNl xyzf_QMixAll = xyzf_QMixNl + xyzf_QMixBZ do f = 1, ncmax ! Advection term !xyzf_Adv = xyzf_AdvScalar(xyzf_QMixNl + xyzf_QMixBZ, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) ! xyzf_Adv(:,:,:,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))) ! numerical diffusion term ! xyzf_Diff = xyzf_NumDiffScalar(xyzf_QMixBl) ! xyzf_NDiff4(:,:,:,f) = & & - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz( xyzf_QMixBl(:,:,:,f) ))))) & & - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz( xyzf_QMixBl(:,:,:,f) ))))) & & - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz( xyzf_QMixBl(:,:,:,f) ))))) xyzf_NDiff2(:,:,:,f) = & & + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz( xyzf_QMixBl(:,:,:,f) ))) & & + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz( xyzf_QMixBl(:,:,:,f) ))) & & + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz( xyzf_QMixBl(:,:,:,f) ))) end do ! sum ! xyzf_DQMixDtNl = xyzf_Orig + xyzf_Adv + xyzf_NDiff4 + xyzf_NDiff2 ! 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_NDiff4(1:nx,1:ny,1:nz,f)) call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_NDiff2', xyzf_NDiff2(1:nx,1:ny,1:nz,f)) end do ! time integration ! xyzf_QMixAl = xyzf_QMixBl + (2.0d0 * DelTimeLong) * xyzf_DQMixDtNl ! Set Margin ! call SetMargin_xyzf(xyzf_QMixAl) ! 負の値を埋める ! call FillNegativeQMix(xyzf_QMixAl) ! Set Margin ! call SetMargin_xyzf(xyzf_QMixAl) !------------------------------ ! tendency of VelX ! ! initializa ! pyz_Orig = pyz_DVelXDtNl ! Advection term !pyz_Adv = pyz_AdvVelX(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) ! pyz_Adv = & & - 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 ) ) ! Numerical diffusion term !pyz_Diff = pyz_NumDiffVelX(pyz_VelXBl) pyz_NDiff4 = & & - NuHm4 * ( pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz( pyz_VelXBl ))))) & & - NuHm4 * ( pyz_dy_pqz(pqz_dy_pyz(pyz_dy_pqz(pqz_dy_pyz( pyz_VelXBl ))))) & & - NuVm4 * ( pyz_dz_pyr(pyr_dz_pyz(pyz_dz_pyr(pyr_dz_pyz( pyz_VelXBl ))))) pyz_NDiff2 = & & + NuHm2 * ( pyz_dx_xyz(xyz_dx_pyz( pyz_VelXBl ))) & & + NuHm2 * ( pyz_dy_pqz(pqz_dy_pyz( pyz_VelXBl ))) & & + NuVm2 * ( pyz_dz_pyr(pyr_dz_pyz( pyz_VelXBl ))) ! sum ! pyz_DVelXDtNl = pyz_Orig + pyz_Adv + pyz_NDiff4 + pyz_NDiff2 call HistoryAutoPut(TimeN, 'VelXAdv', pyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelXNDiff', pyz_NDiff4(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelXNDiff2', pyz_NDiff2(1:nx,1:ny,1:nz)) !------------------------------ ! tendency of VelY ! ! ininitalize xqz_Orig = xqz_DVelYDtNl ! Advection term ! xqz_Adv = xqz_AdvVelY(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) ! xqz_Adv = & & - 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 ) ) ! Numerical diffusion term ! xqz_Diff = xqz_NumDiffVelY(xqz_VelYBl) ! xqz_NDiff4 = & & - NuHm4 * ( xqz_dx_pqz(pqz_dx_xqz(xqz_dx_pqz(pqz_dx_xqz( xqz_VelYBl ))))) & & - NuHm4 * ( xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz( xqz_VelYBl ))))) & & - NuVm4 * ( xqz_dz_xqr(xqr_dz_xqz(xqz_dz_xqr(xqr_dz_xqz( xqz_VelYBl ))))) xqz_NDiff2 = & & + NuHm2 * ( xqz_dx_pqz(pqz_dx_xqz( xqz_VelYBl ))) & & + NuHm2 * ( xqz_dy_xyz(xyz_dy_xqz( xqz_VelYBl ))) & & + NuVm2 * ( xqz_dz_xqr(xqr_dz_xqz( xqz_VelYBl ))) ! sum ! xqz_DVelYDtNl = xqz_Orig + xqz_Adv + xqz_NDiff4 + xqz_NDiff2 call HistoryAutoPut(TimeN, 'VelYAdv', xqz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelYNDiff', xqz_NDiff4(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelYNDiff2', xqz_NDiff2(1:nx,1:ny,1:nz)) !------------------------------ ! tendency of VelZ ! ! Initialization ! xyr_Orig = xyr_DVelZDtNl do f = 1, GasNum xyzf_QMixPerMolWt(:,:,:,f) = xyzf_QMixNl(:,:,:,IdxG(f)) / MolWtWet(IdxG(f)) end do ! Advection term ! xyr_Adv = xyr_AdvVelZ(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) ! xyr_Adv = & & - 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 ) ) ! Numerical diffusion term !xyr_Diff = xyr_NumDiffVelZ(xyr_VelZBl) ! xyr_NDiff4 = & & - NuHm4 * ( xyr_dx_pyr(pyr_dx_xyr(xyr_dx_pyr(pyr_dx_xyr( xyr_VelZBl ))))) & & - NuHm4 * ( xyr_dy_xqr(xqr_dy_xyr(xyr_dy_xqr(xqr_dy_xyr( xyr_VelZBl ))))) & & - NuVm4 * ( xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr( xyr_VelZBl ))))) xyr_NDiff2 = & & + NuHm2 * ( xyr_dx_pyr(pyr_dx_xyr( xyr_VelZBl ))) & & + NuHm2 * ( xyr_dy_xqr(xqr_dy_xyr( xyr_VelZBl ))) & & + NuVm2 * ( xyr_dz_xyz(xyz_dz_xyr( xyr_VelZBl ))) ! Buoyancy due to temperature disturbunce !xyr_BuoyT = xyr_Buoy(xyz_PTempNl) ! xyr_BuoyT = Grav * xyr_avr_xyz( xyz_PTempNl / xyz_PTempBZ) ! Buoyancy due to molecular weight ! xyr_BuoyM = & & + 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 ) ! Buoyancy due to loading ! xyr_BuoyD = & & - Grav * xyr_avr_xyz( sum(xyzf_QMixNl(:,:,:,GasNum+1:ncmax), 4) ) & & / ( 1.0d0 + xyr_QMixBZ ) ! sum ! xyr_DVelZDtNl = xyr_Orig + xyr_Adv + xyr_NDiff4 + xyr_NDiff2 & & + 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_NDiff4(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelZNDiff2', xyr_NDiff2(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)) ! Set Margin ! ! call SetMargin_pyz( pyz_DVelXDtNl ) ! call SetMargin_xqz( xqz_DVelYDtNl ) ! call SetMargin_xyr( xyr_DVelZDtNl ) ! call SetMargin_xyz( xyz_DPTempDtNl ) ! call SetMargin_xyz( xyz_DCDensDtNl ) ! call SetMargin_xyzf(xyzf_DQMixDtNl ) end subroutine Dynamics_Long_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) & ) 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_DExnerDtNs(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) :: xyz_VelLaplaNs(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_VelXPGrad(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_VelYPGrad(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_VelZPGrad(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: pyz_VelXSWF(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xqz_VelYSWF(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyr_VelZSWF(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_ExnerAll(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_Adv(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_NDiff4(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_NDiff2(imin:imax,jmin:jmax,kmin:kmax) ! initialize: Divergence of velocity ! xyz_VelDivNs = & & xyz_dx_pyz( pyz_VelXNs ) & & + xyz_dy_xqz( xqz_VelYNs ) & & + xyz_dz_xyr( xyr_VelZNs ) ! call HistoryAutoPut(TimeN, 'VelDiv', xyz_VelDivNs(1:nx,1:ny,1:nz)) ! xyz_VelLaplaNs = pyz_dx_xyz( xyz_VelDivNs ) + xqz_dy_xyz( xyz_VelDivNs ) + xyr_dz_xyz( xyz_VelDivNs ) ! call HistoryAutoPut(TimeN, 'VelLapla', xyz_VelLaplaNs(1:nx,1:ny,1:nz)) !-------------------------------------- ! VelX ! pyz_VelXSWF = Alpha * pyz_dx_xyz( xyz_VelDivNs ) pyz_VelXPGrad = - pyz_avr_xyz( CpDry * xyz_VPTempBZ ) * pyz_dx_xyz( xyz_ExnerNs ) & & + pyz_VelXSWF pyz_DVelXDtNs = pyz_VelXPGrad ! Time integration ! pyz_VelXAs = pyz_VelXNs + DelTimeShort * (pyz_DVelXDtNl + pyz_DVelXDtNs) call HistoryAutoPut(TimeN, 'VelXPGrad', pyz_VelXPGrad(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelXSWF', pyz_VelXSWF(1:nx,1:ny,1:nz)) ! call HistoryAutoPut(TimeN, 'VelXTndNs', pyz_DVelXDtNs(1:nx,1:ny,1:nz)) ! Set Margin ! call SetMargin_pyz( pyz_VelXAs ) ! (inout) !-------------------------------------- ! VelY ! xqz_VelYSWF = Alpha * xqz_dy_xyz( xyz_VelDivNs ) xqz_VelYPGrad = - xqz_avr_xyz( CpDry * xyz_VPTempBZ ) * xqz_dy_xyz( xyz_ExnerNs ) & & + xqz_VelYSWF xqz_DVelYDtNs = xqz_VelYPGrad ! Time integration ! xqz_VelYAs = xqz_VelYNs + DelTimeShort * (xqz_DVelYDtNl + xqz_DVelYDtNs) call HistoryAutoPut(TimeN, 'VelYPGrad', xqz_VelYPGrad(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelYSWF', xqz_VelYSWF(1:nx,1:ny,1:nz)) ! call HistoryAutoPut(TimeN, 'VelYTndNs', xqz_DVelYDtNs(1:nx,1:ny,1:nz)) ! Set Margin ! call SetMargin_xqz( xqz_VelYAs ) ! (inout) !-------------------------------------- ! Exner function ! ! フラックス項の計算. 4 次精度中心差分を用いて計算 ! xyz_Adv = & & - xyz_avr_pyz(pyz_VelXNs * pyz_c4dx_xyz(xyz_ExnerNs)) & & - xyz_avr_xqz(xqz_VelYNs * xqz_c4dy_xyz(xyz_ExnerNs)) & & - xyz_avr_xyr(xyr_VelZNs * xyr_c4dz_xyz(xyz_ExnerNs)) !& ! & + CpDry / CvDry * GasRDry * xyz_ExnerNs * xyz_VelDivNs ! 数値粘性項の計算 xyz_NDiff4 = & & - NuHh4 * (xyz_dx_pyz(pyz_dx_xyz(xyz_dx_pyz(pyz_dx_xyz(xyz_ExnerNs))))) & & - NuHh4 * (xyz_dy_xqz(xqz_dy_xyz(xyz_dy_xqz(xqz_dy_xyz(xyz_ExnerNs))))) & & - NuVh4 * (xyz_dz_xyr(xyr_dz_xyz(xyz_dz_xyr(xyr_dz_xyz(xyz_ExnerNs))))) xyz_NDiff2 = & & + NuHh2 * (xyz_dx_pyz(pyz_dx_xyz(xyz_ExnerNs))) & & + NuHh2 * (xyz_dy_xqz(xqz_dy_xyz(xyz_ExnerNs))) & & + NuVh2 * (xyz_dz_xyr(xyr_dz_xyz(xyz_ExnerNs))) ! 短い時間ステップでのtendency xyz_DExnerDtNs = xyz_DExnerDtNs + xyz_Adv + xyz_NDiff4 + xyz_NDiff2 !!!!!!!!!!!!!!!!!!!!!!! ! ! test ! !!!!!!!!!!!!!!!!!!!!!!! xyz_DExnerDtNs = 0.0d0 call HistoryAutoPut(TimeN, 'ExnerAdv', xyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'ExnerNDiff', xyz_NDiff4(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'ExnerNDiff2', xyz_NDiff2(1:nx,1:ny,1:nz)) xyz_ExnerAs = xyz_Exner( & & pyz_VelXAs, & & xqz_VelYAs, & & xyr_VelZNs, & & xyz_VelDivNs, & & xyz_ExnerNs, & & xyr_DVelZDtNl, & & xyz_DExnerDtNl, & & xyz_DExnerDtNs & & ) ! Set Margin ! call SetMargin_xyz( xyz_ExnerAs ) ! (inout) ! write(*,*) "+++ Exner +++", xyz_ExnerAs(imin:imax,1,kmin:kmax) !-------------------------------------- ! VelZ ! xyr_VelZSWF = Alpha * xyr_dz_xyz( xyz_VelDivNs ) xyr_VelZPGrad = & & - xyr_avr_xyz(CpDry * xyz_VPTempBZ ) & & * ( & & beta * xyr_dz_xyz( xyz_ExnerAs ) & & + (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerNs ) & & ) & & + xyr_VelZSWF xyr_DVelZDtNs = xyr_VelZPGrad ! Time integration ! xyr_VelZAs = xyr_VelZNs + DelTimeShort * (xyr_DVelZDtNl + xyr_DVelZDtNs) call HistoryAutoPut(TimeN, 'VelZPGrad', xyr_VelZPGrad(1:nx,1:ny,1:nz)) call HistoryAutoPut(TimeN, 'VelZSWF', xyr_VelZSWF(1:nx,1:ny,1:nz)) ! call HistoryAutoPut(TimeN, 'VelZTndNs', xyr_DVelZDtNs(1:nx,1:ny,1:nz)) ! Set Margin ! call SetMargin_xyr( xyr_VelZAs ) ! (inout) end subroutine Dynamics_Short_forcing !!!--------------------------------------------------------------------!!! subroutine DynamicsVI_init() ! !エクスナー関数を陰解法で解く際に必要となる, 係数行列の要素を決め, !LU 分解を行う. ! !暗黙の型宣言禁止 implicit none real(DP) :: DTS ! 短い時間格子 DTS = DelTimeShort !配列の割り付け allocate( A(1:nz) ) allocate( B(1+1:nz) ) allocate( C(1:nz-1) ) allocate( xyz_F1BZ(imin:imax,jmin:jmax,kmin:kmax) ) allocate( xyr_F2BZ(imin:imax,jmin:jmax,kmin:kmax) ) !---------------------------------------------------------------- ! 係数行列と共通して利用される配列の値を決める !---------------------------------------------------------------- !係数行列の計算 ! A, B, C を求める際, F1BZ と F2BZ は X 方向に一様なので. ! nx, ny の値で代表させることとした. xyz_F1BZ = & & ( xyz_VelSoundBZ ** 2.0d0 ) & & / (CpDry * xyz_DensBZ * (xyz_VPTempBZ ** 2.0d0)) xyr_F2BZ = & & xyr_avr_xyz( & & CpDry * xyz_DensBZ * ( xyz_VPTempBZ ** 2.0d0 ) & & ) A(1+1: nz-1) = & & (beta ** 2.0d0) & & * xyz_F1BZ(nx,ny,1+1: nz-1) & & * (DTS ** 2.0d0) & & * ( & & xyr_F2BZ(nx,ny,1+1: nz-1) & & / r_dz(1+1: nz-1) & & + xyr_F2BZ(nx,ny,1 : nz-2) & & / r_dz(1: nz-2) & & ) & & / z_dz(1+1: nz-1) & & + 1.0d0 A(1) = & & (beta ** 2.0d0) & & * xyz_F1BZ(nx,ny,1) & & * xyr_F2BZ(nx,ny,1) & & / r_dz(1) & & * (DTS ** 2.0d0) & & / z_dz(1) & & + 1.0d0 A(nz) = & & (beta ** 2.0d0) & & * xyz_F1BZ(nx,ny,nz) & & * xyr_F2BZ(nx,ny,nz-1) & & / r_dz(nz-1) & & * (DTS ** 2.0d0) & & / z_dz(nz) & & + 1.0d0 B(1+1:nz) = & & - (beta ** 2.0d0) & & * xyz_F1BZ(nx,ny,1:nz-1) & & * xyr_F2BZ(nx,ny,1:nz-1) & & * (DTS ** 2.0d0) & & / ( r_dz(1:nz-1) * z_dz(1:nz-1) ) C(1: nz-1) = & & - ( beta ** 2.0d0 ) & & * xyz_F1BZ(nx,ny,1+1:nz) & & * xyr_F2BZ(nx,ny,1 :nz-1) & & * (DTS ** 2.0d0) & & / ( r_dz(1:nz-1) * z_dz(1+1:nz) ) !---------------------------------------------------------------- ! 係数行列を LU 分解 !---------------------------------------------------------------- !配列の大きさを保管 N = nz !係数行列/改行列の次数, 整合寸法 M = nx * ny !方程式の組数 NUD = 1 !係数行列の上三角部分の帯幅 NLD = 1 !係数行列の下三角部分の帯幅 NAL = NLD !LU 分解の結果 L の整合寸法 NA = NUD + NLD + 1 !配列の割り当て ! allocate( AL1(N), AL2(NAL, N), AU2(NA, N), IP(N) ) allocate( AL1(N), IP(N) ) !LU 分解の実行 ! LAPACK の利用 call ResolvLU_Lapack( ) end subroutine DynamicsVI_init !!!--------------------------------------------------------------------!!! function xyz_Exner( & & pyz_VelXAs, & & xqz_VelYAs, & & xyr_VelZNs, & & xyz_VelDivNs, & & xyz_ExnerNs, & & xyr_DVelZDtNl, & & xyz_DExnerDtNl, & & xyz_DExnerDtNs & & ) ! !陰解法を用いたエクスナー関数の計算. ! !暗黙の型宣言禁止 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_DExnerDtNs & & (imin:imax,jmin:jmax,kmin:kmax) !Z 方向の外力項[t] real(DP), intent(in) :: xyz_ExnerNs & & (imin:imax,jmin:jmax,kmin:kmax) !無次元圧力 real(DP) :: xyz_Exner & & (imin:imax,jmin:jmax,kmin:kmax) !無次元圧力[τ+Δτ] !変数定義 real(DP) :: D1(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: D2(1:nx,1:ny,1:nz) real(DP) :: D(nx*ny,nz) real(DP) :: E(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: F(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: xyz_DivVelNs(imin:imax,jmin:jmax,kmin:kmax) real(DP) :: DTS ! 短い時間格子間隔 integer :: ix, jy, kz ! Initialize DTS = DelTimeShort xyz_Exner = 0.0d0 !行列計算のための係数 E = & & - ( 1.0d0 - beta ) * xyr_dz_xyz( xyz_ExnerNs ) & & + ( Alpha * xyr_dz_xyz( xyz_VelDivNs ) + xyr_DVelZDtNl ) & & / xyr_avr_xyz( CpDry * xyz_VPTempBZ ) F = - beta * xyz_F1BZ * 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_DExnerDtNs + xyz_DExnerDtNl ) * DTS D1 = xyz_ExnerNs & & - (1.0d0 - beta) & & * xyz_F1BZ * 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 ) & & + F D1(:,:,1) = D1(:,:,1) & & - beta * xyz_F1BZ(:,:,1) * (DTS ** 2.0d0) & & * xyr_F2BZ(:,:, 1-1) * E(:,:,1-1) & & / z_dz(1) D1(:,:,nz) = D1(:,:,nz) & & + beta * xyz_F1BZ(:,:,nz) * (DTS ** 2.0d0) & & * xyr_F2BZ(:,:,nz) * E(:,:,nz) & & / z_dz(nz) D2 = D1(1:nx,1:ny,1:nz) do kz = 1, nz do jy = 1, ny do ix = 1, nx D(ix + nx * (jy - 1), kz) = D2(ix,jy,kz) end do end do end do !----------------------------------------------------------- !連立一次方程式の解を求める !------------------------------------------------------------ !解の計算 ! LAPACK 利用 call LinSolv_Lapack( D ) !戻り値を出力 do kz = 1, nz do jy = 1, ny do ix = 1, nx xyz_Exner(ix,jy,kz) = D(ix + nx * (jy - 1 ), kz) end do end do end do end function xyz_Exner !!!--------------------------------------------------------------------!!! subroutine ResolvLU_Lapack( ) ! !実 3 項行列の LU 分解(倍精度). LAPACK 利用 ! !暗黙の型宣言禁止 implicit none !変数定義 integer :: INFO !解のコンディションチェック !変数の初期化 INFO = 0 !解行列の計算. 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 ResolvLU_Lapack !!!--------------------------------------------------------------------!!! subroutine LinSolv_Lapack( X ) ! !LU 分解された実 3 項行列の連立 1 次方程式(倍精度). LAPACK 利用 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(inout) :: X(M, N) !定数/解行列 real(DP) :: TX(N, M) !解行列を転置したもの integer :: NRHS ! integer :: INFO integer :: LDB character(1),parameter :: TRANS = 'N' !変数の初期化 NRHS = M INFO = 0 LDB = N TX = transpose( X ) !解行列の計算. LAPACK を使用. call DGTTRS(TRANS, N, NRHS, C, A, B, AL1, IP, TX, LDB, INFO) !解の出力 X = transpose( TX ) !解のコンディションをチェック. ! if (INFO /= 0) then ! call MessageNotify("Error", "lapack_linear", "INFO is not 0") ! stop ! end if end subroutine LinSolv_Lapack end module DynamicsHEVI