C Copyright 1981-2007 ECMWF
C 
C Licensed under the GNU Lesser General Public License which
C incorporates the terms and conditions of version 3 of the GNU
C General Public License.
C See LICENSE and gpl-3.0.txt for details.
C

       PROGRAM SAMPLE1
C
      IMPLICIT NONE
      INTEGER IPROD
      INTEGER INTV
      REAL REALV
      CHARACTER*20 CHARV
      DIMENSION INTV(4), REALV(4), CHARV(4)
C
      INTEGER JPGRIB, JPBYTES
C
      PARAMETER (JPGRIB = 2000000)
C
C     JPBYTES is the size in bytes on an 'INTEGER'
C     Set JPBYTES = 8 on a 64-bit machine.
C
      PARAMETER (JPBYTES = 4)
C
      INTEGER INGRIB, NEWFLD
      DIMENSION INGRIB(JPGRIB), NEWFLD(JPGRIB)
C
      REAL ZNFELDI, ZNFELDO
      DIMENSION ZNFELDI(1), ZNFELDO(1)
C
      INTEGER IUNIT1, IUNIT2, IREC, INLEN, NEWLEN, IRET, NARGS
      INTEGER*4 J
C
C     Externals
      INTEGER INTOUT, INTF, IARGC

      CHARACTER*128 INFILE, OUTFILE, CARG(4)
C
C **********************************************************************
C
C     Pick up file names from command line.
C
      NARGS = IARGC()
      IF( NARGS.LT.4 ) THEN
        print*,'Usage: interpolation_example -i inputfile -o outputfile'
        STOP
      END IF

      DO 101 J=1,NARGS
      CALL GETARG(J,CARG(J))
 101  CONTINUE

      DO 102 J=1,NARGS,2
        IF(CARG(J).EQ.'-i') THEN
           INFILE=CARG(J+1)
        ELSEIF(CARG(J).EQ.'-o') THEN
           OUTFILE=CARG(J+1)
        ELSE
        print*,'Usage: interpolation_example -i inputfile -o outputfile'
        STOP
        END IF
 102  CONTINUE

C     Define the packing accuracy for the new field(s).
C
      INTV(1) = 24
      IRET = INTOUT('accuracy', INTV, REALV, CHARV)
      IF ( IRET.NE.0 ) THEN
        WRITE(*,*) ' First INTOUT failed.'
        STOP
      ENDIF
C
C     Define the geographical area for the new field(s).
C
      REALV(1) =  60.0
      REALV(2) = -10.0
      REALV(3) =  40.0
      REALV(4) =  15.0
      IRET = INTOUT('area', INTV, REALV, CHARV)
      IF ( IRET.NE.0 ) THEN
        WRITE(*,*) ' Second INTOUT failed.'
        STOP
      ENDIF
C
C     Define the grid interval for the new field(s).
C
      REALV(1) = 1.5
      REALV(2) = 1.5
      IRET = INTOUT('grid', INTV, REALV, CHARV)
      IF ( IRET.NE.0 ) THEN
        WRITE(*,*) ' Third INTOUT failed.'
        STOP
      ENDIF
C
C     Open input and output files.
C
      CALL PBOPEN(IUNIT1, INFILE, 'r', IRET)
      IF ( IRET.NE.0 ) STOP ' PBOPEN failed'
      CALL PBOPEN(IUNIT2, OUTFILE, 'w', IRET)
      IF ( IRET.NE.0 ) STOP ' PBOPEN failed'
      IPROD = 0
C
C     Start of loop through input GRIB-coded fields
C
 200  CONTINUE
        IPROD = IPROD + 1
C
C       Read next product.
C
        CALL PBGRIB(IUNIT1, INGRIB, JPGRIB*JPBYTES, IREC, IRET)
        IF ( IRET.EQ.-1 ) GOTO 900
        IF ( IRET.NE.0 ) STOP ' PBGRIB failed'
C
C       Interpolate.
C
        WRITE(*,*) ' Interpolate product number ', IPROD
        NEWLEN = JPGRIB
        INLEN  = IREC/JPBYTES
        IRET = INTF(INGRIB,INLEN,ZNFELDI,NEWFLD,NEWLEN,ZNFELDO)
        IF ( IRET.NE.0 ) THEN
          WRITE(*,*) ' INTF failed.'
          STOP
        ENDIF
C
C       Write the new product to file.
C
       CALL PBWRITE( IUNIT2, NEWFLD, NEWLEN*JPBYTES, IRET)
        IF ( IRET.LT.(NEWLEN*JPBYTES) ) STOP ' PBWRITE failed'
C
C       Loop back for next product.
C
      GOTO 200
C
C     Closedown.
C
 900  CONTINUE
C
      IPROD = IPROD - 1
      WRITE(*,*) ' All done after ', IPROD, ' products.'
C
C     Close input and output files.
C 
      CALL PBCLOSE(IUNIT1, IRET)
      CALL PBCLOSE(IUNIT2, IRET)
C
      STOP
      END
