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(:) !部分ピボット交換の情報を格納 integer, save :: RegXMin = 1 integer, save :: RegXMax = 1 integer, save :: RegYMin = 1 integer, save :: RegYMax = 1 integer, save :: RegZMin = 1 integer, save :: RegZMax = 1 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(inout), optional :: xyz_PTempAl(imin:imax,jmin:jmax,kmin:kmax) ! real(DP), intent(inout), optional :: xyzf_QMixAl(imin:imax,jmin:jmax,kmin:kmax, 1:nf) 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)) ! if (present(xyz_PTempAl)) then ! xyz_PTempAl = xyz_PTempBl + DelTime * xyz_DPTempDtNl ! end if !------------------------------ ! 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 ! time integration ! ! if (present(xyz_PTempAl)) then ! xyzf_QMixAl = xyzf_QMixBl + DelTime * xyzf_DQMixDtNl ! end if !------------------------------ ! 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 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) xqz_VelYAs = xqz_VelYNs 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 RegXMin = 1 RegXMax = nx RegYMin = 1 RegYMax = ny RegZMin = 1 RegZMax = nz !配列の割り付け allocate( A(RegZMin:RegZMax) ) allocate( B(RegZMin+1:RegZMax) ) allocate( C(RegZMin:RegZMax-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(RegZMin+1: RegZMax-1) = & & (beta ** 2.0d0) & & * xyz_F1BZ(RegXMax,RegYMax,RegZMin+1: RegZMax-1) & & * ( & & xyr_F2BZ(RegXMax,RegYMax,RegZMin+1: RegZMax-1) & & / r_dz(RegZMin+1: RegZMax-1) & & + xyr_F2BZ(RegXMax,RegYMax,RegZMin : RegZMax-2) & & / r_dz(RegZMin: RegZMax-2) & & ) & & * (DTS ** 2.0d0) & & / z_dz(RegZMin+1: RegZMax-1) & & + 1.0d0 A(RegZMin) = & & (beta ** 2.0d0) & & * xyz_F1BZ(RegXMax,RegYMax,RegZMin) & & * xyr_F2BZ(RegXMax,RegYMax,RegZMin) & & / r_dz(RegZMin) & & * (DTS ** 2.0d0) & & / z_dz(RegZMin) & & + 1.0d0 A(RegZMax) = & & (beta ** 2.0d0) & & * xyz_F1BZ(RegXMax,RegYMax,RegZMax) & & * xyr_F2BZ(RegXMax,RegYMax,RegZMax-1) & & / r_dz(RegZMax-1) & & * (DTS ** 2.0d0) & & / z_dz(RegZMax) & & + 1.0d0 B(RegZMin+1:RegZMax) = & & - (beta ** 2.0d0) & & * xyz_F1BZ(RegXMax,RegYMax,RegZMin:RegZMax-1) & & * xyr_F2BZ(RegXMax,RegYMax,RegZMin:RegZMax-1) & & / r_dz(RegZmin:RegZMax-1) & & * (DTS ** 2.0d0) & & / z_dz(RegZMin:RegZMax-1) C(RegZMin: RegZMax-1) = & & - ( beta ** 2.0d0 ) & & * xyz_F1BZ(RegXMax,RegYMax,RegZMin+1:RegZMax) & & * xyr_F2BZ(RegXMax,RegYMax,RegZMin :RegZMax-1) & & / r_dz(RegZmin:RegZMax-1) & & * (DTS ** 2.0d0) & & / z_dz(RegZMin+1:RegZMax) !デバッグ出力ファイルオープン. 情報取得. ! 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 = RegZMax - RegZMin +1 !係数行列/改行列の次数, 整合寸法 M = (RegXMax - RegXMin +1)*(RegYMax - RegYMin +1) !方程式の組数 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(RegXMin:RegXMax,RegYMin:RegYMax,RegZMin:RegZMax) real(DP) :: D((RegXMax-RegXMin+1)*(RegYMax-RegYMin+1),(RegZMax-RegZMin+1)) 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(:,:,RegZMin) = D1(:,:,RegZMin) & & - beta * xyz_F1BZ(:,:,RegZMin) * (DTS ** 2.0d0) & & * xyr_F2BZ(:,:, RegZMin-1) * E(:,:,RegZMin-1) & & / z_dz(RegZMin) D1(:,:,RegZMax) = D1(:,:,RegZMax) & & + beta * xyz_F1BZ(:,:,RegZMax) * (DTS ** 2.0d0) & & * xyr_F2BZ(:,:,RegZMax) * E(:,:,RegZMax) & & / z_dz(RegZMax) D2 = D1(RegXMin:RegXMax,RegYMin:RegYMax,RegZMin:RegZMax) do kz = RegZMin, RegZMax do jy = RegYMin, RegYMax do ix = RegXMin, RegXMax D(ix + (RegXMax-RegXMin+1) * (jy - 1), kz) = D2(ix,jy,kz) end do end do end do !----------------------------------------------------------- !連立一次方程式の解を求める !------------------------------------------------------------ !解の計算 ! LAPACK 利用 call LinSolv_Lapack( D ) !戻り値を出力 do kz = RegZMin, RegZMax do jy = RegYMin, RegYMax do ix = RegXMin, RegXMax xyz_Exner(ix,jy,kz) = D(ix + (RegXMax-RegXMin+1) * (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