Layout of routine rd_coor
SUBROUTINE RD_COOR(NCID)
#ifdef DOC
!**** *RDCOOR* - Reading netCDF file containing point coordinates
! Purpose.
! --------
! initialization surface coordinates
! writing coordinates to output files
! reading the catchment mask
!** Interface.
! ----------
! *CALL* *RDCOOR*
! Explicit arguments :
! --------------------
! INPUT
! NCID integer NETCDF FILE ID
! Implicit arguments :
! --------------------
! Method.
! -------
! Opens a file called 'surfclim' to read relevant fields
!
! Externals.
! ----------
! NETCDF-utilities
! Reference.
! ----------
! Author.
! -------
! Bart vd Hurk, KNMI
! Modifications.
! --------------
! Original : 2000-07-13
! ------------------------------------------------------------------
#endif
......
#include "netcdf.inc"
!* check the dimensions
NIDLAT = NCDID(NCID, 'y', IERR)
NIDLON = NCDID(NCID, 'x', IERR)
CALL NCDINQ(NCID,NIDLAT,CNAME,NILAT,IERR)
CALL NCDINQ(NCID,NIDLON,CNAME,NILON,IERR)
IF(NILON.NE.NLON .OR. NILAT.NE.NLAT) THEN
.....
ENDIF
!* read coordinates
ALLOCATE (ZLAT(NLAT))
ALLOCATE (ZLON(NLON))
ALLOCATE (DZLAT(NLAT))
ALLOCATE (DZLON(NLON))
NVARID = NCVID(NCID, 'lat', IERR)
CALL NCVINQ (NCID, NVARID, CNAME, NVARTYP,
+ NVDIMS, NDIMS,NVATTS, IERR)
IF(NVARTYP.EQ.NCFLOAT)THEN
CALL NCVGT(NCID, NVARID, 1, NLAT, ZLAT, IERR)
DZLAT(:) = ZLAT(:)
ELSE
CALL NCVGT(NCID, NVARID, 1, NLAT, DZLAT, IERR)
ENDIF
NVARID = NCVID(NCID, 'lon', IERR)
CALL NCVINQ (NCID, NVARID, CNAME, NVARTYP,
+ NVDIMS, NDIMS,NVATTS, IERR)
IF(NVARTYP.EQ.NCFLOAT)THEN
CALL NCVGT(NCID, NVARID, 1, NLON, ZLON, IERR)
DZLON(:) = ZLON(:)
ELSE
CALL NCVGT(NCID, NVARID, 1, NLON, DZLON, IERR)
ENDIF
!* -- catchment mask
ALLOCATE (ZREAL(NLALO))
ALLOCATE (DZREAL(NLALO))
ALLOCATE (LMASK(NLALO))
ISTART2(1)=1
ISTART2(2)=1
ICOUNT2(1)=NLON
ICOUNT2(2)=NLAT
NVARID = NCVID(NCID, 'Mask', IERR)
CALL NCVINQ (NCID, NVARID, CNAME, NVARTYP,
+ NVDIMS, NDIMS,NVATTS, IERR)
IF(NVARTYP.EQ.NCFLOAT)THEN
CALL NCVGT(NCID, NVARID, ISTART2, ICOUNT2, ZREAL, IERR)
LMASK(:)=(ABS(ZREAL(:)-1.).LT.0.0001)
ELSE
CALL NCVGT(NCID, NVARID, ISTART2, ICOUNT2, DZREAL, IERR)
LMASK(:)=(ABS(DZREAL(:)-1.).LT.0.0001)
ENDIF
NPOI=COUNT(LMASK)
WRITE(NULOUT,*)'TOTAL NUMBER OF GRID POINTS: ',NLALO
WRITE(NULOUT,*)'NUMBER OF GRID POINTS IN CATCHMENT: ',NPOI
!* -- convert to radians
ALLOCATE (ZGELAT(NLALO))
ALLOCATE (ZGELAM(NLALO))
DO JLA=1,NLAT
DO JLO=1,NLON
ICT=ICT+1
ZGELAT(ICT)=...
ZGELAM(ICT)=...
ENDDO
ENDDO
!* -- Pack the arrays (Fortran90 utility defining a new array that
!* contains a selection of the original array only)
ALLOCATE (GELAT(NPOI))
ALLOCATE (GELAM(NPOI))
GELAT=PACK(ZGELAT,LMASK)
GELAM=PACK(ZGELAM,LMASK)
!* -- write coordinates to output files
IPOS(1)=NPOSGG
IPOS(2)=NPOSEFL
IPOS(3)=NPOSWAT
IPOS(4)=NPOSRES
IPOS(5)=NPOSSUS
IPOS(6)=NPOSSUB
IPOS(7)=NPOSEVA
IPOS(8)=NPOSCLD
DO J=1,JPNCDF
NPOS = IPOS(J)
NVARID = NCVID(NPOS, 'lat', IERR)
ISTART1(1)=1
ICOUNT1(1)=NLAT
CALL NCVPT (NPOS, NVARID, ISTART1, ICOUNT1, DZLAT, IERR)
NVARID = NCVID(NPOS, 'lon', IERR)
ISTART1(1)=1
ICOUNT1(1)=NLON
CALL NCVPT (NPOS, NVARID, ISTART1, ICOUNT1, DZLON, IERR)
ENDDO
END SUBROUTINE RD_COOR