#include "options.inc"

c
c ------------------------------------------------------
c   template to read topography from external data file
c   which is currently "data/etopo5/etopo5.cdf"
c   containing a 5 minute global topgraphy by Smith et al
c   or data/etopo2/etopo2.top containing a 2 min. version
c
c   Notes: Called by prep_topo (in prep_topo.F)
c          defines two subroutines: 
c           "read_rose_dims" and "read_rose"
c          rose should be given in meters
c          Header of subroutines should not be touched.
c          CPP directives are set in prep_template.F
c
c                            c.eden  nov 2002
c ------------------------------------------------------
c

#ifdef use_ETOPO5_topography

      subroutine read_rose(rose,x,y,nx,ny)
      implicit none
#include "netcdf.inc"
      integer :: nx,ny
      real rose(nx,ny),x(nx),y(ny)
      character*80 infile,name
      character*80 dname
      integer  ncid,iret,rid,ndims,latid,lonid
      integer  corner(4), edges(4), dims(4)
#ifdef netcdf_real4
      real (kind=4) , allocatable :: var4(:,:) ,x4(:),y4(:)
#else
      real , allocatable :: var4(:,:) ,x4(:),y4(:)
#endif

      infile='data/etopo5/etopo5.cdf'
      name='rose'

      iret=nf_open(infile,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,name,rid)
      iret=NF_INQ_VARNDIMS(NCID, rid, ndims)  
      if (ndims/=2) then
       print*,' number of dimensions for variable not equal 2'
       print*,ndims
       stop
      endif
      iret=NF_INQ_VARDIMID (NCID,rid,DIMS)
      iret=NF_INQ_DIMNAME (NCID, DIMs(1),dname) 
      iret=NF_INQ_VARID (NCID, dNAME, lonid)
      iret=NF_INQ_DIMNAME (NCID, DIMs(2),dname) 
      iret=NF_INQ_VARID (NCID, dNAME, latid)
      corner=1; edges(1)=nx; edges(2)=ny
      allocate( var4( nx,ny), x4(nx), y4(ny) )
      iret= nf_get_vara_real (ncid,rid ,corner,edges,var4)
      rose=var4
      edges(1)=ny
      iret= nf_get_vara_real (ncid,latid ,corner,edges,y4); y=y4
      edges(1)=nx 
      iret= nf_get_vara_real (ncid,lonid,corner,edges,x4);x=x4
      deallocate(var4,y4,x4)
      call ncclos (ncid, iret)
      end subroutine read_rose




      subroutine read_rose_dims(nx,ny)
      implicit none
#include "netcdf.inc"
      integer nx,ny
      character*80 infile,name
      character*80 dname
      integer  ncid,iret,rid,ndims,latid,lonid
      integer dims(4)

      infile='data/etopo5/etopo5.cdf'
      name='rose'
      print*,' reading topography from file ',
     &      infile(1:len_trim(infile))

      iret=nf_open(infile,NF_nowrite,ncid)
      if (iret/=0) then; print*,nf_strerror(iret);stop;endif
      iret=nf_inq_varid(ncid,name,rid)
      iret=NF_INQ_VARNDIMS(NCID, rid, ndims)  
      if (ndims/=2) then
       print*,' number of dimensions for variable not equal 2'
       print*,ndims
       stop
      endif
      iret=NF_INQ_VARDIMID (NCID,rid,DIMS)
      iret=NF_INQ_DIMNAME (NCID, DIMs(1),dname) 
      iret=NF_INQ_VARID (NCID, dNAME, lonid)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), ndims); nx=ndims
      iret=NF_INQ_DIMNAME (NCID, DIMs(2),dname) 
      iret=NF_INQ_VARID (NCID, dNAME, latid)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ndims); ny=ndims
      call ncclos (ncid, iret)
      end subroutine read_rose_dims

#endif

#ifdef use_ETOPO2_topography

c
c   read data from unformatted fortran data file
c
      subroutine read_rose(rose,x,y,nx,ny)
      implicit none
      integer :: nx,ny
      real rose(nx,ny),x(nx),y(ny)
      character*80 infile
      real (kind=4), allocatable :: rose4(:,:) ,x4(:),y4(:)
      integer :: io,j
      infile='data/etopo2/etopo2.top'
      allocate( rose4(nx,ny) ,x4(nx),y4(ny) )
      call getunit(io, infile, 'u s r ieee')
      read (io)nx,ny
      read (io)x4; x=x4
      read (io)y4; y=y4
      do j=1,ny
        read(io) rose4(:,j)
      enddo
      close(io)
      rose=rose4*-1
      deallocate(rose4,x4,y4)
      end subroutine read_rose


      subroutine read_rose_dims(nx,ny)
      implicit none
      integer nx,ny
      character*80 infile
      integer :: io
      infile='data/etopo2/etopo2.top'
      print*,' reading topography from file ',
     &      infile(1:len_trim(infile))
      call getunit(io, infile, 'u s r ieee')
      read (io)nx,ny
      close(io)
      end subroutine read_rose_dims


#endif
