**************************************************************************
***   readheight.f -- TEMIS -- 18 December 2002
***
***   read elevation data file made from the GTOPO30 data base
**************************************************************************
*
      subroutine readheight(heightfile,mlon,mlat,munit,
     #           height,nlon,nlat,nfac,rlon,rlat,nocean,nland,ierr)
c
c  Subroutine to read the elevation data from the GTOPO30 data base,
c  averaged to a lower resolution.
c
c  INPUT:  heightfile = data file
c          mlon       = max. number of longitudes; array dimension
c          mlat       = max. number of latitudes; array dimension
c          munit      = Fortran unit number used for reading the data file;
c                       should be one of 1,2,3,4,7,...,99
c  OUTPUT: height(mlon,mlat) = array with elevation data; is integer*2 !!
c                              values outside (nlon,nlat) are set to nocean
c          nlon       = actual number of longitudes <= mlon
c          nlat       = actual number of latitudes <= mlat
c          nfac       = resolution specifier: resolution is nfac/120 degrees
c                       in latitude and longitude range  [1]
c          rlon(mlon) = array with longitude values for 1,2,...,nlon  [2]
c          rlat(mlat) = array with latitude values for 1,2,...,nlat   [2]
c          nocean     = value in the array indicating an ocean data point  [3]
c          nland      = part of the nfac*nfac averaged data points that
c                       have to be land for the average data point in the
c                       height array to be defined as land  [4]
c          ierr       = error indicator:
c                       0 = no error occured
c                       1 = invalid unit number munit
c                       2 = data file does not exist
c                       3 = error opening data file
c                       4 = error reading data specifiers
c                       5 = array size of height too small: nlon > mlon
c                       6 = array size of height too small: nlat > mlat
c                       7 = error reading elevation data
c
c  Notes:
c  [1] if nfac=60 then 60 by 60 data points of the original data set are
c      averaged and the resolution is then 60/120=0.5 degrees; the number
c      of data points in the map height is therefore:
c          nlon = 360 * 120/nfac
c          nlat = 180 * 120/nfac
c  [2] these coordinates are at the centre of the grid box, which measures
c      nfac by nfac degrees, and are therefore given by:
c          rlon(i) = -180 + D(i)    i=1,2,...,nlon
c          rlat(i) =  +90 - D(i)    i=1,2,...,nlat
c      where: D(i) = ( 2*i - 1 ) * nfac / 240
c  [3] when averaging data points from the original elevation map, the
c      ocean was given 0 metre height; the values in the data array for
c      an ocean data point is usually -9999
c  [4] if for example nland=1800, then at least 1800/(nfac*nfac) data
c      points in the averaging had to be land for the average data point
c      in the height array to be defined as land, i.e. have a value
c      different from 'nocean'
c
c ------------------
c
      implicit none
c
c  Variables in subroutine argument.
c
      integer mlon,mlat,munit
      integer*2 height(mlon,mlat)
      integer nlon,nlat,nfac,nocean,nland,ierr
      real*8 rlon(mlon),rlat(mlat)
      character heightfile*(*)
c
c  Local variables.
c
      integer ilon,ilat,niostat
      logical lexist
c
c ------------------
c
c  Is the given unit number a valid one?
c
      if ((munit.le.0).or.(munit.ge.100).or.
     #    (munit.eq.5).or.(munit.eq.6)) then
         ierr=1
         return
      endif
c
c  Does the data file exist?
c
      inquire(file=heightfile,exist=lexist)
      if (.not.lexist) then
         ierr=2
         return
      endif
c
c  Can the data file be opened?
c  The data file is assumed to be Fortran unformatted; for reading
c  a formatted file, adapt the 'open' and 'read' lines.
c
      open(UNIT=munit,FILE=heightfile,
     #     FORM="unformatted",STATUS="old",IOSTAT=niostat)
cf    open(UNIT=munit,FILE=heightfile,                 
cf   #     FORM="formatted",STATUS="old",IOSTAT=niostat)
      if (niostat.ne.0) then
         ierr=3
         return
      endif
c
c  Read the data specifiers and check the array size.
c
      read(UNIT=munit,IOSTAT=niostat)nlon,nlat,nfac,nocean,nland
cf    read(UNIT=munit,FMT=*,IOSTAT=niostat)nlon,nlat,nfac,nocean,nland
      if (niostat.ne.0) then
         ierr=4
         close(UNIT=munit)
         return
      endif
c
      if (nlon.gt.mlon) then
         ierr=5
         close(UNIT=munit)
         return
      endif
      if (nlat.gt.mlat) then
         ierr=6
         close(UNIT=munit)
         return
      endif
c
c  Read the elevation data; fill the array elements not in the data
c  file with 'nocean' values, so as to at least have a value there.
c
      do ilat=1,nlat
         read(UNIT=munit,IOSTAT=niostat)
     #        (height(ilon,ilat),ilon=1,nlon)
cf       read(UNIT=munit,FMT=*,IOSTAT=niostat)
cf   #        (height(ilon,ilat),ilon=1,nlon)
         if (niostat.ne.0) then
            ierr=7
            return
         endif
         do ilon=nlon+1,mlon
            height(ilon,ilat)=nocean
         enddo
      enddo
      close(UNIT=munit)
      do ilat=nlat+1,mlat
         do ilon=1,mlon
            height(ilon,ilat)=nocean
         enddo
      enddo
c
c  The longitudes and latitudes of the data points are of the centre
c  of the grid boxes, as written above.
c  Fill the remainder of the arrays with zeros.
c
      do ilon=1,nlon
         rlon(ilon)=-180.d0+dble(2*ilon-1)*dble(nfac)/240.d0
      enddo
      do ilon=nlon+1,mlon
         rlon(ilon)=0.d0
      enddo
c
      do ilat=1,nlat
         rlat(ilat)=90.d0-dble(2*ilat-1)*dble(nfac)/240.d0
      enddo
      do ilat=nlat+1,mlat
         rlat(ilat)=0.d0
      enddo
c
c  No error has occurred reading the elevation data.
c
      ierr=0
c
c ------------------
c
      return
      end
*
**************************************************************************