c ETOPO5 to HDF 
c original MODB routine by
c
C ------------------------------------------------
C ---         Read the ASCII Data              ---
C ---                                          ---
C --- GHER :  R.Schoenauen  (4/1/94)           ---
C ------------------------------------------------
c
c HDF conversion by Space Research Software
c
c
c
c This FORTRAN program reads a MODB gridded data representing bathymetries and elevations 
c extracted from the ETOPO5 data base
c and saves it in HDF format

c The water points (negative values ) are converted to the fill value

      program MODB2HDF
      implicit none

c in FORTRAN programs always declare implicit none

c modb
      integer MAX_DIM
      parameter (MAX_DIM = 800000)

      REAL*4 A(MAX_DIM),VALEX,XMIN,XMAX
      INTEGER IMAX,JMAX,KMAX
      CHARACTER*50 NAME
     
c hdf

      integer*4 sd_id, sds_id, rank
      integer dims(3), start(3), edges(3), stride(3), status
      integer sfstart, sfcreate, sfwdata, sfendacc, sfend
	integer sfsfill, sfsrange
	integer sfsdtstr
	character*20 units,coord


C     DFACC_CREATE and DFNT_FLOAT32 are defined in hdf.h.
      integer*4 DFACC_CREATE, DFNT_FLOAT32
      integer*4 X_LENGTH, Y_LENGTH, Z_LENGTH
      parameter (DFACC_CREATE = 4, DFNT_FLOAT32 = 5)


! annotations

      integer max_char
	character*50 file_name
	character*50 file_label
	character*500 file_desc

	file_name = 'ETOPO5.hdf'

      write (6,'(a)') ' read the ascii data'
      write (6,'(a)') ' -------------------'
      write (6,'(a)')

! read the file


	     name = 'extract2-etopoe10abbd6.txt'
      open (unit=10,file=name,form='formatted')
      read (10,102) imax,jmax,kmax
      read (10,100) valex,xmin,xmax

	if ( imax*jmax .gt. MAX_DIM ) then
	 write (*,*) 'Array too large. Change the MAX_DIM value to ',
     &	 imax*jmax
	 stop
	endif
      call READIT3(a,imax,jmax,kmax,valex,xmin,xmax)
      close (10)

100   format (3e12.5)
101   format (5e12.5)
102   format (3i5)




! write it to hdf

      write (*,*) 'Writing file...'
      
      X_LENGTH = imax
      Y_LENGTH = jmax
	Z_LENGTH = kmax

C     Create and open the file and initiate the SD interface.
      sd_id = sfstart(file_name, DFACC_CREATE)

C     Define the rank and dimensions of the data set to be created.
      rank = 2
      dims(1) = X_LENGTH
      dims(2) = Y_LENGTH
	dims(3) = Z_LENGTH

C     Define the location, pattern, and size of the data set
C     that will be written to.
      start(1) = 0
      start(2) = 0
	start(3) = 0
      edges(1) = X_LENGTH
      edges(2) = Y_LENGTH
	edges(3) = Z_LENGTH
      stride(1) = 1
      stride(2) = 1
	stride(3) = 1

C     Create the data set.
      sds_id = sfcreate(sd_id, 'Elevation', DFNT_FLOAT32, rank, dims)
		
C     Set a fill value.
      status = sfsfill(sds_id, valex )

C     Set range values
      status = sfsrange(sds_id, xmax, xmin )

C     Set attributes
      units = 'm'
	coord = 'Cartesian'
      status = sfsdtstr(sds_id, '', units, '', coord)

C     Write the stored data to the data set.
      status = sfwdata(sds_id, start, stride, edges, a )

C     Terminate access to the array.
      status = sfendacc(sds_id)

C     Terminate access to the SD interface and close the file.
      status = sfend(sd_id)



! Save annotations
      max_char = 400
      file_label = 'ETOPO5 data base'
	file_desc='This HDF file shows MODB (http://modb.oce.ulg.ac.be) '/
     %   /'data converted to the HDF '/
     %   /'file format. Conversion made by Space Research '/
     %   /'and IST/MARETEC	(http://hidrox.ist.utl.pt/maretec). '/
	%   /'This file contains gridded data representing '/
     %   /'bathymetries and elevations extracted from the '/
     %   /'ETOPO5 data base. '/
	%   /'The original water values were converted to the fill value '
    	call save_an(file_name,file_label,file_desc,max_char)



      end


!The first line of the ASCII file contains the dimension of the array : 
! IMAX, JMAX, KMAX. The integers are stored with the FORTRAN format I5. 
!The second line contains the exclusion value, the minimum and the maximum value of this 
! file. These numbers are stored with the FORTRAN format E12.5. 
! The remaining lines contains the data of the array, with 5 numbers per line 
! (except the last line for each I-line). 
! The array is stored in horizontal slices from sea surface to sea bottom and from 
! north to south. So the indexes go from : 
!
!   DO K = KMAX to 1
!       DO J = JMAX to 1
!           DO I = 1 to IMAX
!              read
!           OD
!       OD
!   OD
! 
!              ____________________________
!             /                           /| (imax,jmax,kmax)
!            /        sea surface        / |
!           /                           /  |
!          /__________________________ /   |
!          |                          |    |
!          |                          |    | (imax,jmax,1)        n
!          |                          |   /                      /
!          |                          |  /                      /
!     ^  j |                          | /             w  <-----o-----> e
!  k  |  / |__________________________|/                      /
!     | /                           (imax,1,1)               /
!     |---------->                                          s
!	    i
!
!
!
!



C ----------------------------------------
      SUBROUTINE READIT3 (A,IMAX,JMAX,KMAX,valex,xmin,xmax)
C ----------------------------------------

      REAL*4 A(IMAX,JMAX,KMAX), valex, xmin, xmax, aa
      INTEGER IMAX,JMAX,KMAX
	integer i, j, k


      DO 10 K = KMAX,1,-1
      DO 20 J = JMAX,1,-1

         READ (10,101) (A(I,J,K),I=1,IMAX)

20    CONTINUE
10    CONTINUE


!/////////////////////////////////////////////////////////////
! convert land values to the fill value
!/////////////////////////////////////////////////////////////

      do i = 1, imax
	do j = 1, jmax
	do k = 1, kmax

	 if ( a (i,j,k) .lt.0.0 ) a (i,j,k) = valex
		
	enddo
	enddo
	enddo

!/////////////////////////////////////////////////////////////
! redefine xmin and xmax
!/////////////////////////////////////////////////////////////

      xmax = -1e10
	xmin = 1e10
      do i=1,imax
	do j=1,jmax
	do k=1,kmax
	aa = a(i,j,k)
	if ( aa .lt. valex ) then

	  if ( aa .gt. xmax ) xmax = aa
	  if ( aa .lt. xmin ) xmin = aa

	endif
	enddo
	enddo
	enddo


101   FORMAT (5E12.5)

      END




!/////////////////////////////////////////////////////////////////////
! save_an
! saves annotations 
!/////////////////////////////////////////////////////////////////////

	subroutine save_an(file_name,file_label,file_desc,max_char)

	implicit none

c in FORTRAN programs always declare implicit none

      character*(*) file_name 
	character*(*) file_label 
      character*(*) file_desc  
	integer max_char
   
      integer hopen, hclose, afstart, affcreate, afendaccess
      integer status, file_id, an_id, ann_id
      integer afwriteann

      integer*4 DFACC_WRITE, AN_FILE_LABEL
      integer*4 AN_FILE_DESC
      parameter (DFACC_WRITE = 2, AN_FILE_LABEL = 2)
      parameter (AN_FILE_DESC = 3 )

C     Open the HDF file.
      file_id = hopen( file_name, DFACC_WRITE, 0)

C     Initialize the AN interface and obtain an interface id.
      an_id = afstart(file_id)

! label

C     Create the file label and obtain an annotation id
      ann_id = affcreate(an_id, AN_FILE_LABEL)

C     Write the file label to the file.
      status = afwriteann(ann_id, file_label, max_char)

C     Terminate access to the annotation.
      status = afendaccess(ann_id)

! description

C     Create the file label and obtain an annotation id
      ann_id = affcreate(an_id, AN_FILE_DESC)

C     Write the file label to the file.
      status = afwriteann(ann_id, file_desc, max_char)

C     Terminate access to the annotation.
      status = afendaccess(ann_id)

	
C     Close the file.
      status = hclose(file_id)

      end