Fortranについてです。1㎜もわからないので教えてください。

C言語は少しやっていたのですが、Fortranは全く分かりません。
以下のプログラムについてなのですが、実行結果(プログラミングの下に記載）を理想形（一番下に記載）にしたいです。
プログラミング上部の7行部分を変えるだけでできると聞いたのですが、わかる方教えていただけないでしょうか？

do
write(*,'(a,$)')' {n x}==>'
read(*,*) n,x
write(*,'(9e14.6)') x,fj(n,x)
enddo

end
!}


!==============================================================================
function fj(m,x)
!# Abramowitz function J_m(x), Gauss4.
!# verified(Gauss<->asym) 08.9.15.
!# Accuracy upgraded (09.9.2)
!# Accuracy guaranteed down to 5th figure, see fjtest.doc (11.1.26)
!# xlim= 3.5! (18.1.18)
!------------------------------------------------------------------------------
parameter (ng=4)
dimension zz(ng),ww(ng)

!{ constants
pi= 3.14159265358979324
sqrtpi= 1.77245385090552
b2 = 0.6341764926910400E+00 ! asymptotic coeffcient
b3 = 0.5908179585223422E+00
zz(4)= 0.86113631
zz(3)= 0.33998104
zz(2)= -zz(3)
zz(1)= -zz(4)
ww(4)= 0.34785485
ww(3)= 0.65214515
ww(2)= ww(3)
ww(1)= ww(4)
xlim= 5.0 ! asymptotic expansion when x>xlim, revised 11.1.26
na= 32 ! number of asymptotic expansion
td= 5.0 ! upper limit of improper integral
!}

!{ special cases
if(m==-1)then
if(x<1e-8) then
write(*,*)'!*** invalid x.'
return
else if(x<1.0) then
nx= 20.0/x
else
nx= 32
endif
elseif(m==0)then
if(x<1e-6) then
fj= sqrtpi*0.5 -b2*x
return
else if(x<1e-3) then
fj= sqrtpi*0.5 -b2*x +x*(log(x)+0.5)
return
else if(x<1.0) then
nx= nint(30/x**0.66)
else
nx= 64
endif
else
nx= 64
endif
!}

!{ body
if(x<xlim)then ! Gauss4 integral
dt= td/nx
f= 0.0
do i= 1,nx
t2= dt*i
t1= t2 -dt
tp= (t2 +t1)*0.5
tq= (t2 -t1)*0.5
ff= 0.0
do k= 1,ng
t= tp+tq*zz(k)
ff= ff +exp(-t**2 -x/t)*t**m*ww(k)
enddo
ff= ff*tq
f= f +ff
enddo
fj= f
else ! asymptotic expansion
pow= 2.0/3.0
v= 3.0*(x/2)**pow
fac= sqrt(pi/3.0)*(v/3.0)**(m/2.0)*exp(-v)
f= 1.0
a= 1.0
i= 1
a1= a
a= (3*m**2+3*m-1)/12.0
f= f +a/v**i
do i= 2,na
k= i-2
a0= a1
a1= a
a= -(12*k**2+36*k-3*m**2-3*m+25)*a1 +0.5*(m-2*k)*(2*k+3-m)*(2*k+3+2*m)*a0
a= a/(12*(k+2))
f= f +a/v**i
enddo
fj= fac*f
endif
!}

end

!============================= eof ====================================


// 実行結果 //
現在のプログラムだと、以下のようになるのですが（このまま5.0まで実行していく）
{n x}==>1.0
0
0.000000E+00 0.500000E+00
{n x}==>1
0.1
0.100000E+00 0.426339E+00
{n x}==>0
0.2
0.200000E+00 0.506232E+00
{n x}==>1
0.3
0.300000E+00 0.323817E+00

理想形下のように一気に出したいです。
0.000000E+00 0.500000E+00
0.100000E+00 0.426339E+00
0.200000E+00 0.506232E+00
0.300000E+00 0.323817E+00
　　　　　　　　～
0.470000E+01 0.727832E-02
0.480000E+01 0.679401E-02
0.490000E+01 0.634450E-02
0.500000E+01 0.379055E+00

No.1 回答者： Tacosan 回答日時： 2018/12/07 00:14

C がわかるのなら, Fortran をちょっと調べればできるような気がするんだけど.

No.2 回答者： にんたまくん 回答日時： 2018/12/07 16:20

大学生？。

上の7行以下は、コメント、初期値、 fj(m,x)関数の計算をやっているだけみたいですよ。
数学関係の計算をやっているみたいですね。
C, fortranはやっていないけど、Python, matlabは動かせる人間が見ても、
このプログラムは、慣れた人が書いたようには思えません（失礼ながら）。

何人かによって、何回か、reviseもされている様な気配。。
計算アルゴリズムが分かっているなら、この解読するより、Cでいきなり組んだ方が早いかもしれませんね。

No.3 回答者： にんたまくん 回答日時： 2018/12/12 13:15

No.2 です。

余り良いプログラムではないですけど、以下の様に変更すると、動きました。
ただ、数値を確認すると、質問者の方が言っている理想形の値で、
0.500000E+01 0.379055E+00の時だけ異なった値になっています。
0.500000E+01 0.592706E-02となってしまいました。
他の数値は同じであることまでは確認できました。
これ以上は、7行目以下のアルゴリズムが分からないので、
私が回答できるのはここぐらいまでかと。ご参考まで。

n=1
x=0
do
!write(*,'(a,$)')' {n x}==>'
!read(*,*) n,x
write(*,'(9e14.6)') x,fj(n,x)
x=x+0.1
enddo

end
!}