Output issue with program build

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

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/(nnprev) ip1=2nprev ip2=ip1n ip3=ip2nrem 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=2ifp1 theta=isign6.28318530717959d0/(ifp2/ip1) wpr=-2.d0dsin(0.5d0theta)**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=wrwpr-wiwpi+wr wi=wiwpr+wtempwpi+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

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

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)+2i-1, as i2 is used Also, the vector NN looks as if it should be integer4, rather than real8 ?

John

subroutine fftdata (db,data,n,m,n2,n3,nn) 
integer*4 n,m,n2,n3
real*8    db(n,m),data(n3)
integer*4 nn(2) 
!
integer*4 i,j
!
   do j=1,m 
     do i=1,n 
       data(n2*(j-1)+2*i-1)=db(i,j)
!   Note: this defines every second value of data, which agrees with error below       
     enddo 
   enddo 
   call fourn (data,nn,2,1)    ! <---- looks like one issue is in the subroutine 
   return 
end 

subroutine fourn (data,nn,ndim,isign) 
integer*4 ndim,isign
real*8    data(*)
INTEGER*4 nn(ndim)
! 
real*8    wr,wi,wpr,wpi,wtemp,theta,tempr,tempi, two_pi
integer*4 ntot,idim,nprev,n,nrem,ip1,ip2,ip3,i2rev,       &
          i1,i2,i3,i3rev,ibit,ifp1,ifp2,k1,k2
!
   two_pi = 8 * atan(1.0d0)   ! this could be a better constant
!
   ntot=1 
   do 11 idim=1,ndim 
     ntot=ntot*nn(idim) 
11 continue 
   nprev=1 
   do 18 idim=1,ndim 
     n     =nn(idim)                        ! ???  is nn real*8 or integer*4
     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 : is expected !!
             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)   ! check integer arithmetic in (ifp2/ip1)
       wpr  =-2.d0* SIN (0.5d0*theta)**2            ! use generic SIN
       wpi  = SIN (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)   ! why use sngl !!!!
             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 

Was data defined as complex*8 elsewhere ?

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

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.

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

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:

!read input gravity data
      open(5,file=fin)
      read(5,*) ((ba(i,j),i=1,n),j=1,m) ! in mgal
      close(5)

      z0=dabs(z0)
      delh=z0+zb
      delg=tpig*rho*delh
      ba=ba+delg

      call fftdata(ba,data0,n,m,n2,n3,nn)
      data0=flt*data0*exp(wd*z0)/tpig/rho
      do kitr=1,itr
        print*,' '
        print*,'iteration ', kitr
        print*,' '
        tdata=data0
        dep0=dep
        if(kitr.eq.1) goto10
        do io=2,nmax
          ba=dep**io
          call fftdata(ba,data,n,m,n2,n3,nn)
          data=data*wd**(io-1)/fac(io) *flt
          tdata=tdata-data
          big=0.d0
          do k=1,n3,2
            tt=data(k)**2+data(k+1)**2
            if(tt.gt.big) big=tt
          end do
          if(io.eq.2) then
            big2=big
          else
            bign=big
            print*,big2,bign
            r=bign/big2
          write(*,'(a4,i2,a16,e8.3)')'n=',io,'  Sn/S2 ratio =',dsqrt(r)
          if(r.lt.small2) goto 10
          endif
          end do
        end do
!10      continue
        call fftinv(dep,tdata,n,m,n2,n3,nn)
        dep=dep/dnm
        call deprms(dep-dep0,n,m,rms2)
        if(rms2.gt.rms1) then
          print*,' rms1=',rms1
          print*,' rms2=',rms2
          stop ' diverge !!! and stop '
        elseif(rms2.lt.rms) then
          goto 20
        elseif(kitr.ge.itr) then
          print*,' rms1=',rms1
          print*,' rms2=',rms2
          stop ' too slow converging.. please check it again !'
        else
          rms1=rms2
          print*,' rms=',rms2
        endif
      enddo
20    write(*,100) kitr,rms2
100       format(1x,' inversed after ',i2,' iteration !!!',/, &
         ' RMS is expected to be ',f6.2,'(in meter)')

      dep=dep-z0
      open(5,file=fdep)
      write(5,'(256f10.3)') ((dep(i,j),i=1,n),j=1,1) ! in km <---- output grid starts to fill and then stops
      close(5)

This is a program, published through Computers & Geosciences, and so it is almost certain to have issues. That aside, the basic code compiles but clearly errors are creeping in at run-time. What I am not sure about is why the output grid would start to fill in and then not complete? Really annoying that I cannot upload the whole code to show the lnkages (it is not very big !)

Cheers

Lester

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

      open (unit=25,file=fdep) 
      write (25,'(i6,f10.3)') ((i,dep(i,j),i=1,n),j=1,1) ! in km <---- output grid starts to fill and then stops 
      close (unit=25) 
      write (*,*) 'File ',trim(fdep),' complete'

2560 characters is a long text record to use. Excel can not read it (256 character limit) and I think NOTEPAD is limited to 1024 characters, so you should check the way you confirm the error.

I think I have also indicated, don't use a unit number < 10, as you are asking for trouble.

Also: if(r.lt.small2) goto 10 might be better changed to if (r.lt.small2) EXIT as 10 was the end statement of the outer loop and so should continue itterating. (similarly for goto 20)

John

Does specifying the record length as below help?

open (unit=25,file=fdep, recl=2560)

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

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

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

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=2n, n3=2n*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).

      subroutine deprms(ht,n,m,rms)
      real*8 ht(n,m),rms
      do i=1,n
      do j=1,m
        rms=rms+ht(i,j)**2
      enddo
      enddo
      rms=dsqrt(rms/n/m)
      return
      end

Not sure who wrote DEPRMS - this whole program is published:

Young Hong Shin,
, Kwang Sun Choi, Houze Xu. (2006) Three-dimensional forward and inverse models for gravity fields based on the Fast Fourier Transform. Computers & Geosciences, 32, 727–738

Got the output to work with:

write(26,'(f10.3)') ((dep(i,j),i=1,n),j=1,m) ! in km - j=1,m needed to scan the rest

The result is as expected from the supplied output data. The only tweaks needed to the code are for user input to define the input grid geometry rather than hard-coded into the program as it is now.

Well seems to work on the synthetic data; need to see how it performs on 'real' data!

Lester

A test with real data is not working, throwing up a run-time error like before

Floating-point overflow via calls to FOURN and FFTDATA 004027a0 FOURN [+022c] line 264 004022f0 FFTDATA [+00da] line 189 (call fourn) 00401000 main [+0507] line 55 (call fftdata)

The data are a bit over the top in precision so could well be the issue:

76.564239502 76.6045761108 76.5740585327 76.4089889526 76.0467681885 75.442855835 74.5952987671 73.5648574829 72.4396896362 etc

It gets through some of the iteration

iteration 1

rms= 28115.6448491

iteration 2

   3.340443934368E+24         9.981365633954E+31

n= 3 Sn/S2 ratio =.547E+04 3.340443934368E+24 3.879393911229E+39 n= 4 Sn/S2 ratio =.341E+08 3.340443934368E+24 1.112922427321E+47 n= 5 Sn/S2 ratio =.183E+12 3.340443934368E+24 3.784955186161E+54 n= 6 Sn/S2 ratio =.106E+16 3.340443934368E+24 6.912957069900E+61 n= 7 Sn/S2 ratio =.455E+19

ress RETURN to close window . . .

at most it should be rounded to 2 decimal places I guess

The main thing is that the program will work, it now seems a case of tweaking the precision settings and/or data.

To view big data files use Notepad++ it does not have the limitation of Excel!

Lester

Lester,

Your transform does not look to be converging, by the numbers quoted. May be this free method doesn't work with real data !!

John

Think I will go back and redo the subroutines from scratch rather than try to figure the logic of the published code. With a program that is a computation of a mathematical equation, it is important to document all the key elements. Anyway, it is good practice !

The subroutine FOURN was published in Rick Blakely's book on 'Potential Theory and Gravity and Magnetic Applications'

One of the perils of trying to use published code I guess.

Published ? Did you download it or type it in. If typed, there are lots of possibilities of typing errors and also mixing characters with bad fonts, eg 1,i,I,l (one, lower case i, upper case I or lower case L) or o/0 or 5/S or 8/B or 9/g. The list could go on. You could try typing it in twice ( or get someone else to) then use FC or Plato compare to see the result.

John

It was downloaded code from the journal. I am pretty sure the subroutine FOURN is the same as that shown in the Blakely's book since the variable name declarations are identical. The program code did not cite the original source of the FOURN subroutine code (which they should have done).

Will have another look to see. At least it is a fairly short program!

I do have a Matlab code (again published) which does pretty much the same thing, but looks very cumbersome as the summation is written out in full (n=2 to n=10), so in theory it would not be hard to code it in F95.

I am wondering if there is an issue with the way the existing program is computing the RMS variance for convergence/divergence?

Cheers Lester

I think there could be a problem with initialising (zeroing) variables. The following identifies two I think should be done. There are probably many more.

      subroutine deprms(ht,n,m,rms) 
!
      real*8 ht(n,m),rms 
!
      rms = 0     ! initialise rms
!
      do i=1,n 
       do j=1,m 
        rms=rms+ht(i,j)**2 
       end do 
      end do 
      rms = sqrt ( rms / (n*m) )
      return 
      end


subroutine fftdata (db,data,n,m,n2,n3,nn) 

integer*4 n,m,n2,n3 
real*8    db(n,m),data(n3) 
integer*4 nn(2) 
! 
integer*4 i,j 
! 
   data = 0     ! initialise data
!
   do j=1,m 
     do i=1,n 
       data(n2*(j-1)+2*i-1)=db(i,j) 
     end do 
   end do 

   call fourn (data,nn,2,1)

   return 
end 

Also I would recommend than NN be declared as integer, rather than real. Have a look at how I have initialised data and rms. If this should be done, there are probably others like this.

If these changes are required, you would have to question the reliability of the original published code, as a description of rubbish comes to mind.

John