! Sample program for gt4_history/gt4f90io and ISPACK 2001/02/27 S.Takehiro ! 2003/08/21 M.Odaka ! ! Solving 2-D Boussinesq fluid system ! d\zeta/dt + J(\psi,\zeta) = Ra dT/dx + \nabla\zeta ! dT/dt + J(\psi,T) - d\psi/dx = \nabla T ! \nabla\psi = \zeta ! psi = zeta = T = 0 at y=0,1 ! program benard1 use c2pack use gt4_history integer, parameter :: km=16, lm=8 ! 切断波数の設定(X,Y) integer, parameter :: im=64, jm=16 ! 格子点の設定(X,Y) double precision, dimension(0:jm,0:im-1) :: psi, temp, zeta double precision, dimension(-km:km,lm) :: wpsi, wtemp, wzeta double precision, dimension(0:jm,0:im-1) :: x, y double precision, parameter :: xmin=0.0, xmax=4.0 double precision, parameter :: ymin=0.0, ymax=1.0 double precision, parameter :: Ra=1.0e4 double precision, parameter :: dt=1e-4 integer, parameter :: nt=5000, ndisp=500 integer :: i,j, k,l tinit = 0.0 ! 初期時刻設定 do i=0,im-1 do j=0,jm x(j,i)= (xmax-xmin)/im*i y(j,i)= (ymax-ymin)/jm*j enddo enddo temp = 0.0 temp(jm/2,im/2) = 0.01 call c2initial(im,jm,km,lm,xmax,ymax) wtemp = wsin_g(temp) ! 初期値設定 wpsi = 0.0 ; wzeta = 0.0 ! 初期値設定 temp = temp + 1-y ! データ出力用 call output_gtool4_init ! ヒストリー初期化 do it=1,nt wtemp = wtemp + dt*( -jac_ws(wpsi,wtemp) + dx_ws(wpsi) + nabla_ws(wtemp) ) wzeta = wzeta + & dt*( - jac_ws(wpsi,wzeta) & + Ra*dx_ws(wtemp) + nabla_ws(wzeta) ) wpsi = nabla_inv_ws(wzeta) if(mod(it,ndisp) .eq. 0)then temp = grid_ws(wtemp)+1-y psi = grid_ws(wpsi) call output_gtool4 endif enddo call output_gtool4_close stop contains subroutine output_gtool4_init call HistoryCreate( & ! ヒストリー作成 file='benard_1.nc', title='convection in rotating annulus', & source='Sample program of gt4_history/gt4f90io', & institution='GFD_Dennou Club davis project', & dims=(/'x','y','t'/), dimsizes=(/im,jm+1,0/), & longnames=(/'X-coordinate','Y-coordinate','time '/),& units=(/'1','1','1'/), & origin=real(tinit), interval=real(ndisp*dt) ) call HistoryPut('x',x(1,0:im-1)) ! 変数出力 call HistoryPut('y',y(0:jm,1)) ! 変数出力 call HistoryAddVariable( & ! 変数定義 varname='psi', dims=(/'x','y','t'/), & longname='stream function', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='temp', dims=(/'x','y','t'/), & longname='temperature', units='1', xtype='double') end subroutine output_gtool4_init subroutine output_gtool4 write(6,*) 'it = ',it call HistoryPut('t',real(it*dt)) call HistoryPut('psi',trans_g(psi)) call HistoryPut('temp',trans_g(temp)) end subroutine output_gtool4 subroutine output_gtool4_close call HistoryClose end subroutine output_gtool4_close end program benard1