#! /bin/sh
#
#    Last change: Mihai Dima 2001
#
if [ $# -lt 3 ]
then
  echo " "
  echo "         *****     DETREND     *****      "
  echo " "
  echo " The program fits and removes a linear trend"
  echo " "
  echo "   Syntax: trend <file1> <file2> <file3>"
  echo " "
  echo "         <file1> - input file"
  echo "         <file2> - detrended file"
  echo "         <file3> - coefficients' file"
  echo " "
  exit 0
else
 echo " "
 echo "         *****     DETREND     *****      "
 echo " "
 echo "file1="$1
 echo "file2="$2
 echo "file3="$3
 echo " "
 cp $1 fort.1
fi
cat > detrend.f << EOF
      PROGRAM PROCEC

      PARAMETER(NCHAN=1)
      REAL VALUES
  
      CALL PARLES (NCHAN,NDIM,NSET)
      CALL WORKIN(NDIM,NSET)

      END
C************************************************************
      SUBROUTINE WORKIN (IE,NSET)
      IMPLICIT REAL (A-H,O-Z)
C
C      *TAPE1*  - T21 STANDARD INPUTFILE.
C      *TAPE2*  - T21 STANDARD OUTPUTFILE.
C
      PARAMETER (UNDEF=0.9E+10,UNDEFE=.9*UNDEF)
      LOGICAL MERID
C IF YOU USE THIS PROGRAM NOT ON A CRAY-2 DELETE THE
C FOLLOWING TWO LINES
      PARAMETER(IE77=20000)
      DIMENSION X(IE),SY2(IE),XYM(IE),NCOUNT(IE),XMEAN(IE)
C      INTEGER NCOUNT

      DATA MERID/.FALSE./
      DATA ONE/1.0/

C      PRINT *,'GET INPUT PARAMETER'
C      READ(5, *,END=12) VALUES

C*  CHECK THE INPUT FILE.

      REWIND 1
      READ(1,END=10) ITIME,INAME,ILEV,ISIZE
      IF(ISIZE.GT.IE) GOTO 11
      IE=ISIZE
      REWIND 1
      WRITE(6,6010) ' FIRST RECORD:',ITIME,INAME,ILEV,ISIZE

C
C*  MAIN LOOP.
C

C     MITTELWERT VON TIME = TMEAN
      TMEAN = FLOAT(NSET)/2.0 + 0.5
      DO I=1,IE
        NCOUNT(I)=0
      ENDDO
      ICNT=0
C 1.LOOP CALC MEAN
      DO N=1,NSET
         READ(1,END=15) ITIME,INAME,ILEV,ISIZE
         IF(ISIZE.GT.IE) GOTO 14
         READ(1,END=15) (X(I),I=1,IE)
         ICNT=ICNT+1
         DO I=1, IE
            IF(ICNT.EQ.1) THEN
               XMEAN(I)=0.0         
               IF(X(I).LT.0.98*UNDEF) THEN
                  XMEAN(I)=X(I)         
                  NCOUNT(I)=NCOUNT(I)+1
               ENDIF
            ENDIF
            IF(X(I).LT.0.98*UNDEF) THEN
               XMEAN(I)=XMEAN(I)+X(I)
               NCOUNT(I)=NCOUNT(I)+1
            ENDIF
            IF(N.EQ.NSET) THEN
               IF(NCOUNT(I).GT.0) THEN
                  XMEAN(I)=XMEAN(I)/NCOUNT(I)       
               ELSE
                  XMEAN(I)=UNDEF
               ENDIF
            ENDIF
         ENDDO
      ENDDO
      REWIND 1
C 2.LOOP CALC SX2, SY2, XYM
      XN1=1/FLOAT(NSET-1)
      SX2=0.0
      DO N=1, NSET
         TANO=N-TMEAN
         SX2=SX2+TANO**2
         IF(N.EQ.NSET) SX2=XN1 * SX2
         READ(1,END=15) ITIME,INAME,ILEV,ISIZE
         IF(ISIZE.GT.IE) GOTO 14
         READ(1,END=15) (X(I),I=1,IE)
         DO I=1,IE
            IF(N.EQ.1) THEN
               NCOUNT(I)=0
               SY2(I)=0.0
               XYM(I)=0.0
            ENDIF
            IF(X(I).LT.0.98*UNDEF) THEN
               NCOUNT(I)=NCOUNT(I)+1
               X(I)=X(I)-XMEAN(I)
               XYM(I)=XYM(I)+X(I)*TANO
               SY2(I)=SY2(I)+X(I)**2
            ENDIF
            IF(N.EQ.NSET) THEN
               XN1=1/FLOAT(NCOUNT(I)-1)
               IF(XN1.GT.0) THEN
                  SY2(I)=XN1 * SY2(I)
                  XYM(I)=XN1 * XYM(I)
               ELSE
                  SY2(I)=UNDEF
                  XYM(I)=UNDEF
               ENDIF
            ENDIF
         ENDDO
      ENDDO
      REWIND 1
C 3.LOOP WRITE DETREND VALUES
      DO N=1,NSET
         TREND=N-TMEAN
         READ(1,END=15) ITIME,INAME,ILEV,ISIZE
         IF(ISIZE.GT.IE) GOTO 14
         READ(1,END=15) (X(I),I=1,IE)
         DO I=1,IE
            IF(XYM(I).LT.UNDEF) THEN
               IF(SX2.NE.0) BXY=XYM(I)/SX2
               IF(SX2.EQ.0) BXY=0.0
            ELSE
               BXY=UNDEF
            ENDIF
c            dummy=x(I)
            IF(X(I).LT.0.98*UNDEF) THEN
               X(I)=X(I) - BXY*TREND
C MIST               X(I)=X(I) - BXY*(X(I)-XMEAN(I))
            ELSE
               X(I)=UNDEF
            ENDIF
c            difi=x(I)-dummy
c            if(i.eq.970) print*,'x(I)=',dummy,' xdiff(I)=',difi
         ENDDO
         WRITE(2) ITIME,INAME,ILEV,IE
         WRITE(2) (X(I),I=1,IE)
      ENDDO

      DO I=1, IE
         IF(XYM(I).LT.0.98*UNDEF) THEN
            IF(SX2.NE.0) X(I)=XYM(I)/SX2
            IF(SX2.EQ.0) X(I)=0.0
         ELSE
            X(I)=UNDEF
         ENDIF
      ENDDO
      WRITE(3) 10101,1,0,IE
      WRITE(3) (X(I),I=1,IE)
      DO I=1, IE
         IF(XYM(I).LT.0.98*UNDEF) THEN
            X(I)=UNDEF
            IF(SX2*SY2(I).NE.0) X(I)=XYM(I)/(SQRT(SX2)*SQRT(SY2(I)))
         ELSE
            X(I)=UNDEF
         ENDIF
      ENDDO
      WRITE(3) 10101,2,0,IE
      WRITE(3) (X(I),I=1,IE)

      GOTO 9920

C
C END MAIN LOOP
C

   10 WRITE(6,6010) ' INPUT FILE IS EMPTY '
         GOTO 9910
   11 WRITE(6,6010) ' INPUT GRIDS TOO LARGE '
         GOTO 9910
   12 WRITE(6,6010) ' PARAMETER CARD TOO SHORT OR MISSING '
         GOTO 9910
   13 WRITE(6,6010) ' NRWOS * NCOL .NE. IE '
         GOTO 9910
   14 WRITE(6,6010) ' INVALID RECORD LENGTH:', ISIZE
         GOTO 9910
   15 WRITE(6,6010) ' UNEXPECTED *END OF FILE* DURING READ '
         GOTO 9910
   19 WRITE(6,6010) ' ILAT 1 > ILAT2 ', ILAT1, ILAT2
 9910 CONTINUE
C
C*  ABNORMAL TERMINATION.
C
         PRINT *, '*EXTRA*: FATAL END '
         ONE=SQRT(-ONE)
         STOP
C
C*  NORMAL TERMINATION.
C
 9920 WRITE(6,6010) '  LAST RECORD:',ITIME,INAME,ILEV,ISIZE
      WRITE(6,6010) ' NEW FIELD SIZE:',ISIZE
      WRITE(6,6060) ICNT
      PRINT *, '*EXTRA*: NORMAL '
      STOP
 6010 FORMAT(' *EXTRA*: ', A,/,(6I10))
 6011 FORMAT(' *EXTRA*: ', A)
 6060 FORMAT(' *EXTRA* -- READ ',I10,' RECORDS.')
      END
C*******************************************************************
      SUBROUTINE PARLES (NCHAN,NDIM,NSET)
      IMPLICIT REAL (A-H,O-Z)
C
      INTEGER I1,I2,I3

      REWIND NCHAN
      READ(NCHAN,END=10) ITIME,INAME,ILEV,ISIZE
      REWIND NCHAN

      NSET=0
      NDIM=0  
    1 READ (NCHAN,END=100) I1,I2,I3,I4
      NDIM=MAX0(NDIM,I4)
      READ (NCHAN,END=101) 
      NSET=NSET+1
      GOTO 1
      
 101  PRINT *,' MORE HEADER THEN DATA '     
 100  IF(NSET.LT.3) GOTO 20
      PRINT 1000,NCHAN,NDIM,NSET
      REWIND NCHAN
      RETURN    
 10   PRINT *,'INPUT FILE ON TAPE',NCHAN,' IS EMPTY'
      GOTO 30
 20   PRINT *,'JUST ',NSET,' SETS ARE TOO LESS'
 30   STOP 'MISTAKE'
 1000 FORMAT (1X,' TAPE=',I2,',  NDIM=',I5,',  NSET=',I5)
      END
EOF
f77 detrend.f -o detrend.x
detrend.x $1 $2 $3
cp fort.2 $2
cp fort.3 $3
rm fort.1 fort.2 fort.3
rm detrend.x detrend.f
exit