ISPACK F90 を使うことの利点として以下のような点があります.
以下では簡単なプログラムを例に, ISPACK F90 モジュールの使い方を 説明します.
% frt -o prog n2pack.f90 program.f90 -L/usr/local/lib -lisp
境界条件は x,y 方向とも周期的であるとします.
以下が ISPACK F90 を使ったこの方程式を解くためのプログラムです.
! Sample program for gtool_history/gtool4 and ISPACK 2001/09/23 S.Takehiro
!
! Solving 2-D non-divergent (barotropic) fluid system
! d\zeta/dt + J(\psi,\zeta) + beta d\psi/dx = \nabla^2\zeta,
! \nabla\psi = \zeta
!
program baro
use n2pack ! モジュール使用宣言
integer,parameter :: im=32, jm=32 ! 格子点の設定(X,Y)
integer,parameter :: km=15, lm=15 ! 切断波数の設定(X,Y)
double precision, dimension(0:jm-1,0:im-1) :: psi, zeta ! 格子データ(ψ,ζ)
double precision, dimension(-lm:lm,-km:km) :: spsi,szeta ! スペクトルデータ(ψ,ζ)
double precision, dimension(0:jm-1,0:im-1) :: x, y ! 座標変数(X,Y)
double precision, parameter :: xmin=0.0, xmax=1.0 ! 計算範囲(X)
double precision, parameter :: ymin=0.0, ymax=1.0 ! 計算範囲(Y)
double precision, parameter :: beta=1.0, nu=1e-8 ! パラメター
double precision, parameter :: dt=1e-1 ! 時間積分刻み
integer, parameter :: nt=1000, ndisp=10 ! 時間ステップ数
double precision, parameter :: x1=(xmax-xmin)*0.5 ! 初期渦用パラメター
double precision, parameter :: y1=(ymax-ymin)*0.5 ! 初期渦用パラメター
double precision, parameter :: sigma=0.05, zeta0=0.2 ! 初期渦用パラメター
integer :: i,j, k,l
!-------------------- 座標設定 -----------------------
do i=0,im-1
do j=0,jm-1
x(j,i)= (xmax-xmin)/im*i
y(j,i)= (ymax-ymin)/jm*j
enddo
enddo
call n2initial(im,jm,km,lm,xmax-xmin,ymax-ymin) ! スペクトル変換初期設定
!------------------- 初期値設定 ----------------------
zeta= zeta0*exp(-((x-x1)**2+(y-y1)**2)/(2*sigma**2)) ! 初期渦分布
szeta = s_g(zeta) ! スペクトルデータへ変換
szeta(0,0)=0.0 ! 渦度の平均値は 0.
psi=s_laplainv_s(szeta) ! 流線へ変換
!-------------------- 時間積分 -----------------------
do it=1,nt
szeta = szeta + &
dt*( - s_jac_ss(spsi,szeta) &
-beta*s_dx_s(spsi) + nu*s_lapla_s(szeta) )
spsi = s_laplainv_s(szeta)
if(mod(it,ndisp) .eq. 0)then ! 出力処理
call output
endif
enddo
stop
contains
subroutine output
write(6,*) 'it=', it
write(6,*) g_s(spsi)
end subroutine output
end program baro
実際にはスペクトルデータの時間積分を行っているのですが, その時間発展の式は, n2pack モジュールの関数を使うことで あたかも実空間の方程式をそのまま書き下している形に書くことができます.
このプログラムを実行するには n2pack.f90 と一緒にコンパイルするのが 一番簡単です. プログラムと n2pack.f90 を同じディレクトリにおいて コンパイルし ISPACK ライブラリをリンクすれば OK です. 例えば Linux 上の Fujitsu Fortran ならば以下のようになります.
% frt -o baro.out n2pack.f90 baro.f90 -L/usr/local/lib -lisp-L オプションに続くディレクトリは ISPACK ライブラリの在処を 表しています.