#! /bin/sh
#
#    Last change: Mihai Dima 2001
#
if [ $# -lt 5 ]
then
  echo " "
  echo "         *****     CROSSCOR     *****      "
  echo " "
  echo "    The program cross correlates two time series "
  echo " "
  echo " Syntax: crosscor <file1> <file2> <file3> <file4> <lag>"
  echo " "
  echo " <file1> - first time series (EXTRA file with one record)"
  echo " <file2> - second time series (EXTRA file with one record)"
  echo " <file3> - output cross correlation in EXTRA format"
  echo " <file4> - output cross correlation in TEXT format"
  echo " "
  exit 0
else
 echo " "
 echo "         *****     CROSSCOR     *****      "
 echo " "
 echo "file1="$1
 echo "file2="$2
 echo "file3="$3
 echo "file4="$4
 echo "Maximum lag="$5
 cp $1 fort.1
 cp $2 fort.2
fi
cat > crosscor.f << EOF
      PROGRAM CROSSCOR

      PARAMETER(NCHAN1=1,NCHAN2=1,UNDEF=0.9E+10)

c     The program calculates the correlation coefficients 
C     between two time series (file EXTRA),
c     for different lags,
c     from  -LAGMAX  to LAGMAX
c     The obtained series of correlation coefficients is
c     stored an EXTRA and a TEXT files.
c     NT - length of time series
      
      REAL X(10000),Y(10000)
      REAL C(10000)

      INTEGER IDAT,IVAR,ILEV,ILEN

      INTEGER*4       COUNT
      INTEGER*2       ISTAT,IL      
      INTEGER         ERROR,IV2
      INTEGER         IDIM1,IDIM2,ILEN1,ILEN2

      PRINT *,'GET THE MAXIMUM LAG'
C      print*,'lagmax'
      READ(5, *) lagmax
      
      CALL PARLES (NCHAN1,IDIM1,ILEN1)
      CALL PARLES (NCHAN2,IDIM2,ILEN2)

      IF (IDIM1.NE.IDIM2) THEN
        PRINT *,'The time series have different lengths!'
        STOP
      ENDIF  

      PRINT *,' LAGMAX=',LAGMAX

      LEN=IDIM1
      LENCOR=2*LAGMAX+1

      OPEN(4,FORM='FORMATTED',STATUS='UNKNOWN')

      READ(1) IDAT,IVAR,ILEV,ILEN
      READ(1)(X(I),I=1,LEN)
      READ(2) IDAT,IVAR,ILEV,ILEN
      READ(2)(Y(J),J=1,LEN)

      DO 22 LAG=-LAGMAX,LAGMAX

C  Calculates the coefficient of correlation for each lag
C**********************************************************************

        CALL LAGCORELATION(X,Y,LEN,LAG,COEFCOREL)

        C(LAGMAX+LAG+1)=COEFCOREL

C**********************************************************************

22    CONTINUE

      WRITE(3) 1,1,1,LENCOR
      WRITE(3) (C(I),I=1,LENCOR)

      DO I=1,LENCOR
        WRITE(4,'(I4,1X,E12.5)') I-LAGMAX-1,C(I)
      END DO

      CLOSE(3)
      CLOSE(4)

      PRINT *,'REGULAR END!'

      END
C**********************************************************************
C   This subroutine verified the entering file and
C    read the number of points of grids,
C       as well as the number of grids
C
      SUBROUTINE PARLES (NCHAN,IDIM,ILEN)
      IMPLICIT REAL (A-H,O-Z)
      INTEGER I1,I2,I3
      REAL F

      ILEN=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

      IDIM=0  
    1 READ (NCHAN,END=100) I1,I2,I3,I4
      IDIM=MAX0(IDIM,I4)
      READ (NCHAN,END=101) F
      ILEN=ILEN+1
      GOTO1
      
  101 PRINT *,' MORE HEADER THEN DATA '     
  100 PRINT 1000,NCHAN,IDIM,ILEN
      REWIND NCHAN
      RETURN    
 1000 FORMAT (1X,' TAPE=',I2,',  IDIM=',I5,',  ILEN=',I5)
      END      
C**********************************************************************
C   The procedure calculates the coefficient of correlation between
C      series X and Y of length IDIM, when ILAG>0 Y- series lags
C          X- series with LAG-values, and vice versa
C
      SUBROUTINE LAGCORELATION(XP,YP,IDIM,ILAG,COEFCOREL)
      REAL XP(10000),YP(10000),XW(10000),YW(10000)      
      REAL MEDX,MEDY,SR,SR1,SR2,R
      INTEGER LENREAL,IL,IDEF,IUNDEF

      UNDEF=0.9E+10
      LENREAL=IDIM-ABS(ILAG)
      
C      PRINT *,'  IDIM=',IDIM,'  ILAG=',ILAG
C      PRINT *,'  LENREAL=',LENREAL
      
      IF (ILAG.EQ.0) THEN
        DO IL=1,LENREAL,1
          XW(IL)=XP(IL)
          YW(IL)=YP(IL)
        END DO
      ENDIF

      IF (ILAG.GT.0) THEN
        DO IL=1,LENREAL,1
          XW(IL)=XP(IL)
          YW(IL)=YP(IL+ABS(ILAG))
        END DO
      ENDIF

      IF (ILAG.LT.0) THEN
        DO IL=1,LENREAL,1
          XW(IL)=XP(IL+ABS(ILAG))
          YW(IL)=YP(IL)
        END DO
      ENDIF

      IDEF=0
      IUNDEF=0
      MEDX=0.0
      MEDY=0.0

      DO IL=1,LENREAL,1
        IF ((XW(IL).NE.UNDEF).AND.(YW(IL).NE.UNDEF)) THEN
          MEDX=MEDX+XW(IL)
          MEDY=MEDY+YW(IL)
          IDEF=IDEF+1
        ELSE
          IUNDEF=IUNDEF+1
        ENDIF
      END DO

      IF (IUNDEF.LT.INT(LENREAL/2).AND.(IDEF.NE.0)) THEN
        MEDX=MEDX/REAL(IDEF)
        MEDY=MEDY/REAL(IDEF)
      ELSE
        MEDX=0.0
        MEDY=0.0
      ENDIF

      IDEF=0
      IUNDEF=0
      SR=0.0
      SR1=0.0
      SR2=0.0

      IF ((MEDX.NE.0.0).AND.(MEDY.NE.0.0)) THEN
        DO IL=1,LENREAL,1
          IF ((XW(IL).NE.UNDEF).AND.(YW(IL).NE.UNDEF)) THEN
            SR=SR+(XW(IL)-MEDX)*(YW(IL)-MEDY)
            SR1=SR1+(XW(IL)-MEDX)*(XW(IL)-MEDX)
            SR2=SR2+(YW(IL)-MEDY)*(YW(IL)-MEDY)
            IDEF=IDEF+1
          ELSE
            IUNDEF=IUNDEF+1
          ENDIF
        END DO

        IF ((SR1.EQ.0).OR.(SR2.EQ.0).OR.
     *    (IUNDEF.GT.INT(LENREAL/2))) THEN
          R=0.0
        ELSE
          R=SR/SQRT(SR1*SR2)
        ENDIF
      ELSE
         R=UNDEF
      ENDIF

      COEFCOREL=R
      
      RETURN

      END
EOF
f77 crosscor.f -o crosscor.x
crosscor.x $1 $2 $3 $4<< M
$5
M
cp fort.3 $3
cp fort.4 $4
rm fort.1 fort.2 fort.3 fort.4
rm crosscor.x crosscor.f
exit