#! /bin/sh
#
#    Last change: Mihai Dima 2001
#
if [ $# -lt 5 ]
then
  echo " "
  echo "        *****     LMEAN     *****      "
  echo " "
  echo "   The program calculates long term means "
  echo "            from monthly data"
  echo " "
  echo " Syntax:"
  echo "  "
  echo "lmean <file1> <file2> <lst_mnth> <frst_yr> <lst_yr>"
  echo " "
  exit 0
else
 echo " "
 echo "        *****     LMEAN     *****      "
 echo " "
 echo "file1="$1
 echo "file2="$2
 echo "Last month="$3
 echo "First year="$4
 echo "Last year="$5
 cp $1 fort.1
fi
cat > lmean.f <<EOF
      PROGRAM PLMEAN
      IMPLICIT REAL (A-H,O-Z)
      
      PRINT *,'0****LMEAN(INPUT,OUTPUT,TAPE1,TAPE2)****'
      PRINT *,'   '
      NCHAN=1
      CALL PARLES (NCHAN,NDIM,NSET)
      CALL LMEAN (NDIM)       
      STOP 
      END
C***************************************************************
      SUBROUTINE PARLES (NCHAN,NDIM,NSET)
      IMPLICIT REAL (A-H,O-Z)
C
      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
      RETURN    
 1000 FORMAT (1X,' TAPE=',I2,',  NDIM=',I5,',  NSET=',I5)
      END
C***************************************************************
      SUBROUTINE LMEAN(IE)
      IMPLICIT REAL (A-H,O-Z)
C
      DIMENSION ID(4)
      PARAMETER (UNDEF=0.9E+10,UNDEFE=0.9*UNDEF)
C IF YOU USE THIS PROGRAM NOT ON A CRAY-2 DELETE THE
C THE FOLLOWING TWO LINES
C      PARAMETER (IE77=65000)
C      DIMENSION X(IE), XM(IE,12), CNM(IE,12)
C      POINTER (NPX,X),(NPXM,XM),(NPCNM,CNM)
C      ALLOCATE (NPX,IE)
C      ALLOCATE (NPXM,IE*12)
C      ALLOCATE (NPCNM,IE*12)
C AND ACTIVATE THE FOLLOWING LINES
      PARAMETER (IE77=20000)
      DIMENSION X(IE77),XM(IE77,12),CNM(IE77,12)
C NOW YOU HAVE DONE ALL CHANGES FOR STANDARD FORTRAN.
      DATA ONE /1./
      IMM = 1
      NFYEAR=-9999
      LYEAR=999999
      WRITE(6,'(A)') 
     +   'ENTER FIRST MONTH ON OUTPUT. FIRST YEAR, LAST YEAR: '
      PRINT *, '   MONTH NFYEAR LYEAR ' 
      READ(5,*,ERR=1,END=1) IMM,NFYEAR,LYEAR
    1 IM=IMM
      REWIND 1
      REWIND 2
      READ(1,END=100) ID
      ITESTY = (ID(1)+1.)/10000.
      ITESTY = ITESTY - IFIX(ITESTY/10000.+0.1)*10000
      IF (NFYEAR.EQ.-9999) NFYEAR=ITESTY
      PRINT *, 'FIRST MONTH ON OUTPUT, FIRST YEAR, LAST YEAR: '
      PRINT *, IMM,NFYEAR,LYEAR  
      IF (ITESTY.GE.1000.AND.NFYEAR.LT.1000) GOTO 103
      IF (ITESTY.LT.1000.AND.NFYEAR.GE.1000) GOTO 103
      IEE = ID(4)
      IF(IEE.GT. IE77) GOTO 101
      REWIND 1
      DO 20 M = 1,12
         DO 20 I = 1,IEE
           CNM(I,M) = 0.
   20      XM(I,M) = 0.
C
C*  READ DATA FROM TAPE1 UNTIL *EOF*.
C
   10 CONTINUE
      READ(1,END=2) ID
      READ(1) (X(I),I=1,ID(4))
C
C*  FIND NUMBER OF MONTH ID(1) = YYMMDD.
C
      MM = (ID(1)+1.)/100.
      MM = MM - IFIX(MM/100.+0.1)*100
      ITESTY = (ID(1)+1.)/10000.
      ITESTY = ITESTY - IFIX(ITESTY/10000.+0.1)*10000
      IF(MM.LT.1 .OR. MM.GT.12) GOTO 102
      IF(ITESTY.GE.NFYEAR.AND.ITESTY.LE.LYEAR)THEN
      DO 30 I = 1,IEE
        IF(X(I) .GT. UNDEFE) GOTO 30
        XM(I,MM) = XM(I,MM) + X(I)
        CNM(I,MM) = CNM(I,MM) + 1.
   30   CONTINUE
      ENDIF
      GOTO 10
C
C*  ALL DATA READ. DIVISION BY NUMBER OF SAMPLES. OUTPUT.
C
    2 CONTINUE
      MCUT=0
      DO 60 M=1,12
       DO 61 I=1,IEE
        IF (CNM(I,M).GT.MCUT)MCUT=AINT(CNM(I,M))
   61  CONTINUE
   60 CONTINUE
      MCUT=(MCUT+3)/4
      DO 40 M = 1,12
        DO 41 I = 1,IEE
          IF(CNM(I,M).LT.MCUT) THEN
            XM(I,M) = UNDEF
          ELSE
            XM(I,M) = XM(I,M)/CNM(I,M)
          ENDIF
   41     CONTINUE
        H1 = 0
        H2 = 0
        IM=0
        DO 42 I = 1,IEE
          IF(CNM(I,M).NE.0)THEN
          IM=IM+1
          H1 = H1 + CNM(I,M)
          ENDIF
   42     IF(CNM(I,M).GT.H2) H2 = CNM(I,M)
        IF (IM.GT.0)THEN
          CNM(1,M) = H1/IM
          CNM(2,M) = H2
        ELSE
          CNM(1,M) = 0
          CNM(2,M) = 0
        ENDIF
   40 CONTINUE
      PRINT '(1X,2A,/,12(1I10,2F10.1,/))', ' MONTHS  MEAN SAMPLES',
     &      ' MAX SAMPLES',(M, (CNM(I,M),I=1,2),M=1,12)
      IM = IMM-1
      DO 50 M = 1,12
        IM = IM + 1
        IF(IM.GT.12) IM = 1
        IF(CNM(2,IM).NE.0.) THEN
          WRITE(2) IM*100, ID(2), ID(3), ID(4)
          WRITE(2) (XM(I,IM),I=1,IEE)
        ENDIF
   50 CONTINUE
      STOP
C
C*  ABNORMAL TERMINATION.
C
  100 WRITE(6,6010) '*LMEAN*: FILE EMPTY'
         GOTO 9910
  101 WRITE(6,6010) '*LMEAN*: DIMENSION TOO LARGE >',IE77
         GOTO 9910
  102 WRITE(6,6010) '*LMEAN*: NUMBER OF MONTH WRONG', ID, MM
         GOTO 9910
  103 WRITE(6,6010) '*LMEAN*: FIRST YEAR ON FILE DONT FIT '
     +              , NFYEAR,ITESTY
         GOTO 9910
 9910 CONTINUE
C
C*  CAUSE AN ABORT.
C
         PRINT *,'*LMEAN*: FATAL END '
         CALL MINONE(ONE,XNEGO)
         Z=SQRT(XNEGO)
         STOP
C
 6010 FORMAT(1X,A/1X,5I10)
      END
      SUBROUTINE MINONE(XIN,XOUT)
        XOUT=-XIN
      RETURN
      END
EOF
f77 lmean.f -o lmean.x
lmean.x $1 $2 << M
$3
$4
$5
M
cp fort.2 $2
rm fort.1 fort.2
rm lmean.x lmean.f
exit