#! /bin/sh
#
#    Last change: Mihai Dima 2001
#
if [ $# -lt 4 ]
then
  echo " "
  echo "               *****     CRSRHRT     *****      "
  echo " "
  echo " The program correlates a time series with a series of grids"
  echo " "
  echo "    Syntax: crsrhrt <file1> <file2> <file3> <lag>"
  echo " "
  echo "   <file1> - time series (EXTRA file with one record)"
  echo "   <file2> - time series of grids (EXTRA, transposed file)"
  echo "   <file3> - correlation map"
  echo "   <lag>   - the lag used to perform the correlations"
  echo " "
  exit 0
else
 echo " "
 echo "               *****     CRSRHRT     *****      "
 echo " "
 echo "file1="$1
 echo "file2="$2
 echo "file3="$3
 echo "Lag="$4
 cp $1 fort.1
 cp $2 fort.2
fi
cat > crsrhrt.f << EOF
C
C         The program computes a correlation map 
C         based on a time series 
C         and on a series of maps 
C         (EXTRA transposed files)
C         The correlation is performed for a given lag
C         The time series precede the series of maps.


      PROGRAM CRSRHRT
      PARAMETER (UNDEF=0.9E+10)
      
      REAL X(70000),Y(70000)
      REAL C(70000),XMED

      INTEGER         ID(4)
      INTEGER*4       COUNT
      INTEGER*2       ISTAT,IL      
      INTEGER         ILAG
      INTEGER         NCHAN1,NCHAN2,IDIM1,IDIM2,ILEN1,ILEN2

      
      NCHAN1=1
      CALL PARLES (NCHAN1,IDIM1,ILEN1)

      NCHAN2=2
      CALL PARLES (NCHAN2,IDIM2,ILEN2)

      IF (IDIM1.NE.IDIM2) THEN
        PRINT *,'The time series do not have the same length!'
        STOP
      ENDIF  
      
      LEN=IDIM1

      PRINT *,'The time series precede the series of maps!'
      PRINT *,'Enter the lag for correlation'
      PRINT *, '(INTEGER'
      READ(5,*) ILAG
      PRINT *,"LAG= ",ILAG
     
      REWIND 1
      READ(1) ID
      READ(1)(X(I),I=1,LEN)

C ************************************************************
C  Average is subtracted
C      XMED=0.0
C      DO 21 I=1,LEN
C        XMED=XMED+X(I)
C21    CONTINUE
C      XMED=XMED/REAL(LEN)  
C      DO 23 I=1,LEN
C        X(I)=X(I)-XMED
C23     CONTINUE
C      PRINT *,'OK!'
C ************************************************************
      REWIND 2
      DO 22 I=1,ILEN2

        READ(2) ID
        READ(2)(Y(J),J=1,LEN)
        CALL LAGCORELATION(X,Y,LEN,ILAG,COEFCOREL)
        C(I)=COEFCOREL

22    CONTINUE
                  
      WRITE(3) 1,1,1,ILEN2
      WRITE(3) (C(I),I=1,ILEN2)

      
      PRINT*,' '
      PRINT *,' REGULAR END!'

      END
C*****************************************************************
C    This subroutine verifies the entering file and
C    reads the number of grids points,
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 calculate the coefficient of correlation between the 
C      series X and Y of length IDIM, when ILAG>0 Y-series
C          lags 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)

      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


C*****************************************************************
C    This subroutine verifies the entering file and
C    reads the number of grids points,
C       as well as the number of grids
C
C      SUBROUTINE PARLES (NCHAN,IDIM,ILEN)
C      IMPLICIT REAL (A-H,O-Z)
C      INTEGER I1,I2,I3
C      REAL F
C
C      ILEN=0
C
C      REWIND NCHAN
C      READ(NCHAN,END=10) ITIME,INAME,ILEV,ISIZE
C      GOTO 20
C   10 PRINT *,'INPUT FILE ON TAPE',NCHAN,' IS EMPTY'
C      STOP 'MISTAKE'
C   20 REWIND NCHAN
C
C      IDIM=0  
C    1 READ (NCHAN,END=100) I1,I2,I3,I4
C      IDIM=MAX0(IDIM,I4)
C      READ (NCHAN,END=101) F
C      ILEN=ILEN+1
C      GOTO1
C      
C  101 PRINT *,' MORE HEADER THEN DATA '     
C  100 PRINT 1000,NCHAN,IDIM,ILEN
C      REWIND NCHAN
C      RETURN    
C 1000 FORMAT (1X,' TAPE=',I2,'  IDIM=',I5,'  ILEN=',ILEN)
C      END      
      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 crsrhrt.f -o crsrhrt.x
crsrhrt.x $1 $2 $3 << M
$4
M
cp fort.3 $3
rm fort.1 fort.2 fort.3
rm crsrhrt.x crsrhrt.f
exit