program sokuryo15 implicit none c ========== 変数の定義ここから ========== c -- 空間の分割点数(id:x方向、jd:y方向) integer id,jd parameter (id = 51) parameter (jd = 31) c -- 今の瞬間の濃度の値を入れる配列 real dens_old(id,jd) c -- 次の瞬間の濃度の値を入れる配列 real dens_new(id,jd) c -- x方向の流速値 real uvel(id,jd) c -- y方向の流速値 real vvel(id,jd) c -- 時間と空間を分割するサイズ c (dt:時間微分を近似するための微小時間) c (dx:空間微分を近似するための微小距離) real dx,dt c -- 時刻方向に計算した回数(dt x itime = 計算開始からの時間) integer itime c -- 大気が濃度を拡散させる速さのパラメータ(拡散係数) real diff c -- 拡散物質を放出する点の番号 integer isrc integer jsrc c -- 出力ファイル名のための変数、文字列 integer num character(4) cnum character(40) cfile c -- integer i,j real x,y c ========== 変数の定義ここまで ========== c*********************************************************************** c 1.シミュレーションする空間、時間を分割する間隔等、計算のための条件設定 c*********************************************************************** c 空間の分割間隔dx dx = 2.0e0 c 時間の分割間隔 dt = 0.1e0 c 大気中を物質が分散する速度を与えるパラメータ(拡散係数[m^2/s]) diff = 0.5e0 c -- (追加課題)拡散物質を放出する点の番号を変更した計算も行いましょう c (isrcとjsrcの値を設定) isrc = 5 jsrc = 3 c*********************************************************************** c 2. 最初の瞬間の濃度と流速のデータを与える(データの読み込み) c*********************************************************************** c (課題)速度のデータとして、データファイル"wind.dat"から c x方向速度を配列uvelに、y方向速度を配列vvelに c 読み込む処理を加えなさい。 open(20,file='wind.dat') do j=1,jd do i=1,id read(20,*) uvel(i,j), vvel(i,j) enddo enddo close(20) c -- 濃度の初期状態として、全ての点で0を設定する do i=1,id do j=1,jd dens_new(i,j) = 0.0e0 enddo enddo c*********************************************************************** c 3.−5.の繰り返しループ c*********************************************************************** c -- 次の時刻の状態予測を繰り返すループ do itime=1,2000 c ---- 時間を一段ずらす c (5.次の瞬間が現在となり、さらに次の瞬間の計算をするため3.に戻る) c 前の瞬間に「次の瞬間の濃度(dens_new)」だったもの c ⇒ 「次の瞬間」になったので、「現在の濃度(dens_old)」にする c (課題)ここに、全ての点でこの処理を行うプログラムを加えなさい do i=1,id do j=1,jd dens_old(i,j) = dens_new(i,j) enddo enddo c*********************************************************************** c 3.空間を分割した各点において、与えられた流速と現在の濃度を元に、 c 次の瞬間の各点での濃度を計算する c*********************************************************************** c ---- (課題)ここで、次の瞬間の濃度を計算するためのsubroutine c "iryu_kaku"をcallしなさい。 call iryu_kaku(id,jd,dx,dt,diff,uvel,vvel,dens_old,dens_new) c*********************************************************************** c 3.空間を分割した各点において、与えられた流速と現在の濃度を元に、 c 次の瞬間の各点での濃度を計算する c*********************************************************************** c ---- (課題)ここで、境界条件を処理するためのsubroutine c "kyoukai"をcallしなさい call kyoukai(id,jd,isrc,jsrc,uvel,vvel,dens_old,dens_new) c*********************************************************************** c 6.適当な時間おきに計算で得られた瞬間の濃度を書き出してゆく c*********************************************************************** c ---- ループ回数が100の倍数(100で割って余りが0)の時、出力を実行する if (mod(itime,100).eq.0) then c ---- ファイルを上書きしないよう、出力ファイル名として c ループ回数の数字が入ったファイル名を与える c 以下の処理で、cfileという変数には c "dens"+ループ回数(4桁整数)+".dat" c という、文字列が代入される。 write(cnum,'(i4.4)') itime cfile = 'dens' // cnum // '.dat' c ---- (課題)文字列"cfile"を出力ファイル名として c 全ての点での配列dens_newの値を出力する処理を加えなさい c このとき、gnuplotで2次元分布の表示を行う場合を考えて、 c jd行 × id列 のマトリクスになるように出力すると良い。 open (20,file=cfile) do j=1,jd write(20,100) (dens_new(i,j),i=1,id) enddo close(20) 100 format(51f10.6) endif enddo c*********************************************************************** c 3.−5.の繰り返しループ ここまで c*********************************************************************** stop end c c*********************************************************************** c 境界条件を与えるsubroutine c*********************************************************************** subroutine kyoukai(id,jd,isrc,jsrc,uvel,vvel,dens_old,dens_new) integer id,jd integer isrc,jsrc real dens_new(id,jd) real dens_old(id,jd) real uvel(id,jd) real vvel(id,jd) integer i,j c -- x方向の端の処理(境界条件) do j=2,jd-1 dens_new( 1,j) = dens_new( 2,j) dens_new(id,j) = dens_new(id-1,j) enddo c -- y方向の端の処理(境界条件) do i=2,id-1 dens_new(i, 1) = dens_new(i, 2) dens_new(i,jd) = dens_new(i,jd-1) enddo c -- 建物の処理(境界条件) c 1.「その場所の速度が0」、かつ、「隣の速度が0ではない」場合 c ⇒ 建物の壁 ⇒ 領域の端と同じ条件 c 2.「その場所の速度が0」、かつ、「隣の速度も0」の場合 c ⇒ 建物の中 ⇒ 濃度0 do i=2,id-1 do j=2,jd-1 if(uvel(i,j).eq.0.0.and.vvel(i,j).eq.0.0) then if(uvel(i+1,j).ne.0.0.or.vvel(i+1,j).ne.0.0) then dens_new(i,j) = dens_new(i+1,j) elseif(uvel(i-1,j).ne.0.0.or.vvel(i-1,j).ne.0.0) then dens_new(i,j) = dens_new(i-1,j) elseif(uvel(i,j+1).ne.0.0.or.vvel(i,j+1).ne.0.0) then dens_new(i,j) = dens_new(i,j+1) else dens_new(i,j) = 0.0d0 endif endif enddo enddo c -- 濃度の源 do i=1,id do j=1,jd if (i.ge.(isrc-1).and.i.le.(isrc+1).and. & j.ge.(jsrc-1).and.j.le.(jsrc+1)) then dens_new(i,j) = 1.0d0 endif enddo enddo return end subroutine kyoukai c*********************************************************************** c -- 移流拡散の計算を行うsubroutine c*********************************************************************** c これにより、dens_oldとuvel,vvel,dens等を元に、dens_newが計算される。 subroutine iryu_kaku(id,jd,dx,dt,diff,uvel,vvel,dens_old,dens_new) integer id,jd real dx,dt real diff real dens_new(id,jd) real dens_old(id,jd) real uvel(id,jd) real vvel(id,jd) do i=2,id-1 do j=2,jd-1 dens_new(i,j)=dens_old(i,j) & - dt*uvel(i,j) & *(dens_old(i+1,j)-dens_old(i-1,j))/(2.0d0*dx) & - dt*vvel(i,j) & *(dens_old(i,j+1)-dens_old(i,j-1))/(2.0d0*dx) enddo enddo do i=2,id-1 do j=2,jd-1 dens_new(i,j)=dens_new(i,j) & + dt*diff & *(dens_old(i+1,j)-2.0d0*dens_old(i,j)+dens_old(i-1,j)) & /dx/dx & + dt*diff & *(dens_old(i,j+1)-2.0d0*dens_old(i,j)+dens_old(i,j-1)) & /dx/dx dens_new(i,j)=max(dens_new(i,j),0.0d0) dens_new(i,j)=min(dens_new(i,j),1.0d0) enddo enddo return end subroutine iryu_kaku