#! /bin/sh
#
#    Last change: Mihai Dima 2001
#
if [ $# -lt 3 ]
then
  echo " "
  echo "         *****     SUB     *****      "
  echo " "
  echo "  The program subtracts two files (record by record) "
  echo " "
  echo "     Syntax: sub <file1> <file2> <file3>"
  echo " "
  exit 0
else
 echo " "
 echo "     *****     SUB     *****      "
 echo " "
 echo "file1="$1
 echo "file2="$2
 echo "file3="$3
 cp $1 fort.1
 cp $2 fort.2
fi
cat > sub.f <<EOF
      PROGRAM Psub
      IMPLICIT REAL (A-H,O-Z)
      
      NCHAN=1
      CALL PARLES (NCHAN,NDIM1,NSET)
      NCHAN=2
      CALL PARLES (NCHAN,NDIM2,NSET)
      NDIM=MAX0(NDIM1,NDIM2)
      CALL SUB (NDIM)       
      STOP 
      END
C******************************************************************
      SUBROUTINE PARLES (NCHAN,NDIM,NSET)
      IMPLICIT REAL (A-H,O-Z)
C
      INTEGER I1,I2,I3

      NSET=0

      REWIND NCHAN
      READ(NCHAN,END=10) ITIME,INAME,ILEV,ISIZE
      GOTO 20
   10 PRINT *,'INPUT FILE ON TAPE',NCHAN,' IS EMPTY'
      STOP 'MISTAKE'
   20 REWIND NCHAN

      NDIM=0  
    1 READ (NCHAN,END=100) I1,I2,I3,I4
      NDIM=MAX0(NDIM,I4)
      READ (NCHAN,END=101)
      NSET=NSET+1
      GOTO1
      
  101 PRINT *,' MORE HEADER THEN DATA '     
  100 PRINT 1000,NCHAN,NDIM,NSET
      RETURN    
 1000 FORMAT (1X,' TAPE=',I2,',  NDIM=',I5,',  NSET=',I5,/)
C
      END
C*******************************************************************
      SUBROUTINE SUB (NDIM)
      IMPLICIT REAL (A-H,O-Z)
C
      PARAMETER (NDIM77=20000,UNDEF=0.9E+10, UNDEFE=0.9*UNDEF)
      DIMENSION IDATE(4), JDATE(4)
C IF YOU USE THIS PROGRAM NOT ON A CRAY-2 ACTIVATE THE
C FOLLOWING TWO LINES IN THIS PROGRAM !
C     PARAMETER (NDIM77=20000) 
      DIMENSION X(NDIM77), Y(NDIM77)
C AND DELETE THE FOLLOWING LINE
C      DIMENSION X(NDIM), Y(NDIM)
C      POINTER (NPX,X),(NPY,Y)
C      ALLOCATE (NPX,NDIM)
C      ALLOCATE (NPY,NDIM)
C NOW YOU HAVE DONE ALL CHANGES FOR STANDARD FORTRAN.
      REWIND 1
      REWIND 2
      N = 0
    1 READ(1,END=10) IDATE
      READ(2,END=11) JDATE
      N = N + 1
      IF (N.EQ.1) PRINT 1000,N,IDATE,' FIRST RECORD '
      IF (IDATE(4).NE.JDATE(4)) THEN 
        PRINT *,' ** ERROR WITH FIELD LENGTH **'
        PRINT *,' LEN1 NE LEN2 ',IDATE(4),JDATE(4)
        STOP
      ENDIF
      READ(1) (X(I),I=1,IDATE(4))
      READ(2) (Y(I),I=1,IDATE(4))
      WRITE(3) IDATE
      DO 100 I=1,IDATE(4)
        IF (X(I).LT.UNDEFE.AND.Y(I).LT.UNDEFE) THEN
        X(I)=X(I)-Y(I)
        ELSE
        X(I)=UNDEF
        ENDIF
  100 CONTINUE   
      WRITE (3) (X(I), I=1,IDATE(4))
      GOTO 1
  10  PRINT 1000, N,IDATE,' LAST RECORD '
      PRINT *, ' TOTAL NUMBER OF RECORDS FOUND: ', N
      STOP 'O.K.'
  11  PRINT *, ' TOO FEW RECORDS ON **TAPE2** '
      STOP 'MISTAKE'
 1000 FORMAT(2I10,3I6,A20)
      END
EOF
f77 sub.f -o sub.x
sub.x $1 $2 $3
cp fort.3 $3
rm fort.*
rm sub.x sub.f
exit