forums.silverfrost.com

[arctica]

Hello,

I have compiled and built an old fixed-format Fortran code, which runs without any run-time errors. It reads an input grid, outputs a statistics file but fails to generate the full output computation grid.

1. Reads in gravity data grid OK
2. Computes the depth from gravity inversion OK (stats generated)
3. Writes the inversion grid depths - file opened but limited values written
4. writes the statistics file OK

As far as I can tell all but point 3 are running as expected. The output grid is opened and data is partly written to it but on 256 elements and not the 65536 values to match the input matrix geometry.

Not sure what could be the problem, any ideas what to test/try would be helpful.

Cheers
Lester

[brucebowler]

No where near enough info to provide any sort if intelligent advice. Can you post the relevant bits of code, and the values you get and expect?

How about adding some diagnostic writes to make sure things are behaving as expected at prior points in the code?

Have you run it in the debugger to see if variables have the values you expect?

[arctica]

Hi

I put checkmate on and forced the run-time errors and copied the line codes indicated:

FOURN - in file 3dinvgrav.f at line 262 [+05ce]
tempi=data(i3+1)

FFTDATA - in file 3dinvgrav.f at line 187 [+02bf]
call fourn(data,nn,2,1)

MAIN# - in file 3dinvgrav.f at line 55 [+04be]
call fftdata(ba,data0,n,m,n2,n3,nn)

subroutine fftdata(db,data,n,m,n2,n3,nn)
real*8 db(n,m),data(n3),nn(2)
do j=1,m
do i=1,n
data(n2*(j-1)+2*i-1)=db(i,j)
enddo
enddo
call fourn(data,nn,2,1) <---- looks like one issue is in the subroutine
return
end

subroutine fourn(data,nn,ndim,isign)
real*8 wr,wi,wpr,wpi,wtemp,theta
real*8 nn(ndim),data(*)
ntot=1
do 11 idim=1,ndim
ntot=ntot*nn(idim)
11 continue
nprev=1
do 18 idim=1,ndim
n=nn(idim)
nrem=ntot/(n*nprev)
ip1=2*nprev
ip2=ip1*n
ip3=ip2*nrem
i2rev=1
do 14 i2=1,ip2,ip1
if(i2.lt.i2rev)then
do 13 i1=i2,i2+ip1-2,2
do 12 i3=i1,ip3,ip2
i3rev=i2rev+i3-i2
tempr=data(i3)
tempi=data(i3+1) <------- gives an undef var error here
data(i3)=data(i3rev)
data(i3+1)=data(i3rev+1)
data(i3rev)=tempr
data(i3rev+1)=tempi
12 continue
13 continue
endif
ibit=ip2/2
1 if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then
i2rev=i2rev-ibit
ibit=ibit/2
go to 1
endif
i2rev=i2rev+ibit
14 continue
ifp1=ip1
2 if(ifp1.lt.ip2)then
ifp2=2*ifp1
theta=isign*6.28318530717959d0/(ifp2/ip1)
wpr=-2.d0*dsin(0.5d0*theta)**2
wpi=dsin(theta)
wr=1.d0
wi=0.d0
do 17 i3=1,ifp1,ip1
do 16 i1=i3,i3+ip1-2,2
do 15 i2=i1,ip3,ifp2
k1=i2
k2=k1+ifp1
tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1)
tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2)
data(k2)=data(k1)-tempr
data(k2+1)=data(k1+1)-tempi
data(k1)=data(k1)+tempr
data(k1+1)=data(k1+1)+tempi
15 continue
16 continue
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
17 continue
ifp1=ifp2
go to 2
endif
nprev=n*nprev
18 continue
return
end

So it is tripping over in a couple subroutines - part of the issue may be that this is Fortran 77 code?

No way to upload the whole code, not that it that big, so not easy to illustrate.

Lester

[arctica]

The program is running as the iteration converges so it should write out the final grid result (a bunch of z-values)

iteration 1

rms= 2.06079503989

iteration 2

100534.057240 256.006955673
n= 3 Sn/S2 ratio =.505E-01
100534.057240 2.11354468545
n= 4 Sn/S2 ratio =.459E-02
rms= 3.949157578910E-02

iteration 3

119998.989518 353.642016567
n= 3 Sn/S2 ratio =.543E-01
119998.989518 5.38592140571
n= 4 Sn/S2 ratio =.670E-02
119998.989518 0.214114260575
n= 5 Sn/S2 ratio =.134E-02
inverted after 3 iteration !!!
RMS is expected to be 0.01(in meter)

Press RETURN to close window . . .

I haven't looked in detail at the code, but it would be easier if we could upload whole program codes here, but that is not possible.

Cheers
Lester

[JohnCampbell]

Lester,

I reviewed your code, providing declatrations and indents for DO and IF.
You might like to check what I have done.
fftdata defines every second value, n2*(j-1)+2*i-1, as i*2 is used
Also, the vector NN looks as if it should be integer*4, rather than real*8 ?

John

[arctica]

There is no complex declaration in the code (which is F77)

parameter (n=256,m=128,n2=512,n3=65536) !n2=2*n, n3=2*n*m
implicit real*8 (a-h,o-z)
real*8 dep(n,m),ba(n,m),dep0(n,m)
real*8 data(n3),nn(2),tdata(n3),wd(n3),flt(n3),data0(n3)
real*8 rms,rms1,rms2
character*15 fin,fdep

data nn/n,m/
data zb,z0,delx,dely/-7.d0,5.d0,2.d0,3.d0/ !base and mean depths, and grid separations (km)
data tpig,rho/41.9277d0,1.d0/ ! density (in gram/cubic cm)
data itr,nmax,small/50,20,5.d-3/
data rms,rms1/1.d-2,1.d10/
data fin,fdep/'grav.dat','depth.dat'/
data kp,w1,w2/1,0.3,0.4/ ! power of cosine filter & filter parameter
small2=small**2
dnm=dfloat(n*m)
p=datan(1.d0)*8.d0

call wdata(delx,dely,wd,n,m,n2,n3)
wd=wd*p
call filter(w1,w2,n3,wd,flt)
flt=flt**kp

open(22,file='inv_stats.rlt')
!c read input gravity data
open(25,file=fin, status='OLD')
read(25,*) ((ba(i,j),i=1,n),j=1,m) ! in mgal
close(25)

I wonder if it is an issue with the F77 code compiling under F95 (I don't have F77 to test).

Bit tricky to explain without showing the entire code, which we cannot do here!

Cheers
Lester

[LitusSaxonicum]

It's never an issue with Fortran 77 code, as the entireity of Fortran 77 is in the Fortran 95 standard - even if parts of it are "deprecated". Fortran 77 allows complex.

[JohnCampbell]

Lester,

I'm not sure what your problem is.
* You indicated that point 3 was the main outstanding issue, but no clue to that. I would check that DO loops are itterating.
* The code example you showed, only defined every second value in array "data". This looks a problem for the DO 17 loop. My comment about complex, was based on the way you missed every second value, possibly defining the real component only.
* the do loops 15,16,17 look like a complex way to manipulate parts of the data array. If this is old code, (certainly pre-F90, not sure if pre-F77) there is a good chance that these loops would have steped at least once, but with F90 might not do an itteration.

For example : the following loop would not itterate in F90, but would have in older versions of Fortran.
I1 = 1 ; I2 = 0
DO i = i1,i2
END DO
You might need to check this. There are compiler options to force at least 1 itteration of the loop.

John

[arctica]

Hi John

According to the code it should write out a matrix (grid), the file is opened to write to and values added but it does not complete:

[JohnCampbell]

I'd try the following and make sure you get the end message:

[davidb]

Does specifying the record length as below help?

[arctica]

Hi John,

I added the extra code for the check, and the

write (*,*) 'File ',trim(fdep),' complete'

The comment comes up as File depth.dat complete - so that works

The depth.dat file is only 256 entries. I use different editors to view files like Notepad++ and Surfer's editor/spreadsheet so big files are not a problem.

Lester

[arctica]

David

Adding:

open(26,file=fdep, recl=2560)

Did not change anything much, still only get 256 values. Also tried setting a separate channel number to see if that made a difference but not.

Lester

[JohnCampbell]

recl does not help with a formatted output file. It can make things worse.

Isn't 256 is the correct number of output ?
write (25,'(i6,f10.3)') ((i,dep(i,j),i=1,n),j=1,1)

if you change to
write (25,'(i6,f10.3)') ((i,dep(i,j),i=1,n),j=1,m)
this might be what you want ?

or:
write (25,'(2i6,f10.3)') ((j,i,dep(i,j),i=1,n),j=1,m)

Also, who wrote:
call deprms (dep-dep0,n,m,rms2)

perhaps it might work, but is likely to cause stack problems. This is not old Fortran. I would never use this.

John

[arctica]

Hi John,

I checked one of the example output files and the dimensions are 256 x 128 (nx, ny) - counted the values in the file!

First line of the program:

parameter (n=256,m=128,n2=512,n3=65536) !n2=2*n, n3=2*n*m

Here 'n' and 'm' give the x and y dimensions, and I think n2 and n3 are used as part of the Fourier transform, to avoid boundary effects.

So what it appears to be doing is writing the first line and then stopping. The data are in scan-line orientation, so read from top-left to bottom-right - verified this using GMT (Generic Software Tools).