#! /bin/sh
#
#    Last change: Mihai Dima 2001
#
if [ $# -lt 4 ]
then
  echo " "
  echo "           *****     COMPZLAG     *****      "
  echo " "
  echo "      The program calculates composite maps"
  echo " "
  echo "    Syntax: compoz <file1> <file2> <file3> <lag>"
  echo " "
  echo "   <file1> - time series (EXTRA file with one record)"
  echo "   <file2> - series of grids"
  echo "   <file3> - output composite map"
  echo "   <lag>   - lag for which the composites are calculated"
  echo " "
  exit 0
else
 echo " "
 echo "           *****     COMPZLAG     *****      "
 echo " "
 echo "file1="$1
 echo "file2="$2
 echo "file3="$3
 echo "Lag="$4
 cp $1 fort.1
 cp $2 fort.2
fi
cat > compzlag.f << EOF
C
C    The program calculates composite maps 
C    based on a non transposed time series (file1)  
C    and on a series of maps (file2)
C    The Nr. of records in file1 and file 2 
C    must be same
C    The composition is performed with a delay ILAG 
C    The maps lag the values in the time series with ILAG
C
      PROGRAM COMPZLAG
      PARAMETER(UNDEF=0.9E+10)
      INTEGER H(4),I,J
      REAL F(70000)
      REAL STDV,SN,SMED,SX,XPROC
      REAL X(70000),Y(70000)
      REAL MEDSUP(70000),MEDINF(70000),CPZ(70000)
      REAL CONTORSUP(70000),CONTORINF(70000)

      NCHAN=1
      CALL PARLES (NCHAN,IDIM,NSET)

      XPROC=0.5
C      PRINT *, 'Enter the percent to be multiplied’
C      PRINT *, '         with the standard deviation'
C      PRINT *, '      (p=REAL; 0<=p<1; ex:0.5)'
C      READ(5,*) XPROC
C      PRINT *, 'Enter the delay between time series and maps'
C      PRINT *, '          (COMPLETE,POZITIV)'
      READ(5,*) ILAG

      READ (1) H
      IDIM=H(4)            
      SMED=0.0
      SX=0.0
      SN=0.0 

      REWIND 1
60    READ (1,END=61) H
      READ (1) (F(J),J=1,IDIM)
      DO I=1,H(4)
        IF (F(I).NE.UNDEF) THEN
          SMED=SMED+F(I)
          SN=SN+1.0
        ENDIF
      ENDDO
      GOTO 60

61    CONTINUE     
      IF (SN.NE.0.0) THEN
         SMED=SMED/SN   
      ENDIF
      
      REWIND 1      
50    READ (1,END=51) H
      READ (1) (F(J),J=1,IDIM)
      DO 43 I=1,H(4)
        IF (SMED.NE.0.0) THEN
          SX=SX+(F(I)-SMED)*(F(I)-SMED)
        ELSE
          SX=UNDEF  
        ENDIF
43    CONTINUE
      GOTO 50

51    CONTINUE     
      DO I=1,H(4)
         IF (SX.NE.UNDEF) THEN
            STDV=SQRT(SX/(SN-1.0))
         ELSE
            STDV=UNDEF      
         ENDIF
      ENDDO
C      PRINT *,' STDV=',STDV

      READ (2,END=71) H
      DO J=1,H(4)
        MEDSUP(J)=0.0
        MEDINF(J)=0.0
        CONTORSUP(J)=0.0
        CONTORINF(J)=0.0
      ENDDO

      L=1
      REWIND 1
      REWIND 2
      READ (1,END=71) H
      READ (1,END=71) (X(I),I=1,H(4))

      DO I=1,ILAG
        READ (2) H
        READ (2) (Y(J),J=1,H(4))
      ENDDO

70    READ (2,END=71) H
      READ (2) (Y(I),I=1,H(4))
      IF (X(L).GE.(XPROC*STDV)) THEN
        DO J=1,H(4)
          IF (Y(J).NE.UNDEF) THEN
            MEDSUP(J)=MEDSUP(J)+Y(J)
            CONTORSUP(J)=CONTORSUP(J)+1.0
          ENDIF
        ENDDO
      ENDIF
      IF (X(L).LE.(-XPROC*STDV)) THEN
        DO J=1,H(4)
          IF (Y(J).NE.UNDEF) THEN
            MEDINF(J)=MEDINF(J)+Y(J)
            CONTORINF(J)=CONTORINF(J)+1.0
          ENDIF
        ENDDO
      ENDIF
      L=L+1
      GOTO 70

71    CONTINUE

      DO I=1,H(4)
        IF ((CONTORSUP(I).NE.0.0).AND.(CONTORINF(I).NE.0.0)) THEN
          CPZ(I)=MEDSUP(I)/CONTORSUP(I)-MEDINF(I)/CONTORINF(I)
        ELSE
          CPZ(I)=UNDEF
        ENDIF
      ENDDO

      WRITE(3) H
      WRITE(3) (CPZ(J),J=1,H(4))

      CLOSE(3)
      CLOSE(2)
      CLOSE(1)

      STOP 'REGULAR END!'

      END
C******************************************************************
      SUBROUTINE PARLES (NCHAN,NDIM,NSET)
      IMPLICIT REAL (A-H,O-Z)

      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
      REWIND NCHAN
      RETURN    
 1000 FORMAT (1X,' TAPE=',I2,',  NDIM=',I5,',  NSET=',I5)
      END
C******************************************************************
EOF
f77 compzlag.f -o compzlag.x
compzlag.x $1 $2 $3 << EOF
$4
EOF
cp fort.3 $3
rm fort.1 fort.2 fort.3
rm compzlag.x compzlag.f
exit