#! /bin/sh
#
#    Last change: Mihai Dima 2001
#
if [ $# -lt 3 ]
then
  echo " "
  echo "           *****     COMPOZ     *****      "
  echo " "
  echo "      The program calculates composite maps"
  echo " "
  echo "      Syntax: compoz <file1> <file2> <file3>"
  echo " "
  echo "   <file1> - time series (EXTRA file with one record)"
  echo "   <file2> - series of grids"
  echo "   <file3> - output composite map-dif sup inf"
  echo " "
  exit 0
else
 echo " "
 echo "           *****     COMPOZ     *****      "
 echo " "
 echo "file1="$1
 echo "file2="$2
 echo "file3="$3
 cp $1 fort.1
 cp $2 fort.2
fi
cat > compoz.f << EOF
C
C    The program calculates composite maps 
C    on the basis of a continuous time series (file1)  
C    and a series of maps (file2)
C    Nr.  of registration in file1 and file 2 
C    remains 
C
      PROGRAM COMPOZ
      PARAMETER(UNDEF=0.9E+10,XPROC=1.5)
      INTEGER H(4),I,J
      REAL F(20000)
      REAL STDV,SN,SMED,SX,XPROC
      REAL X(20000),Y(20000)
      REAL MEDSUP(20000),MEDINF(20000),CPZ(20000)
      REAL CONTORSUP(20000),CONTORINF(20000)

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

C      PRINT *, 'Enter percent of multiplication '
C      PRINT *,'         with standard deviation'
C      PRINT *, '      (p=REAL; 0<=p<1; ex:0.5)'
C      READ(5,*) XPROC
c      XPROC=1.0
      PRINT *,'XPROC=',XPROC     
      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        I have next line uncomment
      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))
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)
C         Here I write how many maps at me on the screen
           write(*,*)'nr maps +'
           write(*,*)contorsup(i)
           write(*,*)'nr maps -'
           write(*,*)contorinf(i)
        ELSE
          CPZ(I)=UNDEF
        ENDIF
      ENDDO

      WRITE(3) H
      WRITE(3) (CPZ(J),J=1,H(4))
      WRITE(3) H
      WRITE(3) (MEDSUP(J)/CONTORSUP(J),J=1,H(4))
      WRITE(3) H
      WRITE(3) (MEDINF(J)/CONTORINF(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
      
C      PRINT *,'NCHAN=',NCHAN
 
      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 compoz.f -o compoz.x
compoz.x $1 $2 $3
cp fort.3 $3
rm fort.1 fort.2 fort.3
rm compoz.x compoz.f
exit