Layout of routine rd_forc
      SUBROUTINE RD_FORC
#ifdef DOC
!**** *RD_FORC  * - ROUTINE TO SETUP THE ATMOSPHERIC FORCING DATA
!                  FROM NETCDF INPUT
!     PURPOSE.
!     --------
!        Initialize the common block YOMFORC
!**   INTERFACE.
!     ----------
!        *CALL* *RD_FORC*
!     EXPLICIT ARGUMENTS :  
!     --------------------
!     IMPLICIT ARGUMENTS :
!     --------------------
!        COMMON  YOMFORC
!        COMMON  YOMRIP
!        COMMON  YOMLUN1S
!        COMMON  YOMDYN1S
!        COMMON  YOMGC1S
!        COMMON  YOMLOG1S
!     METHOD.
!     -------
!     EXTERNALS.
!     ----------
!        IYMD2C
!     REFERENCE.
!     ----------
!        ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE IFS
!     AUTHOR.
!     -------
!        JEAN-FRANCOIS MAHFOUF AND PEDRO VITERBO  *ECMWF*
!     MODIFICATIONS.
!     --------------
!        ORIGINAL : 95-03-13
!        BART VD HURK (KNMI): READING NETCDF INPUT
#endif
.......
      IMPLICIT LOGICAL (L)
      CHARACTER CHEADER*400
      CHARACTER*100 CNAME
      REAL*4,ALLOCATABLE :: ZTIMAR(:), ZREAL3(:,:)
      REAL,ALLOCATABLE :: ZREAL3D(:,:)
      INTEGER ISTART3(3),ICOUNT3(3)
#include "netcdf.inc"
!---------------------------------------------------------------------
! CONTENTS OF NETCDF FILE
!  lat:units = "degrees_north";long_name = "latitude";
!  lon:units = "degrees_east";long_name = "longitude";
!  time:units = "seconds";long_name = "Seconds since 19790101:00.00";
!  Tair:units = "K";long_name = "Temperature";
!  Qair:units = "kg/kg";long_name = "Specific humidity";
!  Wind_E:units = "m/s";long_name = "Wind speed u";
!  Wind_N:units = "m/s";long_name = "Wind speed v";
!  SWdown:units = "W/m2";long_name = "Downward shortwave radiation";
!  LWdown:units = "W/m2";long_name = "Downward longwave radiation";
!  Rainf:units = "kg/m2s";long_name = "Rainfall";
!  Snowf:units = "kg/m2s";long_name = "Snowfall";
!  PSurf:units = "Pa";long_name = "Pressure";
!---------------------------------------------------------------------
!* Open forcing file
      NCID = NCOPN('forcing', NCNOWRIT, IERR)
      WRITE(NULOUT,*)'NETCDF-FILE forcing OPENED ON UNIT ',NCID
!* check the dimensions
      NIDLAT = NCDID(NCID, 'y', IERR)
      NIDLON = NCDID(NCID, 'x', IERR)
      NIDTIM = NCDID(NCID, 'tstep', IERR)
      CALL NCDINQ(NCID,NIDLAT,CNAME,NILAT,IERR)
      CALL NCDINQ(NCID,NIDLON,CNAME,NILON,IERR)
      CALL NCDINQ(NCID,NIDTIM,CNAME,NITIM,IERR)
      IF(NILON.NE.NLON .OR. NILAT.NE.NLAT) THEN
.....
      ENDIF
      ALLOCATE (ZTIMAR(NITIM))
!     1. READ DATA PROPERTIES
      ZPHISTA=.....
      ZDTFORC=.....
      ZUV=.....
      IFYYYY=.....
      IFMM=.....
      IFDD=.....
      IFTIM=0
      ISTYYYY=0
      ISTMM=0
      ISTDD=0
      ISTTIM=0
      INSTFC=0
!* namelist input
      REWIND(NULNAM)
      READ(NULNAM,NAMFORC)
      DTIMFC=ZDTFORC
      RALT = ZPHISTA
      RZUV = ZUV
!     2. READ TIME DATA
      NVARID = NCVID(NCID, 'time', IERR)
      CALL NCVGT1(NCID, NVARID, 1, ZTIMAR(1), IERR)
      CALL NCVGT1(NCID, NVARID, NITIM, ZTIMAR(NITIM), IERR)
......
!     2.1 calculate first and last time step to read
!
!     Note labeling convention:
!     *** time steps in forcing are labeled by the END TIME of the
!         forcing interval
!     *** required start of read interval is labeled by the START TIME of the
!         first interval
!     *** forcing start time RTSTFC is labeled by the CENTRE of the
!         first forcing time step
      IRDST=....
      NSTPFC=....
!* Check if start of model run > start of forcing
......
!* check if end of model run < end of forcing
......
!     3. read forcing data
      ISTART3(1)=1
      ISTART3(2)=1
      ISTART3(3)=IRDST
      ICOUNT3(1)=NILON
      ICOUNT3(2)=NILAT
      ICOUNT3(3)=NSTPFC
      ALLOCATE (ZREAL3(NLALO,NSTPFC))
      ALLOCATE (ZREAL3D(NLALO,NSTPFC))
!* wind speed u
      NVARID = NCVID(NCID, 'Wind_E', IERR)
      CALL NCVGT(NCID, NVARID, ISTART3, ICOUNT3, ZREAL3, IERR)
      ZREAL3D(:,:)=ZREAL3(:,:)
      DO JT=1,NSTPFC
        UFI(:,JT)=PACK(ZREAL3D(:,JT),LMASK)
      ENDDO
      WRITE(NULOUT,*) " VALUES AT FIRST TIME STEP:"
!* diagnose values at first time step
      CALL MINMAX('WIND U',ZREAL3D(1,1),NLON,NLAT,LMASK,NULOUT)
!* wind speed v
      NVARID = NCVID(NCID, 'Wind_N', IERR)
      CALL NCVGT(NCID, NVARID, ISTART3, ICOUNT3, ZREAL3, IERR)
      ZREAL3D(:,:)=ZREAL3(:,:)
      DO JT=1,NSTPFC
        VFI(:,JT)=PACK(ZREAL3(:,JT),LMASK)
      ENDDO
      CALL MINMAX('WIND V',ZREAL3D(1,1),NLON,NLAT,LMASK,NULOUT)
!* interpolate wind to correct height
      IF(ZUV.NE.RALT)THEN
        DO JL=1,NPOI
          ZFAC=LOG(RALT/VFZ0F(JL))/LOG(ZUV/VFZ0F(JL))
          UFI(JL,1:NSTPFC)=ZFAC*UFI(JL,1:NSTPFC)
          VFI(JL,1:NSTPFC)=ZFAC*VFI(JL,1:NSTPFC)
        ENDDO
        WRITE(NULOUT,*)'SUFCDF: WIND DATA INTERPOLATED FROM ',ZUV,
     *     ' TO ',RALT,' M'
      ENDIF
!* temperature
      NVARID = NCVID(NCID, 'Tair', IERR)
      CALL NCVGT(NCID, NVARID, ISTART3, ICOUNT3, ZREAL3, IERR)
      ZREAL3D(:,:)=ZREAL3(:,:)
      DO JT=1,NSTPFC
        TFI(:,JT)=PACK(ZREAL3D(:,JT),LMASK)
      ENDDO
      CALL MINMAX('TEMP',ZREAL3D(1,1),NLON,NLAT,LMASK,NULOUT)
!* specific humidity
      NVARID = NCVID(NCID, 'Qair', IERR)
.....(etc)
!* surface pressure
      NVARID = NCVID(NCID, 'PSurf', IERR)
.....(etc)
!* shortwave radiation
      NVARID = NCVID(NCID, 'SWdown', IERR)
.....(etc)
!* longwave radiation
      NVARID = NCVID(NCID, 'LWdown', IERR)
.....(etc)
!* rainfall (all large scale)
      NVARID = NCVID(NCID, 'Rainf', IERR)
.....(etc)
!* snowfall (all large scale)
      NVARID = NCVID(NCID, 'Snowf', IERR)
.....(etc)
      WRITE(NULOUT,*) " FORCING DATA READ FOR ",NSTPFC," FORCING STEPS"
      RETURN
      END SUBROUTINE RD_FORC