!= Module DynFunc_3D ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! Version:: $Id: dynamicalcore.f90,v 1.2 2011-05-04 17:55:15 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 DynamicalCore ! !陽開放を用いた力学過程の各項の計算モジュール. !具体的には以下の項を計算するための関数を格納する. ! * 移流項 ! * 浮力項 ! * 気圧傾度力項 ! !モジュール読み込み 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 方向の物理領域の上限 & nf ! 物質数 use constants,only: CpDry, &! 乾燥成分の比熱 & MolWtDry, &! 乾燥成分の分子量 & Grav ! 重力加速度 use moistset, only: SpcWetSymbol, &! & GasNum, &! & IdxG, &! & MolWtWet ! 湿潤成分の分子量 use timeset, only: DelTimeShort, DelTimeLong 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, &!基本場の温位 & xyzf_QMixBZ, & & xyr_QMixBZPerMolWt, & & xyr_QMixBZ 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 !暗黙の型宣言禁止 implicit none !属性の指定 private real(DP), save :: beta = 5.0d-1 !クランクニコルソン法なら 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 :: Nu = 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 !数値粘性の係数 (鉛直方向) !public public Dynamics_Init public DynamicsVI_Init public DynamicsLong public DynamicsShortHEVI contains subroutine Dynamics_Init(cfgfile) !暗黙の型宣言禁止 implicit none character(STRING), intent(in) :: cfgfile !NAMELIST ファイル real(DP) :: DelXMin, DelYMin, DelZMin real(DP) :: AlphaSound = 1.0e-7 !音波減衰項の係数 real(DP) :: AlphaNDiff = 1.0d-4 real(DP) :: NDiffRatio = 1.0d0 !速度に対する粘性を上げる場合は数字を 1 以上にする. integer :: unit !装置番号 integer :: l NAMELIST /dynamicalcore_nml/ AlphaSound, AlphaNDiff, NDiffRatio !ファイルオープン. 情報取得. call FileOpen(unit, file=cfgfile, mode='r') read(unit, NML=dynamicalcore_nml) close(unit) !------------------------------------------------------------------- ! 音波減衰項の減衰率 Min(DelX, DelZ) ** 2.0 に比例 ! DelXMin = minval(x_dx) DelYMin = minval(y_dy) DelZMin = minval(z_dz) Nu = AlphaSound * ( Min(DelXMin, DelZMin) ** 2.0d0 ) / DelTimeShort ! CReSS マニュアルでは, 2 次精度中心差分の場合, Alpha < 1/8 くらいが適当と述べている. ! deepconv では NuH として 500 以上の値が入っていたそうだ, NuHh = AlphaNDiff * ( DelXMin ** 2.0d0 ) / DelTimeLong NuVh = AlphaNDiff * ( DelZMin ** 2.0d0 ) / DelTimeLong NuHm = NuHh * NDiffRatio NuVm = NuVh * NDiffRatio if (myrank == 0) then call MessageNotify( "M", "dynamicalcore_init", "Nu = %f", d=(/Nu/) ) call MessageNotify( "M", "dynamicalcore_init", "NuHh = %f", d=(/NuHh/) ) call MessageNotify( "M", "dynamicalcore_init", "NuVh = %f", d=(/NuVh/) ) call MessageNotify( "M", "dynamicalcore_init", "NuHm = %f", d=(/NuHm/) ) call MessageNotify( "M", "dynamicalcore_init", "NuVm = %f", d=(/NuVm/) ) end if call HistoryAutoAddVariable( & & varname='PTempAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of potential temperature', & & units='K.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='PTempNDiff',& & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of potential temperature',& & units='K.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='CDensAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of cloud density', & & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='CDensNDiff',& & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of cloud density',& & units='K.m-3.s-1', & & xtype='double') do l = 1, nf 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='double') 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='double') end do call HistoryAutoAddVariable( & & varname='VelXAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of velocity (x)', & & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='VelXNDiff',& & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of velocity (x)',& & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='VelXPGrad', & & dims=(/'x','y','z','t'/), & & longname='Pressure gradient term of velocity (x)', & & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='VelYAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of velocity (y)', & & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='VelYNDiff',& & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of velocity (y)',& & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='VelYPGrad', & & dims=(/'x','y','z','t'/), & & longname='Pressure gradient term of velocity (y)', & & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='VelZAdv', & & dims=(/'x','y','z','t'/), & & longname='Advection term of velocity (z)', & & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='VelZNDiff',& & dims=(/'x','y','z','t'/), & & longname='Numerical diffusion term of Velocity (z)',& & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='VelZBuoyT',& & dims=(/'x','y','z','t'/), & & longname='Buoyancy (Temperature)',& & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='VelZBuoyM',& & dims=(/'x','y','z','t'/), & & longname='Buoyancy (MolWt)',& & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='VelZBuoyD',& & dims=(/'x','y','z','t'/), & & longname='Buoyancy (Drag)',& & units='K.m-3.s-1', & & xtype='double') call HistoryAutoAddVariable( & & varname='VelZPGrad', & & dims=(/'x','y','z','t'/), & & longname='Pressure gradient term of velocity (z)', & & units='K.m-3.s-1', & & xtype='double') end subroutine Dynamics_Init subroutine DynamicsLong( & & Time, DelTime, & ! (in) & pyz_VelXBl, pyz_VelXNl, & ! (in) & xqz_VelYBl, xqz_VelYNl, & ! (in) & xyr_VelZBl, xyr_VelZNl, & ! (in) & xyz_PTempBl, xyz_PTempNl, & ! (in) & xyz_ExnerBl, xyz_ExnerNl, & ! (in) & xyzf_QMixBl, xyzf_QMixNl, & ! (in) & xyz_KmBl, xyz_KmNl, & ! (in) & xyz_CDensBl, xyz_CDensNl, & ! (in) & pyz_DVelXDtNl, & ! (inout) & xqz_DVelYDtNl, & ! (inout) & xyr_DVelZDtNl, & ! (inout) & xyz_DPTempDtNl, & ! (inout) & xyzf_DQMixDtNl, & ! (inout) & xyz_DKmDtNl, & ! (inout) & xyz_DCDensDtNl & ! (inout) & ) implicit none real(DP), intent(in) :: Time, DelTime real(DP), intent(in) :: pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_PTempBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_PTempNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_ExnerBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_ExnerNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyzf_QMixBl(imin:imax,jmin:jmax,kmin:kmax, 1:nf) real(DP), intent(in) :: xyzf_QMixNl(imin:imax,jmin:jmax,kmin:kmax, 1:nf) real(DP), intent(in) :: xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(in) :: xyz_KmNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: pyz_DVelXDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: xqz_DVelYDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: xyz_DPTempDtNl(imin:imax,jmin:jmax,kmin:kmax) real(DP), intent(inout) :: xyzf_DQMixDtNl(imin:imax,jmin:jmax,kmin:kmax, 1:nf) real(DP), intent(inout) :: xyz_DKmDtNl(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) :: 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_NDiff(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_NDiff(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_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_Orig(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) real(DP) :: xyzf_Orig(imin:imax,jmin:jmax,kmin:kmax, 1:nf) real(DP) :: xyzf_Adv(imin:imax,jmin:jmax,kmin:kmax, 1:nf) real(DP) :: xyzf_NDiff(imin:imax,jmin:jmax,kmin:kmax, 1:nf) real(DP) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, 1:nf) real(DP) :: xyzf_QMixPerMolWt(imin:imax,jmin:jmax,kmin:kmax, 1:nf) integer :: f !------------------------------ ! tendency of cloud density ! ! initialize xyz_Orig = xyz_DCDensDtNl ! フラックス項の計算. 4 次精度中心差分を用いて計算 ! xyz_Adv = & & - xyz_dx_pyz(pyz_VelXNl * pyz_avr_xyz(xyz_CDensNl)) & & - xyz_dy_xqz(xqz_VelYNl * xqz_avr_xyz(xyz_CDensNl)) & & - xyz_dz_xyr(xyr_VelZNl * xyr_avr_xyz(xyz_CDensNl)) ! 数値粘性項の計算 xyz_NDiff = & & + NuHh * (xyz_dx_pyz(pyz_dx_xyz(xyz_CDensBl))) & & + NuHh * (xyz_dy_xqz(xqz_dy_xyz(xyz_CDensBl))) & & + NuVh * (xyz_dz_xyr(xyr_dz_xyz(xyz_CDensBl))) xyz_DCDensDtNl = xyz_Orig + xyz_Adv + xyz_NDiff call HistoryAutoPut(Time, 'CDensAdv', xyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(Time, 'CDensNDiff', xyz_NDiff(1:nx,1:ny,1:nz)) !---------------------------------- ! 拡散係数 ! ! Initialize ! xyz_Orig = xyz_DKmDtNl ! Advection term ! xyz_Adv = & & - xyz_avr_pyz(pyz_VelXNl * pyz_dx_xyz(xyz_KmNl)) & & - xyz_avr_xqz(xqz_VelYNl * xqz_dy_xyz(xyz_KmNl)) & & - xyz_avr_xyr(xyr_VelZNl * xyr_dz_xyz(xyz_KmNl)) ! Numerical diffusion term ! xyz_NDiff = & & NuHm * (xyz_dx_pyz(pyz_dx_xyz( xyz_KmBl ))) & & + NuHm * (xyz_dy_xqz(xqz_dy_xyz( xyz_KmBl ))) & & + NuVm * (xyz_dz_xyr(xyr_dz_xyz( xyz_KmBl ))) xyz_DKmDtNl = xyz_Orig + xyz_Adv + xyz_NDiff !------------------------------------- ! 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_dx_xyz(xyz_PTempAll)) & & - xyz_avr_xqz(xqz_VelYNl * xqz_dy_xyz(xyz_PTempAll)) & & - xyz_avr_xyr(xyr_VelZNl * xyr_dz_xyz(xyz_PTempAll)) ! numerical diffusion term ! xyz_Num = xyz_NumDiffScalar( xyz_PTempBl) ! xyz_NDiff = & & NuHh * (xyz_dx_pyz(pyz_dx_xyz( xyz_PTempBl ))) & & + NuHh * (xyz_dy_xqz(xqz_dy_xyz( xyz_PTempBl ))) & & + NuVh * (xyz_dz_xyr(xyr_dz_xyz( xyz_PTempBl ))) ! sum ! xyz_DPTempDtNl = xyz_Orig + xyz_Adv + xyz_NDiff ! output ! call HistoryAutoPut(Time, 'PTempAdv', xyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(Time, 'PTempNDiff', xyz_NDiff(1:nx,1:ny,1:nz)) !------------------------------ ! tendency of mixing ratio ! ! initialize xyzf_Orig = xyzf_DQMixDtNl xyzf_QMixAll = xyzf_QMixNl + xyzf_QMixBZ do f = 1, nf ! 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_dx_xyz(xyzf_QMixAll(:,:,:,f))) & & - xyz_avr_xqz(xqz_VelYNl * xqz_dy_xyz(xyzf_QMixAll(:,:,:,f))) & & - xyz_avr_xyr(xyr_VelZNl * xyr_dz_xyz(xyzf_QMixAll(:,:,:,f))) ! numerical diffusion term ! xyzf_Diff = xyzf_NumDiffScalar(xyzf_QMixBl) ! xyzf_NDiff(:,:,:,f) = & & NuHh * (xyz_dx_pyz(pyz_dx_xyz( xyzf_QMixBl(:,:,:,f) ))) & & + NuHh * (xyz_dy_xqz(xqz_dy_xyz( xyzf_QMixBl(:,:,:,f) ))) & & + NuVh * (xyz_dz_xyr(xyr_dz_xyz( xyzf_QMixBl(:,:,:,f) ))) end do ! sum ! xyzf_DQMixDtNl = xyzf_Orig + xyzf_Adv + xyzf_NDiff ! output ! do f = 1, nf call HistoryAutoPut(Time, trim(SpcWetSymbol(f))//'_Adv', xyzf_Adv(1:nx,1:ny,1:nz,f)) call HistoryAutoPut(Time, trim(SpcWetSymbol(f))//'_NDiff', xyzf_NDiff(1:nx,1:ny,1:nz,f)) end do !------------------------------ ! 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_dx_pyz( pyz_VelXNl ) ) & & - pyz_avr_pqz( pqz_avr_xqz( xqz_VelYNl ) * pqz_dy_pyz( pyz_VelXNl ) ) & & - pyz_avr_pyr( pyr_avr_xyr( xyr_VelZNl ) * pyr_dz_pyz( pyz_VelXNl ) ) ! Numerical diffusion term !pyz_Diff = pyz_NumDiffVelX(pyz_VelXBl) pyz_NDiff = & & NuHm * ( pyz_dx_xyz( xyz_dx_pyz( pyz_VelXBl ) ) ) & & + NuHm * ( pyz_dy_pqz( pqz_dy_pyz( pyz_VelXBl ) ) ) & & + NuVm * ( pyz_dz_pyr( pyr_dz_pyz( pyz_VelXBl ) ) ) ! sum ! pyz_DVelXDtNl = pyz_Orig + pyz_Adv + pyz_NDiff call HistoryAutoPut(Time, 'VelXAdv', pyz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(Time, 'VelXNDiff', pyz_NDiff(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_dx_xqz( xqz_VelYNl ) ) & & - xqz_VelYNl * xqz_avr_xyz( xyz_dy_xqz( xqz_VelYNl ) ) & & - xqz_avr_xqr( xqr_avr_xyr( xyr_VelZNl ) * xqr_dz_xqz( xqz_VelYNl ) ) ! Numerical diffusion term ! xqz_Diff = xqz_NumDiffVelY(xqz_VelYBl) ! xqz_NDiff = & & NuHm * ( xqz_dx_pqz( pqz_dx_xqz( xqz_VelYBl ) ) ) & & + NuHm * ( xqz_dy_xyz( xyz_dy_xqz( xqz_VelYBl ) ) ) & & + NuVm * ( xqz_dz_xqr( xqr_dz_xqz( xqz_VelYBl ) ) ) ! sum ! xqz_DVelYDtNl = xqz_Orig + xqz_Adv + xqz_NDiff call HistoryAutoPut(Time, 'VelYAdv', xqz_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(Time, 'VelYNDiff', xqz_NDiff(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_dx_xyr( xyr_VelZNl ) ) & & - xyr_avr_xqr( xqr_avr_xqz( xqz_VelYNl ) * xqr_dy_xyr( xyr_VelZNl ) ) & & - xyr_VelZNl * xyr_avr_xyz( xyz_dz_xyr( xyr_VelZNl ) ) ! Numerical diffusion term !xyr_Diff = xyr_NumDiffVelZ(xyr_VelZBl) ! xyr_NDiff = & & NuHm * ( xyr_dx_pyr( pyr_dx_xyr( xyr_VelZBl ) ) ) & & + NuHm * ( xyr_dy_xqr( xqr_dy_xyr( xyr_VelZBl ) ) ) & & + NuVm * ( 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:nf), 4) ) & & / ( 1.0d0 + xyr_QMixBZ ) ! sum ! ! xyr_DVelZDtNl = xyr_Orig + xyr_Adv + xyr_NDiff + xyr_BuoyT + xyr_BuoyM + xyr_BuoyD xyr_DVelZDtNl = xyr_Orig + xyr_Adv + xyr_NDiff + xyr_BuoyT call HistoryAutoPut(Time, 'VelZAdv', xyr_Adv(1:nx,1:ny,1:nz)) call HistoryAutoPut(Time, 'VelZNDiff', xyr_NDiff(1:nx,1:ny,1:nz)) call HistoryAutoPut(Time, 'VelZBuoyT', xyr_BuoyT(1:nx,1:ny,1:nz)) call HistoryAutoPut(Time, 'VelZBuoyM', xyr_BuoyM(1:nx,1:ny,1:nz)) call HistoryAutoPut(Time, 'VelZBuoyD', xyr_BuoyD(1:nx,1:ny,1:nz)) end subroutine DynamicsLong subroutine DynamicsShortHEVI( & & Time, DelTime, & ! (in) & pyz_VelXNs, & ! (in) & xqz_VelYNs, & ! (in) & xyr_VelZNs, & ! (in) & xyz_ExnerNs, & ! (in) & pyz_DVelXDtNl, & ! (in) & xqz_DVelYDtNl, & ! (in) & xyr_DVelZDtNl, & ! (in) & xyz_DExnerDtNs, & ! (in) & pyz_VelXAs, & ! (out) & xqz_VelYAs, & ! (out) & xyr_VelZAs, & ! (out) & xyz_ExnerAs & ! (out) & ) real(DP), intent(in) :: Time, DelTime 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(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_DivVel(imin:imax,jmin:jmax,kmin:kmax) ! initialize: Divergence of velocity ! xyz_DivVel = & & xyz_dx_pyz( pyz_VelXNs ) & & + xyz_dy_xqz( xqz_VelYNs ) & & + xyz_dz_xyr( xyr_VelZNs ) !-------------------------------------- ! VelX ! pyz_DVelXDtNs = & & - pyz_avr_xyz( CpDry * xyz_PTempBZ ) & & * ( pyz_dx_xyz( xyz_ExnerNs ) - pyz_dx_xyz( Nu * xyz_DivVel ) ) ! Time integration ! pyz_VelXAs = pyz_VelXNs + DelTime * (pyz_DVelXDtNl + pyz_DVelXDtNs) call HistoryAutoPut(Time, 'VelXPGrad', pyz_DVelXDtNs(1:nx,1:ny,1:nz)) !-------------------------------------- ! VelY ! xqz_DVelYDtNs = & & - xqz_avr_xyz( CpDry * xyz_PTempBZ ) & & * ( xqz_dy_xyz( xyz_ExnerNs ) - xqz_dy_xyz( Nu * xyz_DivVel ) ) ! Time integration ! xqz_VelYAs = xqz_VelYNs + DelTime * (xqz_DVelYDtNl + xqz_DVelYDtNs) call HistoryAutoPut(Time, 'VelYPGrad', xqz_DVelYDtNs(1:nx,1:ny,1:nz)) !-------------------------------------- ! Exner function ! xyz_ExnerAs = xyz_Exner( & & pyz_VelXNs, & & pyz_VelXAs, & & xqz_VelYNs, & & xqz_VelYAs, & & xyr_VelZNs, & & xyz_ExnerNs, & & xyr_DVelZDtNl, & & xyz_DExnerDtNs & & ) !-------------------------------------- ! VelZ ! xyr_DVelZDtNs = & & - xyr_avr_xyz(CpDry * xyz_PTempBZ ) & & * ( & & beta * xyr_dz_xyz( xyz_ExnerAs ) & & + (1.0d0 - beta) * xyr_dz_xyz( xyz_ExnerNs ) & & - xyr_dz_xyz( Nu * xyz_DivVel ) & & ) ! Time integration ! xyr_VelZAs = xyr_VelZNs + DelTime * (xyr_DVelZDtNl + xyr_DVelZDtNs) call HistoryAutoPut(Time, 'VelZPGrad', xyr_DVelZDtNs(1:nx,1:ny,1:nz)) end subroutine DynamicsShortHEVI !!!--------------------------------------------------------------------!!! subroutine DynamicsVI_init() ! !エクスナー関数を陰解法で解く際に必要となる, 係数行列の要素を決め, !LU 分解を行う. ! !暗黙の型宣言禁止 implicit none real(DP) :: DTS ! 短い時間格子 DTS = DelTimeShort !配列の割り付け allocate( A(1:nz) ) allocate( B(2: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 方向に一様なので. ! RegXMax, RegYMax の値で代表させることとした. 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(2:nz-1) = & & (beta ** 2.0d0) & & * xyz_F1BZ(1,1,2:nz-1) & & * ( & & xyr_F2BZ(1,1,2:nz-1) / r_dz(2:nz-1) & & + xyr_F2BZ(1,1,1:nz-2) / r_dz(1:nz-2) & & ) & & * (DTS ** 2.0d0) & & / z_dz(2: nz-1) & & + 1.0d0 A(1) = & & (beta ** 2.0d0) & & * xyz_F1BZ(1,1,1) & & * xyr_F2BZ(1,1,1) / r_dz(1) & & * (DTS ** 2.0d0) & & / z_dz(1) & & + 1.0d0 A(nz) = & & (beta ** 2.0d0) & & * xyz_F1BZ(1,1,nz) & & * xyr_F2BZ(1,1,nz-1) / r_dz(nz-1) & & * (DTS ** 2.0d0) & & / z_dz(nz) & & + 1.0d0 B(2:nz) = & & - (beta ** 2.0d0) & & * xyz_F1BZ(1,1,1:nz-1) & & * xyr_F2BZ(1,1,1:nz-1) / r_dz(1:nz-1) & & * (DTS ** 2.0d0) & & / z_dz(1:nz-1) C(1: nz-1) = & & - ( beta ** 2.0d0 ) & & * xyz_F1BZ(1,1,2:nz) & & * xyr_F2BZ(1,1,1:nz-1) / r_dz(1:nz-1) & & * (DTS ** 2.0d0) & & / z_dz(1:nz) !デバッグ出力ファイルオープン. 情報取得. ! call FileOpen(unit, file="dynimpfunc_log", mode='w') ! write(unit,*) "kz, A_k" ! write(unit,*) "------------------------------------" ! do kz = RegZmin, RegZmax ! write(unit,*) kz, A(kz) ! end do ! write(unit,*) "------------------------------------" ! write(unit,*) "kz, B_k" ! write(unit,*) "------------------------------------" ! do kz = RegZmin+1, RegZmax ! write(unit,*) kz, B(kz) ! end do ! write(unit,*) "------------------------------------" ! write(unit,*) "kz, C_k" ! write(unit,*) "------------------------------------" ! do kz = RegZmin, RegZmax-1 ! write(unit,*) kz, C(kz) ! end do ! write(unit,*) "------------------------------------" ! close(unit) !---------------------------------------------------------------- ! 係数行列を 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) ) !LU 分解の実行 ! LAPACK の利用 call ResolvLU_Lapack( ) end subroutine DynamicsVI_init !!!--------------------------------------------------------------------!!! function xyz_Exner( & & pyz_VelXNs, & & pyz_VelXAs, & & xqz_VelYNs, & & xqz_VelYAs, & & xyr_VelZNs, & & xyz_ExnerNs, & & xyr_DVelZDtNl, & & xyz_DExnerDtNs & & ) ! !陰解法を用いたエクスナー関数の計算. ! !暗黙の型宣言禁止 implicit none !入出力変数 real(DP), intent(in) :: pyz_VelXNs & & (imin:imax,jmin:jmax,kmin:kmax) !速度 u [τ] real(DP), intent(in) :: pyz_VelXAs & & (imin:imax,jmin:jmax,kmin:kmax) !速度 u [τ+Δτ] real(DP), intent(in) :: xqz_VelYNs & & (imin:imax,jmin:jmax,kmin:kmax) !速度 v [τ] 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) :: xyr_DVelZDtNl & & (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(1:nx*ny,1: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 DTS = DelTimeShort !変数の初期化 xyz_Exner = 0.0d0 !速度の収束を計算 xyz_DivVelNs = xyz_dx_pyz( pyz_VelXNs ) & & + xyz_dy_xqz( xqz_VelYNs ) & & + xyz_dz_xyr( xyr_VelZNs ) !行列計算のための係数 E = xyr_dz_xyz( Nu * xyz_DivVelNs ) & & - ( 1.0d0 - beta ) * xyr_dz_xyz( xyz_ExnerNs ) & & + 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 ) & & - xyr_dz_xyz( Nu * xyz_DivVelNs ) & & ) & & + xyr_DVelZDtNl * DTS & & ) & & ) & & + xyz_DExnerDtNs * 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(:,:,0) * E(:,:,0) / 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 DynamicalCore