C                                                                       00010001
C***********************************************************************
C                 COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991            *
C***********************************************************************
CTITLESAMEDIA -- SMOOTH THE ARRAY WITH A MOVING MEDIAN FILTER           00020001
CA    AUTHOR       S. BICKEL                                            00030001
CA    DESIGNER     J. J. LEE                                            00040001
CA    LANGUAGE     FORTRAN                                              00050001
CA    SYSTEM       IBM / CRAY                                           00060001
CA    WRITTEN      AUGUST, 1990                                         00070001
C     REVISED      12-21-91 JJC - MODIFIED TO MEET EDP STANDARDS.
C                                                                       00080001
CA                                                                      00090001
CA    CALL SAMEDIA(NY, Y, NF, X, IM)                                    00260000
CA                                                                      00110001
CA         INPUT   NY     = THE NUMBER OF ELEMENTS IN Y               I400120001
CA         IN/OUT  Y      = INPUT ARRAY TO BE SMOOTHED                R400130001
CA         INPUT   NF     = THE LENGTH OF THE RUNNING MEDIAN FILTER   I4
CA         INPUT   X      = WORK REAL ARRAY WITH NF ELEMENTS          R400140001
CA         INPUT   IM     = WORK INTEGER INDEX ARRAY WITH NF ELEMENTS I400180001
C                                                                       00280000
C  PURPOSE: TO SMOOTH THE ARRAY, Y, WITH A MOVING MEDIAN FILTER.        00290000
C**************************************************************
C                                                             *
C SUBROUTINES AND FUNCTIONS CALLED FROM THIS ROUTINE          *
C                                                             *
C SAHPSRT                                                     *
C                                                             *
C**************************************************************
C                                                                       00410000
      SUBROUTINE SAMEDIA(NY,Y,NF,X,IM)                                  00260000
C
      IMPLICIT INTEGER (A-Z)                                            00420000
C
      INTEGER NY,NF,IM(NF),LOCATE                                       00430000
      INTEGER N,N2,N1,N3,J,L1,IP,K,IS,L,JL,JU,JM                        00440000
      REAL Y(NY),X(NF)                                                  00460000
C
      N = MIN ( NF, NY )
      N2 = N / 2
      IF (N2 .EQ. 0) RETURN
      DO 100 J = 1, N
         X(J) = Y(J)
         IM(J) = J
 100  CONTINUE
      CALL SAHPSRT ( N, X, IM )
      DO 120 J = 1, N
         Y(J) = IM(J)
         IM(J) = J
 120  CONTINUE
      CALL SAHPSRT ( N, Y, IM )
      N1 = N2 + 1
      N3 = N2 + N1
      Y(N - N2) = X(N1)
      IF (N3 .GT. N) Y(N - N2) = .5 * (Y(N - N2) + X(N2) )
      DO 140 J = 1, N - N1
         Y(J) = Y(N - N2)
 140  CONTINUE
C                                                                       00640000
C    INITIALIZATION COMPLETE; BEGIN MAIN LOOP. FIND THE INDEX FOY Y(J). 00650000
C                                                                       00660000
      DO 240 J = N + 1, NY
         JL = 0
         JU = N + 1
 160     CONTINUE
         IF (JU - JL .GT. 1) THEN
C          JM=ISHFT(JU+JL,-1)                                           00710000
            JM = (JU + JL )  / 2
            IF (Y(J) .GT. X(JM)) THEN
               JL = JM
            ELSE
               JU = JM
            ENDIF
            GO TO 160
         ENDIF
         L1 = JL + 1
C                                                                       00800000
C    Y(J) HAS BEEN LOCATED. NOW UPDATA THE INDEX ARRAY                  00810000
C                                                                       00820000
         IP = IM(1)
         IS = IP
         IF (IP .GE. L1) IP = IP + 1
         DO 180 K = 1, N - 1
            IF (IM(K + 1) .GE. L1) THEN
               IM(K) = IM(K + 1) + 1
            ELSE
               IM(K) = IM(K + 1)
            ENDIF
            IF (IM(K) .GT. IP) IM(K) = IM(K) - 1
 180     CONTINUE
         IM(N) = L1
         IF (L1 .GT. IP) IM(N) = L1 - 1
C                                                                       00950000
C    IM(J) HAS BEEN UPDATED. NEXT UPDATA THE X ARRAY                    00960008
C                                                                       00970000
         L = L1
         IF (IS .LT. L) L = L - 1
         IF (IS .LT. L) THEN
            DO 200 K = IS, L - 1
               X(K) = X(K + 1)
 200        CONTINUE
         ELSE IF (IS .GT. L) THEN
            DO 220 K = IS, L + 1, -1
               X(K) = X(K - 1)
 220        CONTINUE
         ENDIF
         X(L) = Y(J)
C                                                                       01080000
C    X(J) HAS BEEN UPDATED.                                             01090000
C    NOW SET THE OUTPUT EQUAL TO THE MEDIAN OF THE X ARRAY              01100000
C                                                                       01110000
         Y(J - N2) = X(N1)
         IF (N3 .GT. N) Y(J - N2) = .5 * (Y(J - N2) + X(N2) )
 240  CONTINUE
C                                                                       01140000
C    MAIN LOOP COMPLETE; COMPLETE THE END OF THE ARRAY                  01150000
C                                                                       01160000
      DO 260 J = NY - N2 + 1, NY
         Y(J) = Y(NY - N2)
 260  CONTINUE
      IF (N3 .GT. N) THEN
         DO 280 J = NY, 2, -1
            Y(J) = (Y(J) + Y(J - 1) )  / 2.
 280     CONTINUE
      ENDIF
C
      RETURN
      END