Inclassable déclaration à (1) fortran

À la ligne 99, la formule de gsurf(iel) me donne l'erreur:

unclassifiable statement at (1)

où l'on est au début de la ligne 99, lors de la compilation du programme. Toutes les suggestions sur la façon de résoudre ce problème?

program gravity

implicit none
real(8) Lx,Ly,sx,sy,xsphere,ysphere,r,A,rho1,rho2,dx,G1
integer np,nel,nelx,nely,i,nnx,nny,j,counter,nsurf,iel
real(8),dimension(:),allocatable :: xcgrid
real(8),dimension(:),allocatable :: ycgrid
real(8),dimension(:),allocatable :: xgrid
real(8),dimension(:),allocatable :: ygrid
real(8),dimension(:),allocatable :: rho
real(8),dimension(:),allocatable :: xsurf
real(8),dimension(:),allocatable :: ysurf
real(8),dimension(:),allocatable :: gsurf

nnx=101.
nny=101.
Lx=100.
Ly=100.
nelx=nnx-1.
nely=nny-1.
nel=nelx*nely
np=nnx*nny
sx=Lx/nelx
sy=Ly/nely
xsphere=50.
ysphere=50.
r=12.
nsurf=7  !number of gravimeters
G1=6.6738480*10**(-11) !m^3 kg^-1 s^-2

dx=Lx/(nsurf-1.)

!==========================================================

allocate(xgrid(np))
allocate(ygrid(np))

counter=0
do i=1,nnx
    do j=1,nny
    counter=counter+1   
    xgrid(counter)=dble(i-1)*sx
    ygrid(counter)=dble(j-1)*sy
    end do
end do

call write_two_columns(np,xgrid,ygrid,'grid_init.dat')
!==========================================================

allocate(xcgrid(np))
allocate(ycgrid(np))


counter=0
do i=1,nnx-1
    do j=1,nny-1
    counter=counter+1   
    xcgrid(counter)=dble(i-1)*sx+0.5*sx
    ycgrid(counter)=dble(j-1)*sy+0.5*sy
    end do
end do

call write_two_columns(np,xcgrid,ycgrid,'gridc_init.dat')
!==========================================================

allocate(rho(nel))

rho1=3000. !kg/m^3
rho2=3200. !kg/m^3

do i=1,nel  
    if (sqrt((xsphere-xcgrid(i))**2)+((ysphere-ycgrid(i))**2)<r) then
    rho(i)=3200.
    else 
    rho(i)=3000.
    end if
end do


call write_three_columns(nel,xcgrid,ycgrid,rho,'inclusion.dat')
!==========================================================

allocate(xsurf(nsurf))
allocate(ysurf(nsurf))

do i=1,nsurf
xsurf(i)=(i-1)*dx
ysurf(i)=ly
end do

call write_two_columns(nsurf,xsurf,ysurf,'surf_init.dat')
!==========================================================

allocate(gsurf(nel))

do i=1,nsurf
xsurf(i)=(i-1)*dx
ysurf(i)=ly
    do iel=1,nel
    gsurf(iel)=2.*G1*(((rho(iel)-rho1)*(y(iel)-ygrid))/((x(iel)-xgrid)**2.+(y(iel)-ygrid))**2.)))*sx*sy
    end do
end do

call write_two_columns (nel,ysurf,xsurf,gsurf,'gravity.dat')

deallocate(xgrid)
deallocate(ygrid)
deallocate(xcgrid)
deallocate(ycgrid)
deallocate(xsurf)
deallocate(ysurf)

end program"