!c Description: !c 非静力学モデル deepconv/arare !c !c Current Code Owner: !c sugiyama@gfd-dennou.org !c !c Copyright (C) SUGIYAMA Ko-ichiro, 2004, All rights reserved program arare use gt4_history use gridset, only: gridset_init, DimXMin, DimXMax, DimZMin, DimZMax, & & RegZMin, RegZMax use fileset, only: fileset_init use timeset, only: timeset_init, NstepTShort, DelTShort use physprm, only: physprm_init use expname, only: expname_init use basicprm, only: basicprm_init use linlib, only: linlib_init use if_display implicit none real(8), allocatable :: ss_Exner_A(:,:) real(8), allocatable :: fs_VelX_A(:,:) real(8), allocatable :: sf_VelZ_A(:,:) real(8), allocatable :: ss_Exner_B(:,:) real(8), allocatable :: fs_VelX_B(:,:) real(8), allocatable :: sf_VelZ_B(:,:) real(8), allocatable :: fs_Fx_B(:,:) real(8), allocatable :: sf_Fz_B(:,:) real(8), allocatable :: ss_VelDiv_B(:,:) real(8) :: Time integer :: t ! integer, parameter :: NStepMax = 50 ! integer :: NStepInterval type(gt_history) :: ss_reg, fs_reg, sf_reg type(gt_history) :: ss_dim, fs_dim, sf_dim !!! !!!初期化 !!! !--- I/O ファイル名の初期化 call fileset_init !--- 格子点情報の初期化 call gridset_init !--- 時間積分の設定の初期化 call timeset_init !--- 物理量の初期化 call physprm_init !--- 実験名の初期化 call expname_init !--- 基本場の値を設定 call basicprm_init !--- 行列計算ライブラリの初期化. !--- 引数はスカラーの z 方向の計算領域の大きさ call linlib_init(RegZMax - RegZMin) !--- ファイル出力を行う時間間隔を決定 ! if (NStepTShort > NStepMax) then ! NStepInterval = NStepTShort / NStepMax ! else ! NStepInterval = 1 ! end if !--- 出力ファイルの定義 call OpenRegFile call OpenDimFile !!! !!! CFL 条件のチェック !!! call cflchk !!! !!!配列の確保, 初期化 !!! call ArareAlloc call ArareArrayInit !--- 初期擾乱の設定 call ArareDisturb !!! !!!計算 !!! !--- とりあえず長い時間ステップは無視 fs_Fx_B = 0.0d0 sf_Fz_B = 0.0d0 !--- 最初の時間 Time = 0.0d0 call OutputRegFile call OutputDimFile !--- 短いタイムステップでの積分 do t = 1, NstepTShort Time = Time + DelTShort write(*,*) "Time: ", Time !--- 速度の収束 call VelDiv(fs_VelX_B, sf_VelZ_B, ss_VelDiv_B) !--- 速度 u の計算 call VelXHeVi(ss_Exner_B, ss_VelDiv_B, fs_Fx_B, fs_VelX_B, fs_VelX_A) !--- エクスナー関数の計算 call ExnerHeVi( & & ss_VelDiv_B, sf_Fz_B, fs_VelX_A, sf_VelZ_B, ss_Exner_B, ss_Exner_A) !--- 速度 w の計算 call VelZHeVi( & & ss_Exner_B, ss_Exner_A, ss_VelDiv_B, sf_Fz_B, sf_VelZ_B, sf_VelZ_A) !--- ループを回すための処置 ss_Exner_B = ss_Exner_A fs_VelX_B = fs_VelX_A sf_VelZ_B = sf_VelZ_A !--- ファイル出力 ! if (mod(t, NStepInterval) == 0) then call OutputRegFile call OutputDimFile ! end if end do !--- 出力ファイルのクローズ call CloseRegFile call CloseDimFile contains subroutine ArareAlloc allocate(ss_Exner_A(DimXMin:DimXMax, DimZMin:DimZMax), & & fs_VelX_A(DimXMin:DimXMax, DimZMin:DimZMax) , & & sf_VelZ_A(DimXMin:DimXMax, DimZMin:DimZMax), & & ss_Exner_B(DimXMin:DimXMax, DimZMin:DimZMax), & & fs_VelX_B(DimXMin:DimXMax, DimZMin:DimZMax), & & sf_VelZ_B(DimXMin:DimXMax, DimZMin:DimZMax), & & fs_Fx_B(DimXMin:DimXMax, DimZMin:DimZMax), & & sf_Fz_B(DimXMin:DimXMax, DimZMin:DimZMax), & & ss_VelDiv_B(DimXMin:DimXMax, DimZMin:DimZMax) ) end subroutine ArareAlloc subroutine ArareArrayInit use basicprm, only: ss_ExnerBasicZ ss_Exner_A = ss_ExnerBasicZ ss_Exner_B = ss_Exner_A fs_VelX_A = 0.0d0 fs_VelX_B = 0.0d0 sf_VelZ_A = 0.0d0 sf_VelZ_B = 0.0d0 fs_Fx_B = 0.0d0 sf_Fz_B = 0.0d0 ss_VelDiv_B = 0.0d0 Time = 0.0d0 end subroutine ArareArrayInit subroutine ArareDisturb use gridset, only: s_X, s_Z, Xmin, Xmax, Zmin, Zmax, & & DimXMin, DimXMax, DimZMin, DimZMax use physprm, only: pi integer :: i, k real(8) :: mu real(8) :: sigma !--- エクスナー関数の擾乱 mu = ( Zmax - Zmin ) * 1.0d0 / 5.0d0 sigma = ( Zmax - Zmin ) / 20.0d0 do k = DimZMin, DimZMax ss_Exner_A(:,k) = ss_Exner_A(:,k) & & + 10.0d0 * ss_Exner_A(:,k) & & * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) & & / (dsqrt(2.0d0 * pi) * sigma) end do ! mu = ( Xmax - Xmin ) * 4.0d0 / 5.0d0 ! sigma = ( Xmax - Xmin ) / 20.0d0 ! ! do i = DimXMin, DimXMax ! ss_Exner_A(i,:) = ss_Exner_A(i,:) & ! & + 100.0d0 * ss_Exner_A(i,:) & ! & * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) & ! & / (dsqrt(2.0d0 * pi) * sigma) ! end do ss_Exner_B = ss_Exner_A end subroutine ArareDisturb !--- gtool4 出力関連 subroutine OpenRegFile use gridset, only: RegXMin, RegXMax, RegZMin, RegZMax, & & NX, NZ, f_X, f_Z, s_X, s_Z use expname, only: exptitle, expsrc, expinst use fileset, only: fs_RegFile, sf_RegFile, ss_RegFile ! ヒストリー作成 call HistoryCreate( & & file = fs_RegFile, & & title = exptitle, & & source = expsrc, & & institution = expinst, & & dims=(/'x','z','t'/), & & dimsizes=(/NX+1, NZ, 0/), & & longnames=(/'X-coordinate', & & 'Z-coordinate', & & 'Time '/), & & units=(/'m','m','s'/), origin=0.0, & & interval=0.0, & & history = fs_reg) ! 変数出力 call HistoryPut('x', f_X(RegXMin: RegXMax), fs_reg) call HistoryPut('z', s_Z(RegXMin+1: RegXMax), fs_reg) ! ヒストリー作成 call HistoryCreate( & & file = sf_RegFile, & & title = exptitle, & & source = expsrc, & & institution = expinst, & & dims=(/'x','z','t'/), & & dimsizes=(/NX, NZ+1, 0/), & & longnames=(/'X-coordinate', & & 'Z-coordinate', & & 'Time '/), & & units=(/'m','m','s'/), & & origin=0.0, & & interval=0.0, & & history = sf_reg) ! 変数出力 call HistoryPut('x', s_X(RegXMin+1: RegXMax), sf_reg) call HistoryPut('z', f_Z(RegXMin: RegXMax), sf_reg) ! ヒストリー作成 call HistoryCreate( & & file = ss_RegFile, & & title = exptitle, & & source = expsrc, & & institution = expinst, & & dims=(/'x','z','t'/), & & dimsizes=(/NX, NZ, 0/), & & longnames=(/'X-coordinate', & & 'Z-coordinate', & & 'Time '/), & & units=(/'m','m','s'/), origin=0.0, & & interval=0.0, & & history = ss_reg) ! 変数出力 call HistoryPut('x', s_X(RegXMin+1: RegXMax), ss_reg) call HistoryPut('z', s_Z(RegXMin+1: RegXMax), ss_reg) ! 無次元圧力 call HistoryAddVariable( & & varname='Exner', dims=(/'x','z','t'/), & & longname='nondimensional pressure', units='1',& & xtype='double', & & history = ss_reg) ! 速度 call HistoryAddVariable( & & varname='VelX', dims=(/'x','z','t'/), & & longname='zonal velocity', & & units='m/s', xtype='double', & & history = fs_reg) ! 速度 call HistoryAddVariable( & & varname='VelZ', dims=(/'x','z','t'/), & & longname='vertical velocity', & & units='m/s', xtype='double', & & history = sf_reg) end subroutine OpenRegFile subroutine OutputRegFile use gridset, only: RegXMin, RegXMax, RegZMin, RegZMax call HistoryPut('t', Time, ss_reg) call HistoryPut('t', Time, fs_reg) call HistoryPut('t', Time, sf_reg) call HistoryPut('Exner', & & ss_Exner_A(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg) call HistoryPut('VelX', & & fs_VelX_A(RegXMin:RegXMax, RegZMin+1:RegZMax), fs_reg) call HistoryPut('VelZ', & & sf_VelZ_A(RegXMin+1:RegXMax, RegZMin:RegZMax), sf_reg) end subroutine OutputRegFile subroutine CloseRegFile call HistoryClose( ss_reg ) call HistoryClose( fs_reg ) call HistoryClose( sf_reg ) end subroutine CloseRegFile !--- gtool4 出力関連 subroutine OpenDimFile use gridset, only: f_X, f_Z, s_X, s_Z use expname, only: exptitle, expsrc, expinst use fileset, only: fs_DimFile, sf_DimFile, ss_DimFile integer:: N N = size(f_X, 1) ! ヒストリー作成 call HistoryCreate( & & file = fs_DimFile, & & title = exptitle, & & source = expsrc, & & institution = expinst, & & dims=(/'x','z','t'/), & & dimsizes=(/N, N, 0/), & & longnames=(/'X-coordinate', & & 'Z-coordinate', & & 'Time '/), & & units=(/'m','m','s'/), origin=0.0, & & interval=0.0, & & history = fs_dim) ! 変数出力 call HistoryPut('x', f_X, fs_dim) call HistoryPut('z', s_Z, fs_dim) ! ヒストリー作成 call HistoryCreate( & & file = sf_DimFile, & & title = exptitle, & & source = expsrc, & & institution = expinst, & & dims=(/'x','z','t'/), & & dimsizes=(/N, N, 0/), & & longnames=(/'X-coordinate', & & 'Z-coordinate', & & 'Time '/), & & units=(/'m','m','s'/), & & origin=0.0, & & interval=0.0, & & history = sf_dim) ! 変数出力 call HistoryPut('x', s_X, sf_dim) call HistoryPut('z', f_Z, sf_dim) ! ヒストリー作成 call HistoryCreate( & & file = ss_DimFile, & & title = exptitle, & & source = expsrc, & & institution = expinst, & & dims=(/'x','z','t'/), & & dimsizes=(/N, N, 0/), & & longnames=(/'X-coordinate', & & 'Z-coordinate', & & 'Time '/), & & units=(/'m','m','s'/), origin=0.0, & & interval=0.0, & & history = ss_dim) ! 変数出力 call HistoryPut('x', s_X, ss_dim) call HistoryPut('z', s_Z, ss_dim) ! 無次元圧力 call HistoryAddVariable( & & varname='Exner', dims=(/'x','z','t'/), & & longname='nondimensional pressure', units='1',& & xtype='double', & & history = ss_dim) ! 速度 call HistoryAddVariable( & & varname='VelX', dims=(/'x','z','t'/), & & longname='zonal velocity', & & units='m/s', xtype='double', & & history = fs_dim) ! 速度 call HistoryAddVariable( & & varname='VelZ', dims=(/'x','z','t'/), & & longname='vertical velocity', & & units='m/s', xtype='double', & & history = sf_dim) end subroutine OpenDimFile subroutine OutputDimFile call HistoryPut('t', Time, ss_dim) call HistoryPut('t', Time, fs_dim) call HistoryPut('t', Time, sf_dim) call HistoryPut('Exner', ss_Exner_A, ss_dim) call HistoryPut('VelX', fs_VelX_A, fs_dim) call HistoryPut('VelZ', sf_VelZ_A, sf_dim) end subroutine OutputDimFile subroutine CloseDimFile call HistoryClose( ss_dim ) call HistoryClose( fs_dim ) call HistoryClose( sf_dim ) end subroutine CloseDimFile end program arare