program readvh1data

! read in 2D+time netcdf dataset
! compile with:
! ifort -o readnc readnc.f90 -I/ncsu/netcdftools/include -L/ncsu/netcdftools/lib -lnetcdf
      
use netcdf !netcdf modules
      
IMPLICIT NONE

CHARACTER(LEN=8) :: filename
CHARACTER(LEN=16) :: coord1, coord2, coord3, coord4
CHARACTER(LEN=16), DIMENSION(:), ALLOCATABLE :: varname


INTEGER :: ncid, ncstat

INTEGER :: n, i, j, ndim, nvar, nv, nslice
INTEGER :: imax, jmax, nmax
INTEGER, DIMENSION(4) :: data_start, data_size

REAL, DIMENSION(:,:,:), ALLOCATABLE :: zro, zpr, zux, zuy
REAL, DIMENSION(:), ALLOCATABLE :: zxc, zyc, time

write(6,*) 'Read a 2D+time netCDF data file.'
write(6,fmt="('Input name of netCDF data file: ')")
read (5,"(a8)") filename

!######################################################################
! OPEN DATA FILE
ncstat = nf90_open(filename, nf90_nowrite, ncid)

! retrieve dimensions and coordinates
 ncstat = nf90_inquire(ncid, ndim, nvar)
 ncstat = nf90_inquire_dimension(ncid, 1, coord1, imax)
 ncstat = nf90_inquire_dimension(ncid, 2, coord2, jmax)
 ncstat = nf90_inquire_dimension(ncid, 3, coord3, nmax)

print *, imax, jmax, nmax, ndim

 nvar = nvar - ndim  ! first ndim variables are the coordinate arrays
 ALLOCATE( varname(nvar) )
 do nv = 1, nvar
   ncstat = nf90_inquire_variable(ncid, nv+ndim, varname(nv))
   print *, nv, varname(nv)
 enddo


 ! allocate coordinate arrays
 ALLOCATE( zxc(imax) )
 ALLOCATE( zyc(jmax) )
 ALLOCATE( time(nmax) )

 ! allocate data arrays
 ALLOCATE(zro(imax,jmax,nmax) )
 ALLOCATE(zpr(imax,jmax,nmax) )
 ALLOCATE(zux(imax,jmax,nmax) )
 ALLOCATE(zuy(imax,jmax,nmax) )


 ! read coordinates arrays
 ncstat = nf90_get_var(ncid, 1, zxc)
 ncstat = nf90_get_var(ncid, 2, zyc)
 ncstat = nf90_get_var(ncid, 3, time)
  

 ! read variable arrays
 ncstat = nf90_get_var(ncid, ndim+1, zro); if (ncstat /= nf90_NoErr ) print *, NF90_STRERROR(ncstat)
 ncstat = nf90_get_var(ncid, ndim+2, zpr)
 ncstat = nf90_get_var(ncid, ndim+3, zux)
 ncstat = nf90_get_var(ncid, ndim+4, zuy)
 
! CLOSE netCDF FILE
ncstat = nf90_close(ncid)
!############################################################################

print *, zro(1,1,1), zpr(1,1,1)

END PROGRAM readvh1data