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 2 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 c c This example reads two files containing the horizontal and vertical c components of the velocity and saves the two data sets in one HDF file program MODB2HDF implicit none c in FORTRAN programs always declare implicit none c modb real*4 u(167,63,34) real*4 v(167,63,34) real*4 valex,xmin,xmax integer imax,jmax,kmax character*50 name real*4 aa 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 real*4 uu, vv ! annotations integer max_char character*50 file_name character*50 file_label character*500 file_desc file_name = 'MODB_ex2_uv.hdf' write (6,'(a)') ' read the ascii data' write (6,'(a)') ' -------------------' write (6,'(a)') ! read u velocity name = 'summedoo.u' open (unit=10,file=name,form='formatted') read (10,102) imax,jmax,kmax read (10,100) valex,xmin,xmax call readit(u,imax,jmax,kmax) close (10) ! read v velocity name = 'summedoo.v' open (unit=10,file=name,form='formatted') read (10,102) imax,jmax,kmax read (10,100) valex,xmin,xmax call readit(v,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) ! /////// check it do i = 1, imax do j = 1, jmax do k = 1, kmax uu = u (i,j,k) vv = v (i,j,k) if ( uu.gt.valex. and. vv.eq.valex) then write (*,*) 'redefinition at ', i, j, k v (i,j,k) = 0.0 u (i,j,k) = 0.0 endif if ( vv.gt.valex. and. uu.eq.valex) then write (*,*) 'redefinition at ', i, j, k u (i,j,k) = 0.0 v (i,j,k) = 0.0 endif enddo enddo enddo !///////////////////////////////////////////////////////////// ! convert to cm/s !///////////////////////////////////////////////////////////// do i = 1, imax do j = 1, jmax do k = 1, kmax uu = u (i,j,k) vv = v (i,j,k) if ( uu . gt. valex ) u (i,j,k) = uu * 100.0 if ( vv . gt. valex ) v (i,j,k) = vv * 100.0 enddo enddo enddo !///////////////////////////////////////////////////////////// ! redefine xmin and xmax for u !///////////////////////////////////////////////////////////// xmax = -1e10 xmin = 1e10 do i=1,imax do j=1,jmax do k=1,kmax aa = u(i,j,k) if ( aa .gt. valex ) then if ( aa .gt. xmax ) xmax = aa if ( aa .lt. xmin ) xmin = aa endif enddo enddo enddo ! 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, 'U velocity', & 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 ! status = sfsdtstr(sds_id, label, units, display, coord) units = 'cm/s' coord = 'Cartesian' status = sfsdtstr(sds_id, '', units, '', coord) C Write the stored data to the data set. status = sfwdata(sds_id, start, stride, edges, u ) 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) !///////////////////////////////////////////////////////////// ! redefine xmin and xmax for v !///////////////////////////////////////////////////////////// xmax = -1e10 xmin = 1e10 do i=1,imax do j=1,jmax do k=1,kmax aa = v(i,j,k) if ( aa .gt. valex ) then if ( aa .gt. xmax ) xmax = aa if ( aa .lt. xmin ) xmin = aa endif enddo enddo enddo C Create another data set. sds_id = sfcreate(sd_id, 'V velocity', & 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 Write the stored data to the data set. status = sfwdata(sds_id, start, stride, edges, v ) 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.' file_desc = 'This HDF file shows MODB data converted to the HDF '/ % /'file format. Conversion made by Space Research '/ % /'and IST/MARETEC (http://hidrox.ist.utl.pt/maretec). '/ % /'This file shows geostrofic velocities in summer. '/ % /'Units are cm/s.'/ % /'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