!---------------------------------------------------------------------- ! Copyright (c) 2002--2026 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! Sample program for SPML(based on gt4f90io and ISPACK) ! ! 2002/08/19 S.Takehiro ! 2004/01/26 M.Odaka ! 2004/06/03 M.Odaka ! 2004/06/17 M.Odaka ! 2005/03/16 S.Takehiro ! 2026/07/02 S.Takehiro esc_module => et_module ! ! Solving 2-D Boussinesq fluid system (TT problem) ! d\zeta/dt + J(\psi,\zeta) = PrRt dT/dx + PrRc dC/dx + Pr\nabla\zeta ! dT/dt + J(\psi,T) - d\psi/dx = \nabla T ! dC/dt + J(\psi,C) - d\psi/dx = Tau \nabla C ! \nabla\psi = \zeta ! psi = zeta = T = C = 0 at y=0,1 ! program ddfcnv0 use et_module use gtool_history implicit none !---- 空間解像度設定 ---- integer, parameter :: km=21, lm=21 ! 切断波数の設定(X,Y) integer, parameter :: im=64, jm=32 ! 格子点の設定(X,Y) !---- 変数 ---- real(8) :: yx_Psi(0:jm,0:im-1) ! 格子データ(流線) real(8) :: yx_Temp(0:jm,0:im-1) ! 格子データ(温度/温度擾乱) real(8) :: yx_Comp(0:jm,0:im-1) ! 格子データ(組成/組成擾乱) real(8) :: yx_Zeta(0:jm,0:im-1) ! 格子データ(渦度) real(8) :: et_Psi(-km:km,0:lm) ! スペクトルデータ(流線) real(8) :: et_TempA(-km:km,0:lm) ! スペクトルデータ(温度擾乱,t+1) real(8) :: et_CompA(-km:km,0:lm) ! スペクトルデータ(組成擾乱,t+1) real(8) :: et_ZetaA(-km:km,0:lm) ! スペクトルデータ(渦度,t+1) real(8) :: et_TempB(-km:km,0:lm) ! スペクトルデータ(温度擾乱,t) real(8) :: et_CompB(-km:km,0:lm) ! スペクトルデータ(組成擾乱,t) real(8) :: et_ZetaB(-km:km,0:lm) ! スペクトルデータ(渦度,t) real(8) :: e_TempBndry(-km:km,2) ! 境界値 real(8) :: e_CompBndry(-km:km,2) ! 境界値 real(8) :: e_ZetaBndry(-km:km,2) ! 境界値 !---- 座標変数など ---- real(8), parameter :: xmin=0.0, xmax=2.0 ! 領域範囲(X 方向) real(8), parameter :: ymin=0.0, ymax=1.0 ! 領域範囲(Y 方向) !---- 時間積分パラメター ---- real(8), parameter :: dt=5d-6 ! 時間ステップ間隔 ! integer, parameter :: nt=5000, ndisp=250 ! 時間積分数, 表示ステップ integer, parameter :: nt=200000, ndisp=4000 ! 時間積分数, 表示ステップ !---- 物理パラメター ---- real(8), parameter :: Ra=-1.2e4 ! 温度レイリー数 real(8), parameter :: Rc= 1.0e4 ! 組成レイリー数 real(8), parameter :: Pr=1.0 ! プランドル数 real(8), parameter :: Tau=0.1 ! 組成拡散/熱拡散の比 integer :: it=0 ! DO 変数 !---------------- 座標値の設定 --------------------- call et_initial(im,jm,km,lm,xmin,xmax,ymin,ymax) ! ISPACK初期化 !------------------- 初期値設定 ---------------------- yx_Temp = 0.0 ; yx_Temp(jm/2,im/2) = 0.01 ! 温度擾乱場 yx_Comp = 0.0 yx_Psi = 0.0 yx_Zeta = 0.0 et_Psi = et_yx(yx_Psi) et_TempA = et_yx(yx_Temp) ; et_TempB = et_TempA et_CompA = et_yx(yx_Comp) ; et_CompB = et_CompA et_ZetaA = et_yx(yx_Zeta) ; et_ZetaB = et_ZetaA e_TempBndry = 0.0D0 e_CompBndry = 0.0D0 e_ZetaBndry = 0.0D0 call output_gtool_init ! ヒストリー初期化 call output_gtool ! 初期値出力 !------------------- 時間積分 ---------------------- do it=1,nt et_TempA = et_TempB + & dt*( - et_Jacobian_et_et(et_Psi,et_TempB) & + et_Dx_et(et_Psi) & + et_Lapla_et(et_TempB) ) et_CompA = et_CompB + & dt*( - et_Jacobian_et_et(et_Psi,et_CompB) & + et_Dx_et(et_Psi) & + Tau*et_Lapla_et(et_CompB) ) et_ZetaA = et_ZetaB + & dt*( - et_Jacobian_et_et(et_Psi,et_ZetaB) & + Pr*Ra*et_Dx_et(et_TempB) & + Pr*Rc*et_Dx_et(et_CompB) & + Pr*et_Lapla_et(et_ZetaB) ) call et_Boundaries(et_TempA,e_TempBndry,cond='DD') call et_Boundaries(et_CompA,e_CompBndry,cond='DD') call et_boundaries(et_ZetaA,e_ZetaBndry,cond='DD') et_Psi = et_LaplaInv_et(et_ZetAa) et_TempB = et_TempA et_CompB = et_CompA et_ZetaB = et_ZetaA if(mod(it,ndisp) .eq. 0)then ! 出力 call output_gtool endif enddo call output_gtool_close stop contains subroutine output_gtool_init call HistoryCreate( & ! ヒストリー作成 file='ddfcnv-et-0.nc', title='Double diffusive convection', & source='Sample program of spmodel', & institution='GFD_Dennou Club davis/spmodel 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) ) call HistoryPut('x',x_X) ! 変数出力 call HistoryAddattr('x','topology','circular') ! 周期属性 call HistoryAddattr('x','modulo',xmax-xmin) ! 周期属性 call HistoryPut('y',y_Y) ! 変数出力 call HistoryAddVariable( & ! 変数定義 varname='psi', dims=(/'x','y','t'/), & longname='stream function', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='zeta', dims=(/'x','y','t'/), & longname='voticity', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='temp', dims=(/'x','y','t'/), & longname='temperature', units='1', xtype='double') call HistoryAddVariable( & ! 変数定義 varname='comp', dims=(/'x','y','t'/), & longname='composition', units='1', xtype='double') end subroutine output_gtool_init subroutine output_gtool if ( Ra .ge. 0 )then yx_Temp = yx_et(et_TempA) + 1-yx_Y ! 出力は全温度場 else yx_Temp = yx_et(et_TempA) + yx_Y ! 出力は全温度場 endif if ( Rc .ge. 0 ) then yx_Comp = yx_et(et_CompA) + 1-yx_Y ! 出力は全組成場 else yx_Comp = yx_et(et_CompA) -yx_Y ! 出力は全組成場 endif yx_Psi = yx_et(et_Psi) yx_Zeta = yx_et(et_ZetaA) write(6,*) 'it = ',it call HistoryPut('psi',transpose(yx_Psi)) call HistoryPut('zeta',transpose(yx_Zeta)) call HistoryPut('temp',transpose(yx_Temp)) call HistoryPut('comp',transpose(yx_Comp)) end subroutine output_gtool subroutine output_gtool_close call HistoryClose end subroutine output_gtool_close end program