c MODB to HDF 
c original MODB routine by
c
c ------------------------------------------------
c ---         read the ascii data              ---
c ---                                          ---
c ---   this program is ok for data arrays     ---
c ---   with fixed dimensions (167,63,34)      ---
c ---                                          ---
c --- gher :  r.schoenauen  (4/1/94)           ---
c ------------------------------------------------
c
c HDF conversion by Space Research Software
c
c Example 1 
c
c This FORTRAN program reads a MED2 MODB scalar data file
c and saves it in HDF format
c Dimension scales are used to save latitude, longitude and depth
c Text annotations are added to the HDF file

      program MODB2HDF
      implicit none

c in FORTRAN programs always declare implicit none

c modb

      real*4 a(167,63,34)
      real*4 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)

!     dimensions
      real*4 x_pos(167),y_pos(63),z_pos(34)
	integer*4 dim_index, dim_id, count
      integer sfdimid, sfsdscale, sfsdmname
      real*4 dim_scale(167)
	real*4 alat, along
	integer i, j, k
	real*4 z

! annotations

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

	file_name = 'MODB_ex1.hdf'



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

! read the file 'autmedoo.sal'

      name = 'autmedoo.sal' 
      open (unit=10,file=name,form='formatted')
      read (10,102) imax,jmax,kmax
      read (10,100) valex,xmin,xmax
      call readit(a,imax,jmax,kmax)
      close (10)

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

!/////////////////////////////////////////////////////////////
! read and define dimensions
!/////////////////////////////////////////////////////////////


! the point (i,j,k) is located at (-5.5 + i*0.25 , i=1,167) long 
! and (30.0 + j*0.25 , j=1,63) lat at a depth given by the kth line of
! the file med2.zed 

      do i = 1, imax
	 along = -5.5 + float(i)*0.25 
	 x_pos(i) = along
      enddo

	do j = 1, jmax
	 alat = 30.0 + float(j)*0.25
	 y_pos(j) = alat
      enddo

      open (unit=11,file='MED2.zed',form='formatted')
	do k = 1, kmax
	 read (11,*) z
	 write (*,*) z
	 z_pos(k) = z
	enddo
	close (11)



! 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 = 3
      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, 'Salinity', 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 = 'psu'
	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     For each dimension of the array data set, 
      do dim_index = 1, 3

C      - select the dimension id, 
         dim_id = sfdimid(sds_id, dim_index-1)

C      - set the dimension names and values. 
         if (dim_index .eq. 1) then
         status = sfsdmname(dim_id, 'X_Axis') 
	   count = X_LENGTH
	   do i = 1, X_LENGTH
	    dim_scale(i) = x_pos(i) 
	   enddo
         end if

         if (dim_index .eq. 2) then
         status = sfsdmname(dim_id, 'Y_Axis') 
	   count = Y_LENGTH
  	   do i = 1, Y_LENGTH
	    dim_scale(i) = y_pos(i) 
	   enddo
         end if

         if (dim_index .eq. 3) then
         status = sfsdmname(dim_id, 'Z_Axis')
	   count = Z_LENGTH 
	   do i = 1, Z_LENGTH
	    dim_scale(i) = z_pos(i) 
	   enddo
         end if


C     -  write the dimension scale to the data set 
         status = sfsdscale(dim_id, count, DFNT_FLOAT32, dim_scale)

      enddo ! dimensions


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 = 300
      file_label = 'MODB data. AUTUMN salinity'
	file_desc = 'This HDF file shows MODB data converted to the HDF '/
     %   /'file format. Conversion made by Space Reserach '/
     %   /'and IST/MARETEC	(http://hidrox.ist.utl.pt/maretec). '/
	%   /'This file shows the salinity in autumn. '/
     %   /'MODB URL: http://modb.oce.ulg.ac.be' 
    	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 READIT (A,IMAX,JMAX,KMAX)
C ----------------------------------------

      REAL*4 A(IMAX,JMAX,KMAX)
      INTEGER IMAX,JMAX,KMAX

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

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

20    CONTINUE
      write(*,*) 'k = ', k
10    CONTINUE

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
