Fortranの一次元浅水流方程式の計算プログラム

fortranで一次元浅水流方程式を行うプログラムを作成したいのですが、初学者すぎてわからないことがあり、以下の主プログラムとサブルーチンを組み合わせると求めたいプログラムが完成するらしいのですが、どうすればいいかわかりません。

～１次元浅水流方程式を行う主プログラム～

use subprogs
implicit none
real(8), allocatable :: x(:), u(:), v(:), w(:), h(:)
real(8), allocatable :: p(:), q(:), pn(:), qn(:)
real(8) dt, xl, dx, fr, gr, ho, dh
integer n, itr, itrmax, i, pintv
call set_init(n, dt, xl, dx, fr, gr, itrmax, pintv, ho, dh, &
x,u,v, w,h,p, a,pn, an)!変数値,初期条件などの設定
call print_uh(x,h,u, gr,n,1) !出カファイルを開き,初期条件を出力
do itr = 1, itrmax
call chk_cno(n, dx, dt, v, w)

doi=2,n-1!内部格子点のpn,qn を特性曲線法により定める

call cm1d(pn, p, i, v, dt, dx)
call cm1d(qn, q, i, w, dt, dx)
enddo
call bc_thru(pn, qn, p, 9,v,w, n, dt,dx)!境界のpnとqnを定める
call pq2uhvw(pn,qn, gr,u,h,v,w, n)! pn,anからu,h,v,wを定める
p(:)=pn(:)!全格子点のpを更新する
q(:)= qn(:)!全格子点のgを更新する
if (mod(itr, pintv) ==0) call print_uh(x,h, u, gr, n,0)!出力
enddo
call print_uh(x, h, u, gr, n, -1)
!出力ファイルをcloseする
deallocate(x, u, v, w, h, p, q, pn,qn)!割付解除する
end program main

～サブルーチンset_initの演算の主要部分～

open(10, file = 'data.d')!入力ファイルを開く
read(10, *) n, itrmax, pintv!格子点数n,反復回数,出力間隔
read(10, *) dt, xl, fr, gr, h0
close(10)


dx = xl / dble(n - 1) !d×の設定
dh = 0.1d0 * h0 !初期水面形設定のためのdhを設定
re (x(n), u(n), v(n), w(n), h(n), p(n), q(n), pn(n),qn(n))
x(:)=(/(dx*dble(i-1),i=1, n)/)!座標を設定
u(:)= sqrt(gr* h0)* fr!フルード数frを用いて流速を設定

h(:) = h0 + dh *exp(-(x(:)-0.5d0*xl)**2/ 1.0d2)!初期水面形
call uh2pqvw(u,h,gr,p,q,v,w,n)!u,hからp,q,v,wを定める

～サブルーチンprint_uhの演算の主要部分～

subroutine print_uh(x, h, u, gr, n, fopen)
!....(変数の宣言文)
if (fopen == 1) then  !fopen=1ならファイルをopenする
open(20, file = 'uh.d')
else if (fopen == -1)then   !fopen=-1ならファイルをcloseしてreturn
close (20)
return
endif
do i=1,n!各格子点上のx,h,u,フルード数の出力
write(20,'(10e16.8)')x(i), h(i), u(i), u(i) I sqrt(gr * h(i))
enddo
write(20,*)    !gnuplotによる描画のため空行を入れておく
end subroutine print_uh

～サブルーチンchk_cnoの演算の主要部分～

cno1 = maxval(abs(v(:))* dt / dx)
cno2 = maxval(abs(w(:)) x dt I dx)
cno = max(cno1,cno2)
if (cno >= 1) then
write(*, *) 'stop, cno >= 1, cno = ', cno
stop
endif

～サブルーチンcmldの演算の主要部分～

cno=v(i)*dt/dx　　!クーラン数を計算
f(cno>= 0.0d0) then　　　!特性曲線の出発点位置に応じて空
pn(i)= p(i)-cno*(p(i)-p(i-1))
else
pn(i)=p(i)-cno*(p(i+1)-p(i))
endif

～サブルーチンbc_thruの演算の主要部分～

do i = 1, n, n - 1 !このループでは,iは1とnという値のみをとる
pn(i)=p(i) !前ステップの値をデフォルト値として設定
qn(i)= q(i) !同上
cno1 =v(i)*dt/ dx!クーラン数を計算
cno2 = w(i) * dt/dx!同上
if (i == 1) then !上流端(特性速度が負ならp,q算
if(cno1 <0.0d0)pn(i)=p(i)-cno1 *(p(i+1)-p(i))
if (cno2 < 0.0d0) qn(i) = q(i)- cno2*(q(i+1)-q(i))
elseif(i== n)then! 下流端(特性速度が正ならpn,qnを計算)
if (cno1 > 0.0d0)pn(i)=p(i)-cno1*(p(i )-p(i-1))
if (cno2 > 0.0d0) qn(i)= q(i)-cno2*(q(i)-q(i-1))
endif
enddo