Fortran code to read unformatted binary files 1

[Space_Fede]

Hello everyone,
I am trying to study the data of the gravity anomaly provided by the National Geospacial-Intelligence Agency.

In the following link you can find the Fortran code that should read the binary files in which the data I am interested in, are stored.
Link

When I try to compile the code, several errors appear. I am a Fortran newbie so any help would be appreciated.

For the record, I tried to compile it both with gfortran (in wsl) and cgywin.

Thanks in advance

 

[JamesONeil]

@mikrom yes I did but no way to find these files anymore

 

[mikrom]

@JamesONeil,
I still have the Fortran Source, but I have already deleted the 3 huge files.

 

[JamesONeil]

@mikrom I still have the 3 huge binary data files so if you have the source one I will be glad ! Right now I am trying to compile it after removing the shell part but I have a lot of "unclassifiable statement error" and it seems it is because of the fixed format in Fortran. So if you still have your source I am interested, thanks !

 

[mikrom]

@JamesONeil,
here is the extracted fortan code I mentioned above (I don't have the original with the shell part)

read_3files.f

c-----------------------------------------------------------------------
      program readmin1
      implicit real*8(a-h,o-z)
      parameter(exclud=9999.d0,
     $   dlatg=2.5d0/60.d0,dlong=2.5d0/60.d0,nrowsg=04320,ncolsg=08640)
         real*4 data(ncolsg)
      dimension sdat(ncolsg,2),stati(22)
c-----------------------------------------------------------------------
c
c   Read-in Gravity Anomaly 2.5x2.5 minute point value file.
c
      write(6,6001)
 6001 format(///15x,'Statistics of Gravity Anomaly Values  (mGal)',/)
      do i = 1, nrowsg
        read(1) (  data(j),j=1,ncolsg)
        rlat = 90.d0 - (i-1.d0)*dlatg - dlatg*0.5d0
        do j = 1, ncolsg
          clon = (j-1.d0)*dlong + dlong*0.5d0
c
c   Print a few values.
c
          if(i.eq.   1.and.j.lt.20) then  !  top row
            write(6,6101) rlat,clon,data(j)
 6101       format(5x,2f15.10,f15.5)
          endif  !  top row
          sdat(j,1) = data(j)
          sdat(j,2) = 1.d0
        enddo  !  j
        call stats(90.d0,0.d0,i,dlatg,dlong,nrowsg,ncolsg,sdat,exclud,0,
     $             stati)
      enddo  !  i
c-----------------------------------------------------------------------
c
c   Read-in xi Deflection of Vertical 2.5x2.5 minute point value file.
c
      write(6,6002)
 6002 format(///15x,'Statistics of xi DOV Values  (arc-second)',/)
      do i = 1, nrowsg
        read(2) (  data(j),j=1,ncolsg)
        rlat = 90.d0 - (i-1.d0)*dlatg - dlatg*0.5d0
        do j = 1, ncolsg
          clon = (j-1.d0)*dlong + dlong*0.5d0
c
c   Print a few values.
c
          if(i.eq.   1.and.j.lt.20) then  !  top row
            write(6,6101) rlat,clon,data(j)
          endif  !  top row
          sdat(j,1) = data(j)
          sdat(j,2) = 1.d0
        enddo  !  j
        call stats(90.d0,0.d0,i,dlatg,dlong,nrowsg,ncolsg,sdat,exclud,0,
     $             stati)
      enddo  !  i
c-----------------------------------------------------------------------
c
c   Read-in eta Deflection of Vertical 2.5x2.5 minute point value file.
c
      write(6,6003)
 6003 format(///15x,'Statistics of eta DOV Values  (arc-second)',/)
      do i = 1, nrowsg
        read(3) (  data(j),j=1,ncolsg)
        rlat = 90.d0 - (i-1.d0)*dlatg - dlatg*0.5d0
        do j = 1, ncolsg
          clon = (j-1.d0)*dlong + dlong*0.5d0
c
c   Print a few values.
c
          if(i.eq.   1.and.j.lt.20) then  !  top row
            write(6,6101) rlat,clon,data(j)
          endif  !  top row
          sdat(j,1) = data(j)
          sdat(j,2) = 1.d0
        enddo  !  j
        call stats(90.d0,0.d0,i,dlatg,dlong,nrowsg,ncolsg,sdat,exclud,0,
     $             stati)
      enddo  !  i
c-----------------------------------------------------------------------
      stop
      end
C
C
C
      SUBROUTINE STATS(TOPLAT,WSTLON,I,GRDN,GRDE,NROWS,NCOLS,DATA,
     $                 EXCLUD,ISIG,STAT)
C-----------------------------------------------------------------------
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*20 DLABEL(14)
      CHARACTER*20 SLABEL( 8)
      DIMENSION DATA(NCOLS,2)
      DOUBLE PRECISION DEXCLUD,DG,SD,STAT(22)
      SAVE
      DATA PI/3.14159265358979323846D+00/
      DATA DLABEL/'    Number of Values','  Percentage of Area',
     $            '       Minimum Value',' Latitude of Minimum',
     $            'Longitude of Minimum','       Maximum Value',
     $            ' Latitude of Maximum','Longitude of Maximum',
     $            '     Arithmetic Mean','  Area-Weighted Mean',
     $            '      Arithmetic RMS','   Area-Weighted RMS',
     $            '   Arithmetic S.Dev.','Area-Weighted S.Dev.'/
      DATA SLABEL/'       Minimum Sigma',' Latitude of Minimum',
     $            'Longitude of Minimum','       Maximum Sigma',
     $            ' Latitude of Maximum','Longitude of Maximum',
     $            'Arithmetic RMS Sigma','Area-wghtd RMS Sigma'/
C-----------------------------------------------------------------------
      IF(I.EQ.1) THEN
      DEXCLUD=EXCLUD
      DTR=PI/180.D0
      FOURPI=4.D0*PI
      DPR=GRDN*DTR
      DLR=GRDE*DTR
      CAREA=2.D0*DLR*SIN(DPR/2.D0)
      DO 10 K=1,22
      STAT(K)=0.D0
   10 CONTINUE
      STAT( 3)= DEXCLUD
      STAT( 6)=-DEXCLUD
      STAT(15)= DEXCLUD
      STAT(18)= 0.D0
      ENDIF
C-----------------------------------------------------------------------
      DLAT=TOPLAT-(I-1.D0)*GRDN-GRDN/2.D0
      COLATC=(90.D0-DLAT)*DTR
      AREA=CAREA*SIN(COLATC)
C-----------------------------------------------------------------------
      DO 110 J=1,NCOLS
      DLON=WSTLON+(J-1.D0)*GRDE+GRDE/2.D0
      DG=DATA(J,1)
      SD=DATA(J,2)
      IF(DG.LT.DEXCLUD) THEN
C-----------------------------------------------------------------------
      STAT( 1)=STAT( 1)+1.D0
      STAT( 2)=STAT( 2)+AREA
      IF(DG.LE.STAT( 3)) THEN
      STAT( 3)=DG
      STAT( 4)=DLAT
      STAT( 5)=DLON
      ENDIF
      IF(DG.GE.STAT( 6)) THEN
      STAT( 6)=DG
      STAT( 7)=DLAT
      STAT( 8)=DLON
      ENDIF
      STAT( 9)=STAT( 9)+DG
      STAT(10)=STAT(10)+DG*AREA
      STAT(11)=STAT(11)+DG**2
      STAT(12)=STAT(12)+DG**2*AREA
      IF(SD.LE.STAT(15)) THEN
      STAT(15)=SD
      STAT(16)=DLAT
      STAT(17)=DLON
      ENDIF
      IF(SD.GE.STAT(18)) THEN
      STAT(18)=SD
      STAT(19)=DLAT
      STAT(20)=DLON
      ENDIF
      STAT(21)=STAT(21)+SD**2
      STAT(22)=STAT(22)+SD**2*AREA
C-----------------------------------------------------------------------
      ENDIF
  110 CONTINUE
C-----------------------------------------------------------------------
      IF(I.NE.NROWS) RETURN
      IF(STAT(1).GT.0.D0) THEN
      STAT( 9)=STAT( 9)/STAT( 1)
      STAT(10)=STAT(10)/STAT( 2)
      STAT(11)=SQRT(STAT(11)/STAT( 1))
      STAT(12)=SQRT(STAT(12)/STAT( 2))
      STAT(13)=SQRT(STAT(11)**2-STAT( 9)**2)
      STAT(14)=SQRT(STAT(12)**2-STAT(10)**2)
      STAT(21)=SQRT(STAT(21)/STAT( 1))
      STAT(22)=SQRT(STAT(22)/STAT( 2))
      STAT( 2)=STAT( 2)/FOURPI*100.D0
      ELSE
      DO 120 J=3,22
      STAT(J)=DEXCLUD
  120 CONTINUE
      ENDIF
C=======================================================================
      NUM=INT(STAT(1))
      WRITE(6,6001) DLABEL(1),NUM
 6001 FORMAT(/5X,A20,3X,I11)
      DO 210 K=2,14
      WRITE(6,6002) DLABEL(K),STAT(K)
 6002 FORMAT(5X,A20,3X,F15.3)
  210 CONTINUE
      WRITE(6,6003)
 6003 FORMAT(' ')
      IF(ISIG.EQ.1) THEN
      DO 220 K=1,8
      WRITE(6,6002) SLABEL(K),STAT(K+14)
  220 CONTINUE
      ENDIF
C=======================================================================
      RETURN
      END

 

[JamesONeil]

Hello, so I created a file containing your source and renamed it .f90 because I had some problems when compiling it that may come from a "fixed format". Anyway I tried to correct what I could but I don't understand these errors :

gfortran read_3files.f90 -o readmin1
read_3files.f90:6:0:

    6 | c-----------------------------------------------------------------------
      |
Error: Unclassifiable statement at (1)
read_3files.f90:88:64:

   88 |       DATA DLABEL/'    Number of Values','  Percentage of Area',
      |                                                                1
Error: Syntax error in DATA statement at (1)
read_3files.f90:89:6:

   89 |      $            '       Minimum Value',' Latitude of Minimum',
      |      1
Error: Invalid character in name at (1)
read_3files.f90:90:6:

   90 |      $            'Longitude of Minimum','       Maximum Value',
      |      1
Error: Invalid character in name at (1)
read_3files.f90:91:6:

   91 |      $            ' Latitude of Maximum','Longitude of Maximum',
      |      1
Error: Invalid character in name at (1)
read_3files.f90:92:6:

   92 |      $            '     Arithmetic Mean','  Area-Weighted Mean',
      |      1
Error: Invalid character in name at (1)
read_3files.f90:93:6:

   93 |      $            '      Arithmetic RMS','   Area-Weighted RMS',
      |      1
Error: Invalid character in name at (1)
read_3files.f90:94:6:

   94 |      $            '   Arithmetic S.Dev.','Area-Weighted S.Dev.'/
      |      1
Error: Invalid character in name at (1)
read_3files.f90:95:64:

   95 |       DATA SLABEL/'       Minimum Sigma',' Latitude of Minimum',
      |                                                                1
Error: Syntax error in DATA statement at (1)
read_3files.f90:96:6:

   96 |      $            'Longitude of Minimum','       Maximum Sigma',
      |      1
Error: Invalid character in name at (1)
read_3files.f90:97:6:

   97 |      $            ' Latitude of Maximum','Longitude of Maximum',
      |      1
Error: Invalid character in name at (1)
read_3files.f90:98:6:

   98 |      $            'Arithmetic RMS Sigma','Area-wghtd RMS Sigma'/
      |      1
Error: Invalid character in name at (1)
read_3files.f90:174:3:

  174 | C=======================================================================
      |   1
Error: Invalid character in name at (1)
read_3files.f90:189:3:

  189 | C=======================================================================
      |   1
Error: Invalid character in name at (1)

Thank you for your help

 

[mikrom]

It's not fortran 90 so you should not use the extension .f90
It's fortran 77, so if you are using gfortran compiler then use the extension .f like me - I used the file name read_3files.f
I compiled it with the command

$ gfortran read_3files.f -o readmin1

 

[JamesONeil]

Ok so I created a new file .f and did exactly this command :

gfortran read_3files.f -o readmin1
read_3files.f:1:1:

    1 | program readmin1
      | 1
Error: Non-numeric character in statement label at (1)

I am using notepad++ to visualize the file so I don't know if there could be a problem with indentation or line but it doesn't seem that it is the problem...

 

[mikrom]

As I posted the source the formatter indented the keyword program on the first column. But it must be in the 7-th column like the keyword implicit. Look again at the source I posted above, now I inserted a comment line on the first line and it should be now ok.

 

[JamesONeil]

Ok great thanks a lot @mikrom now it compiles !
@Space_Fede do you still have the 3 files :
Dg01_cnt2.5x2.5_EGM08_to2190_WGS84_ell_nh
xi_cnt2.5x2.5_EGM08_to2190_WGS84_ell_nh
eta_cnt2.5x2.5_EGM08_to2190_WGS84_ell_nh
because I thought I had them but turns out not and because the link doesn't work anymore and it is not available on NGA anymore it is kind of a pity ^^' But I really need them for my research though...

 

[Space_Fede]

Hi James!
I have them! They are already converted in Little Endian. Unfortunately the files are too large to be attached in the post.
Here it is a we transfer link:

https://we.tl/t-dQTFHbHW6x

 

[JamesONeil]

Great ! Thanks @Space_Fede ! So now I can compile and get the results My LAST question and after I will stop bothering both of you ^^ So this program is compiling only a few values of the files we have and of course it is our choice to add more. So I guess that the values nrowsg and ncolsg have to be changed since they are used in the loops. In fact maybe only nrowsg because ncolsg is not supposed to be more than it is. My question is to what number format corresponds nrowsg =04320 ?
I just want to know how to modify the code in order for example to get values from a window of latitudes and longitudes. From the code given you get only 19 values with latitude 89.979... and different longitudes starting at 0.0208...

      program readmin1
      implicit real*8(a-h,o-z)
      parameter(exclud=9999.d0,
     $   dlatg=2.5d0/60.d0,dlong=2.5d0/60.d0,nrowsg=04320,ncolsg=08640)
         real*4 data(ncolsg)
      dimension sdat(ncolsg,2),stati(22)
c-----------------------------------------------------------------------
c
c   Read-in Gravity Anomaly 2.5x2.5 minute point value file.
c
      write(6,6001)
 6001 format(///15x,'Statistics of Gravity Anomaly Values  (mGal)',/)
      do i = 1, nrowsg
        read(1) (  data(j),j=1,ncolsg)
        rlat = 90.d0 - (i-1.d0)*dlatg - dlatg*0.5d0
        do j = 1, ncolsg
          clon = (j-1.d0)*dlong + dlong*0.5d0
c
c   Print a few values.
c
          if(i.eq.   1.and.j.lt.20) then  !  top row
            write(6,6101) rlat,clon,data(j)
 6101       format(5x,2f15.10,f15.5)
          endif  !  top row
          sdat(j,1) = data(j)
          sdat(j,2) = 1.d0
        enddo  !  j
        call stats(90.d0,0.d0,i,dlatg,dlong,nrowsg,ncolsg,sdat,exclud,0,
     $             stati)
      enddo  !  i

Thanks a lot for your time

 

[mikrom]

From the code given you get only 19 values with latitude 89.979... and different longitudes starting at 0.0208...

Yes, it's because, they want to print only few values that's when i=1 and j is less than 20, i.e.
1 <= j < 20. i is row index and j is column index, so they print only first 19 columns from first row.

c
c   Print a few values.
c
          [highlight #FCE94F]if(i.eq.1.and.j.lt.20) then  !  top row[/highlight]
            write(6,6101) rlat,clon,data(j)
 6101       format(5x,2f15.10,f15.5)
          endif  !  top row

If you comment out or delete the IF with the ENDIF, then you will get all values.

 

[mikrom]

Here is the minimal program to read all values from the first file
read_data.f

c-----------------------------------------------------------------------
      program read_data
      implicit real*8(a-h,o-z)
      parameter(exclud=9999.d0,
     $   dlatg=2.5d0/60.d0,dlong=2.5d0/60.d0,nrowsg=04320,ncolsg=08640)
         real*4 data(ncolsg)
      dimension sdat(ncolsg,2),stati(22)
c-----------------------------------------------------------------------
c
c   Read-in Gravity Anomaly 2.5x2.5 minute point value file.
c
      write(6,6001)
 6001 format(15x,'Statistics of Gravity Anomaly Values  (mGal)')
      do i = 1, nrowsg
        read(1) (  data(j),j=1,ncolsg)
        rlat = 90.d0 - (i-1.d0)*dlatg - dlatg*0.5d0
        do j = 1, ncolsg
          clon = (j-1.d0)*dlong + dlong*0.5d0
c
c   Print a few values.
c
c         if(i.eq.1.and.j.lt.20) then  !  top row
            write(6,6101) rlat,clon,data(j)
 6101       format(5x,2f15.10,f15.5)
c         endif  !  top row
          sdat(j,1) = data(j)
          sdat(j,2) = 1.d0
        enddo  !  j
c       call stats(90.d0,0.d0,i,dlatg,dlong,nrowsg,ncolsg,sdat,exclud,0,
c    $             stati)
      enddo  !  i
C=======================================================================
      RETURN
      END

Compile it with

gfortran read_data.f -o read_data

and run it redirecting the output into a text file

./read_data > data.txt

Now you can observe all data. But the text file is very huge, it has 1.9 GB and 4320*8640 = 37324800 lines.
I'm looking only at begin and end of the file:

$ head data.txt
               Statistics of Gravity Anomaly Values  (mGal)
       89.9791666667   0.0208333333        6.07722
       89.9791666667   0.0625000000        6.07550
       89.9791666667   0.1041666667        6.07377
       89.9791666667   0.1458333333        6.07204
       89.9791666667   0.1875000000        6.07032
       89.9791666667   0.2291666667        6.06860
       89.9791666667   0.2708333333        6.06687
       89.9791666667   0.3125000000        6.06515
       89.9791666667   0.3541666667        6.06343

$ tail data.txt
      -89.9791666667 359.6041666667      -32.69642
      -89.9791666667 359.6458333333      -32.69638
      -89.9791666667 359.6875000000      -32.69634
      -89.9791666667 359.7291666667      -32.69630
      -89.9791666667 359.7708333333      -32.69625
      -89.9791666667 359.8125000000      -32.69622
      -89.9791666667 359.8541666667      -32.69617
      -89.9791666667 359.8958333333      -32.69613
      -89.9791666667 359.9375000000      -32.69609
      -89.9791666667 359.9791666667      -32.69605

 

[mikrom]

I just want to know how to modify the code in order for example to get values from a window of latitudes and longitudes.

Then modify this condition, from which rows and columns you want to get the data:

c
c   Print a few values.
c
          if([highlight #FCE94F]i.eq.1.and.j.lt.20[/highlight]) then  !  top row
            write(6,6101) rlat,clon,data(j)
 6101       format(5x,2f15.10,f15.5)
          endif  !  top row

 

[JamesONeil]

Hello @mikrom, ok I understand know how the data are read and selected. Thanks a lot for your help and I think I will finally be able to be ok on my own with these files to do what I want ^^ I appreciated that you were not enough about some questions that might sound stupid for somebody who knows how to work with fortran but that are important for total beginners.
Have a nice day !