! Sample program for gtool_history/gtool4 and ISPACK 2001/02/27 S.Takehiro ! ! 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 benard2 use c2pack use gtool_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 type(GT_HISTORY) :: hst_psi, hst_temp integer :: i,j, k,l 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='psi.nc', title='convection in rotating annulus', & source='Sample program of gtool_history/gtool4', & 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=0.0, interval=real(ndisp*dt), & history=hst_psi ) call HistoryPut('x',x(1,0:im-1),hst_psi) ! 変数出力 call HistoryPut('y',y(0:jm,1),hst_psi) ! 変数出力 call HistoryAddVariable( & ! 変数定義 varname='psi', dims=(/'x','y','t'/), & longname='stream function', units='1', xtype='double',& history=hst_psi ) call HistoryCreate( & ! ヒストリー作成 file='temp.nc', title='convection in rotating annulus', & source='Sample program of gtool_history/gtool4', & institution='GFD_Dennou Club davis project', & dims=(/'1','1','1'/), dimsizes=(/im,jm+1,0/), & longnames=(/'X-coordinate','Y-coordinate','time '/),& units=(/'m','m','s'/), & origin=0.0, interval=real(ndisp*dt), & history=hst_temp ) call HistoryPut('x',x(1,0:im-1),hst_temp) ! 変数出力 call HistoryPut('y',y(0:jm,1),hst_temp) ! 変数出力 call HistoryAddVariable( & ! 変数定義 varname='temp', dims=(/'x','y','t'/), & longname='temperature', units='1', xtype='double',& history=hst_temp) end subroutine output_gtool4_init subroutine output_gtool4 write(6,*) 'it = ',it call HistoryPut('psi',trans_g(psi),hst_psi) call HistoryPut('temp',trans_g(temp),hst_temp) end subroutine output_gtool4 subroutine output_gtool4_close call HistoryClose(hst_psi) call HistoryClose(hst_temp) end subroutine output_gtool4_close end program benard2