C @(#)smosub3.for	13.1.1.1 (ESO-DMD) 06/02/98 18:11:31
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
      SUBROUTINE NORMIT(A,NDIM)
C
      IMPLICIT NONE
C
      INTEGER  NDIM,N
      REAL     A(NDIM),R
      DOUBLE PRECISION   SUM,DR
C
      SUM = 0.D0
      DO 1000, N=1,NDIM
         SUM = SUM + A(N)
1000  CONTINUE
C 
      IF (SUM.GT.10.D-20) THEN
         DR = 1.D0 / SUM
         R = SNGL(DR)
         DO 2000, N=1,NDIM
            A(N) = A(N) * R
2000     CONTINUE
      ENDIF
C 
      RETURN
      END

      SUBROUTINE TEMPLB(FLAG,A,B,NPIX,WEITS,RADIUS,CUTS,STAT)
C
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C	
C.IDENTIFICATION
C  subroutine TEMPLB			version 1.50	860610
C  K. Banse				ESO - Garching
C
C.KEYWORDS
C  template matching, neighbourhood
C
C.PURPOSE
C replace pixels by weighted sum of nw*nw matrix (nw = 2*radius + 1)
C
C.ALGORITHM
C  filter image by taking the weighted sum W(x,y) of all points in neighbourhood
C  of x,y with given radius, and replace original value I(x,y) by W(x,y)
C	
C.INPUT/OUTPUT
C  call as TEMPLATE(A,B,NPIX,WEITS,RADIUS,STAT)
C
C input par:
C  FLAG:	I*4 		function flag:
C				= 1 for Gauss filter
C				= 2 for digital filter
C				= 3 for max filter
C				= 4 for min filter
C  A:		R*4 array	input array
C  B:		R*4 array	output array
C  NPIX:	I*4 array	no. of pixels per line, no. of lines
C  WEITS:	R*4 array	weights for NWX*NWY matrix, ordered as
C				1 2 3 ... NWX
C				NWX+1 ... 2NWX
C				............
C				(NWY-1)*NWX+1 ... NWX*NWY
C  RADIUS:	I*4 array	radius of template matrix in x,y
C	
C  output par:
C  CUTS:	R*4 array	cut values of result frame
C  STAT:	I*4		return status, should be = 0, else problems...
C	
C --------------------------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER     FLAG,NPIX(*),RADIUS(*)
      INTEGER     NX,NY,DIAMX,DIAMY,RADYNX
      INTEGER     N,NP,NL,BINDX,STAT
      INTEGER     WINDX,OFFSET,OFF
      INTEGER     NRDXP1,RADXP1,NXRADX,LEFT
C
      REAL        A(*),B(*),WEITS(*),CUTS(*)
      REAL        TEST
C 
      DOUBLE PRECISION  DTEST
C
C  init variables + test size of template
      BINDX = 1                                !index for frame B
      NX = NPIX(1)
      NY = NPIX(2)
      DIAMX = RADIUS(1) * 2
      DIAMY = RADIUS(2) * 2
      IF ((DIAMX.GE.NX).OR.(DIAMY.GE.NY)) THEN        !data matrix at least NW
         STAT = 1
         RETURN
      ELSE
         STAT = 0
      ENDIF
C
      RADYNX = RADIUS(2) * NX                  !all following constants for
      RADXP1 = RADIUS(1) + 1                   !speed reasons...
      NXRADX = NX - RADIUS(1)
      NRDXP1 = NXRADX + 1
      LEFT = RADYNX + RADIUS(1)            !lower left corner of template matrix
      OFFSET = RADYNX
C
C  branch according to filter type
      IF (FLAG.GE.3) GOTO 5000
C
C*******
C
C  main loop over inner lines
C
C*******
C
      DO 2000, NL=RADIUS(2)+1,NY-RADIUS(2)            !step through each line
C
C  step through each pixel of data frame
C  i.e., OFFSET+RADIUS(1)+1  >>>  OFF SET+NX-RADIUS(1)
C 
         DO 1500, NP=OFFSET+RADXP1,OFFSET+NXRADX 
            DTEST = 0.D0                              !clear sum
            WINDX = 1                                !index for weight matrix
C
C  sum + weigh inner part
            DO 1000, OFF=NP-LEFT,NP+LEFT,NX        !step through lines of window
               DO 800, N=OFF,OFF+DIAMX
                  DTEST = DTEST + ( A(N) * WEITS(WINDX) )
                  WINDX = WINDX + 1
800            CONTINUE
1000        CONTINUE
C
            B(BINDX) = DTEST                     !fill output frame
            BINDX = BINDX + 1
1500     CONTINUE
C
C  update main offset
         OFFSET = OFFSET + NX
C
2000  CONTINUE
C
C finally update the cuts
      CALL CHCUTS(B,BINDX,CUTS)
      RETURN
C
C  here for min/max filter
C 
5000  IF (FLAG.EQ.3) THEN
C 
C  max filter
         DO 7000, NL=RADIUS(2)+1,NY-RADIUS(2)            !step through each line
C
C  step through each pixel of data frame
C  i.e., OFFSET+RADIUS(1)+1  >>>  OFF SET+NX-RADIUS(1)
C
            DO 6000, NP=OFFSET+RADXP1,OFFSET+NXRADX
               TEST = A(NP-LEFT)                        !take lower left pixel
C
C  find max in inner part
               DO 5500, OFF=NP-LEFT,NP+LEFT,NX    
                  DO 5400, N=OFF,OFF+DIAMX
                     IF (A(N).GT.TEST) TEST = A(N)
5400              CONTINUE
5500           CONTINUE
C
               B(BINDX) = TEST                     !fill output frame
               BINDX = BINDX + 1
6000        CONTINUE
C
            OFFSET = OFFSET + NX                   !update main offset
7000     CONTINUE
C 
      ELSE
C 
C  min filter
         DO 9000, NL=RADIUS(2)+1,NY-RADIUS(2)            !step through each line
C
C  step through each pixel of data frame
C  i.e., OFFSET+RADIUS(1)+1  >>>  OFF SET+NX-RADIUS(1)
C
            DO 8800, NP=OFFSET+RADXP1,OFFSET+NXRADX
               TEST = A(NP-LEFT)                        !take lower left pixel
C
C  find min in inner part
               DO 8600, OFF=NP-LEFT,NP+LEFT,NX   
                  DO 8400, N=OFF,OFF+DIAMX
                     IF (A(N).LT.TEST) TEST = A(N)
8400              CONTINUE
8600           CONTINUE
C
               B(BINDX) = TEST                     !fill output frame
               BINDX = BINDX + 1
8800        CONTINUE
C
            OFFSET = OFFSET + NX                   !update main offset
9000     CONTINUE
      ENDIF
C
C finally update the cuts
      CALL CHCUTS(B,BINDX,CUTS)
      RETURN
C
      END

      SUBROUTINE CHCUTS(A,ASZ1,CUTS)
C 
      INTEGER     ASZ1
      INTEGER     N,NDIM
C 
      REAL        A(*),CUTS(*)
      REAL        TEST
C 
      NDIM = ASZ1 - 1
      IF (NDIM.LT.3) RETURN
C 
      IF (A(1).GT.A(2)) THEN
         CUTS(1) = A(2)
         CUTS(2) = A(1)
      ELSE
         CUTS(1) = A(1)
         CUTS(2) = A(2)
      ENDIF
C 
      DO 1000, N=3,NDIM
         TEST = A(N)
         IF (TEST.LT.CUTS(1)) THEN
            CUTS(1) = TEST
         ELSE IF (TEST.GT.CUTS(2)) THEN
            CUTS(2) = TEST
         ENDIF
1000  CONTINUE
C 
      RETURN
      END

      SUBROUTINE EXPAX(IMNOA,OFFPIX,OUTFR,INBUF,RSTAT)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine EXPX
C  K. Banse				ESO - Garching
C  1.00    910517
C
C.KEYWORDS
C  border regions
C
C.PURPOSE
C  expand image by duplicating border lines + columns
C
C.ALGORITHM
C  get the no. of columns to expand at left + at right
C  get the no. of lines to expand in the beginning + end
C  flip columns in x, lines in y and add to image
C
C.INPUT/OUTPUT
C  call via
C   EXPAX(IMNOA,OFFPIX,OUTFR,INBUF,RSTAT)
C
C  input par:
C  IMNOA	integer         file no. of input frame
C  OFFPIX	integer		offset for 1. pixel to map
C  OUTFR	char. array	output frame
C  INBUF	integer array	xl,xr,yb,ye: no. of cols (left/right)
C					     no. of lines (beginning/end)
C  output par:
C  RSTAT	integer         return status, 0 = o.k.
C
C.VERSIONS
C  1.00   creation
C	
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      INBUF(4),RSTAT,OFFPIX
      INTEGER      IAV,N,NN,LOOPLM,RMAIND,CHUNKI,CHUNKO
      INTEGER      NAXISA,MOFF,LOFF,NOLIN
      INTEGER*8    PNTRX,PNTRY
      INTEGER      IMNOA,IMNOC,IMNOX,IMNOY
      INTEGER      TRUSIZ,SIZE,STAT
      INTEGER      NPIXA(2),NPIXC(2)
      INTEGER      UNI(1),NULO,MADRID(1)
C
      CHARACTER*(*) OUTFR
      CHARACTER    CUNITA*64,IDENTA*72
C 
      REAL         CUTVAL(4)
C
      DOUBLE PRECISION    STEPA(2),STARTA(2)
      DOUBLE PRECISION    STARTC(2)
C
      COMMON  /VMR/  MADRID
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      DATA    NPIXA  /1,1/,   NPIXC /1,1/
      DATA    STARTA /0.0,0.0/,  STEPA  /1.0,1.0/
      DATA    CUNITA /' '/,   IDENTA /' '/
C 
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C 
      RSTAT = 0
      CALL STDRDI(IMNOA,'NAXIS',1,1,IAV,NAXISA,UNI,NULO,STAT)
      IF (NAXISA.GT.2) NAXISA = 2                                !max 2-dim 
      CALL STDRDI(IMNOA,'NPIX',1,NAXISA,IAV,NPIXA,UNI,NULO,STAT)
      CALL STDRDD(IMNOA,'START',1,NAXISA,IAV,STARTA,UNI,NULO,STAT)
      CALL STDRDD(IMNOA,'STEP',1,NAXISA,IAV,STEPA,UNI,NULO,STAT)
      CALL STDRDC(IMNOA,'IDENT',1,1,72,IAV,IDENTA,UNI,NULO,STAT)
      N = (NAXISA+1) * 16
      CALL STDRDC(IMNOA,'CUNIT',1,1,N,IAV,CUNITA,UNI,NULO,STAT)
      CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTVAL,UNI,NULO,STAT)
C
C  get size of input + output frame
      SIZE = NPIXA(1) * NPIXA(2)
      IF ( (INBUF(1).GT.NPIXA(1)) .OR.
     +     (INBUF(2).GT.NPIXA(1)) .OR.
     +     (INBUF(1).LT.0) .OR.
     +     (INBUF(2).LT.0) ) THEN
         RSTAT = 1
         RETURN
      ENDIF
C 
      NPIXC(1) = NPIXA(1) + INBUF(1) + INBUF(2)
      IF (NAXISA.GT.1) THEN
         IF ( (INBUF(3).GT.NPIXA(2)) .OR.
     +        (INBUF(4).GT.NPIXA(2)) .OR.
     +        (INBUF(3).LT.0) .OR.
     +        (INBUF(4).LT.0) ) THEN
            RSTAT = 2
            RETURN
         ENDIF
         NPIXC(2) = NPIXA(2) + INBUF(3) + INBUF(4)
      ENDIF
C
C create result frame
      IF (OFFPIX.LE.0) THEN
         TRUSIZ = NPIXC(1) * NPIXC(2)
         CALL STFCRE(OUTFR,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               TRUSIZ,IMNOC,STAT)
C 
         CALL STDWRI(IMNOC,'NAXIS',NAXISA,1,1,UNI,STAT)
         CALL STDWRI(IMNOC,'NPIX',NPIXC,1,NAXISA,UNI,STAT)
         STARTC(1) = STARTA(1) - (INBUF(1)*STEPA(1))
         IF (NAXISA.GT.1)
     +      STARTC(2) = STARTA(2) - (INBUF(3)*STEPA(2))
         CALL STDWRD(IMNOC,'START',STARTC,1,NAXISA,UNI,STAT)
         CALL STDWRD(IMNOC,'STEP',STEPA,1,NAXISA,UNI,STAT)
         CALL STDWRC(IMNOC,'IDENT',1,IDENTA,1,72,UNI,STAT)
         N = (NAXISA+1) * 16
         CALL STDWRC(IMNOC,'CUNIT',1,CUNITA,1,N,UNI,STAT)
         CALL STDWRR(IMNOC,'LHCUTS',CUTVAL,1,4,UNI,STAT)
      ELSE
         CALL STFOPN(OUTFR,D_R4_FORMAT,0,F_IMA_TYPE,IMNOC,STAT)
      ENDIF
C 
C  get buffersize of at least one line but no more than total size
      NN = 40000                                !first guess ...
      IF (NN.LT.NPIXA(1)) NN = NPIXA(1)
      IF (NN.GT.SIZE) NN = SIZE
      NOLIN = NN / NPIXA(1)                     !buffer = NOLIN * NPIXA(1)
      CHUNKI = NOLIN * NPIXA(1)
      CHUNKO = NOLIN * NPIXC(1)
      IF (NAXISA.GT.1) THEN
         LOOPLM = NPIXA(2) / NOLIN              !get no. of passes through loop
         RMAIND = NPIXA(2) - (LOOPLM*NOLIN)     !and remainder
      ELSE
         LOOPLM = 1
         RMAIND = 0
      ENDIF
C 
      CALL STFCRE('expaxi',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            CHUNKI,IMNOX,STAT)
      CALL STFMAP(IMNOX,F_X_MODE,1,CHUNKI,IAV,PNTRX,STAT)
      CALL STFCRE('expaxo',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            CHUNKO,IMNOY,STAT)
      CALL STFMAP(IMNOY,F_X_MODE,1,CHUNKO,IAV,PNTRY,STAT)
C 
C  now we loop through input frame
      MOFF = 1 + OFFPIX                        !point to 1. pixel in input fr.
      IF (NAXISA.GT.1) THEN                    !point to 1. pixel in output fr.
         LOFF = (INBUF(3)*NPIXC(1)) + 1       
      ELSE
         LOFF = 1
      ENDIF
C 
      DO 2000, N=1,LOOPLM
         CALL STFGET(IMNOA,MOFF,CHUNKI,IAV,MADRID(PNTRX),STAT)    !data in
         CALL EXPA1(MADRID(PNTRX),MADRID(PNTRY),NOLIN,NPIXA(1),
     +              NPIXC(1),INBUF)
         CALL STFPUT(IMNOC,LOFF,CHUNKO,MADRID(PNTRY),STAT)          !data out
         MOFF = MOFF + CHUNKI
         LOFF = LOFF + CHUNKO
2000  CONTINUE
C 
      IF (RMAIND.GT.0) THEN
         CHUNKI = RMAIND * NPIXA(1)
         CHUNKO = RMAIND * NPIXC(1)
         CALL STFGET(IMNOA,MOFF,CHUNKI,IAV,MADRID(PNTRX),STAT) 
         CALL EXPA1(MADRID(PNTRX),MADRID(PNTRY),RMAIND,NPIXA(1),
     +              NPIXC(1),INBUF)
         CALL STFPUT(IMNOC,LOFF,CHUNKO,MADRID(PNTRY),STAT)  
      ENDIF
C 
C  for 1-dim frames we are finished
      CALL STFCLO(IMNOX,STAT)
      CALL STFCLO(IMNOY,STAT)
      IF (NAXISA.EQ.1) GOTO 8000
C 
C  get buffersize of INBUF(3/4)*2 lines
      IF (INBUF(3).GE.INBUF(4)) THEN
         NN = INBUF(3)
      ELSE
         NN = INBUF(4)
      ENDIF
      IF (NN.EQ.0) GOTO 8000                !nothing to do in Y
C 
      NN = (2*NN + 1) * NPIXC(1)
      CALL STFCRE('expaxz',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            NN,IMNOX,STAT)
      CALL STFMAP(IMNOX,F_X_MODE,1,NN,IAV,PNTRX,STAT)
C  
      NN = 2*INBUF(3) + 1
      CHUNKI = NN * NPIXC(1)
      IF (CHUNKI.GT.0) THEN
         CALL STFGET(IMNOC,1,CHUNKI,IAV,MADRID(PNTRX),STAT) 
         CALL EXPA2(1,MADRID(PNTRX),NPIXC,INBUF(3))
         CALL STFPUT(IMNOC,1,CHUNKI,MADRID(PNTRX),STAT)  
      ENDIF
C  
      NN = 2*INBUF(4) + 1
      CHUNKI = NN * NPIXC(1)
      IF (CHUNKI.GT.0) THEN
         MOFF = ((NPIXC(2)-NN)*NPIXC(1)) + 1
         CALL STFGET(IMNOC,MOFF,CHUNKI,IAV,MADRID(PNTRX),STAT) 
         CALL EXPA2(2,MADRID(PNTRX),NPIXC,INBUF(4))
         CALL STFPUT(IMNOC,MOFF,CHUNKI,MADRID(PNTRX),STAT)  
      ENDIF
      CALL STFCLO(IMNOX,STAT)
C 
8000  CALL STFCLO(IMNOC,STAT)
      RETURN
C 
      END

      SUBROUTINE EXPA1(X,Y,NOLIN,NPIXA,NPIXC,INBUF)
C 
      IMPLICIT NONE
C 
      REAL     X(*),Y(*)
C 
      INTEGER  NOLIN,NPIXA,NPIXC,INBUF(*)
      INTEGER  N,NN,IDX,KDX,IOFF,KOFF
C 
C  loop over all lines
      IOFF = 0
      KOFF = 0
      DO 5000, NN=1,NOLIN
C 
C fill first Y pixels with first flipped X pixels
         IF (INBUF(1).GT.0) THEN
            IDX = 1 + IOFF
            KDX = INBUF(1) + KOFF + 1
            DO 1000, N=1,INBUF(1)
               Y(IDX) = X(KDX)
               IDX = IDX + 1
               KDX = KDX - 1
1000        CONTINUE
         ELSE
            IDX = 1 + IOFF
            KDX = 1 + KOFF
         ENDIF
C 
C copy center Y pixels from all X pixels
         DO 2000, N=1,NPIXA
            Y(IDX) = X(KDX)
            IDX = IDX + 1
            KDX = KDX + 1
2000     CONTINUE
C 
C fill last Y pixels with last flipped X pixels
         KDX = KDX - 1
         IF (INBUF(2).GT.0) THEN
            DO 3000, N=1,INBUF(2)
               KDX = KDX - 1
               Y(IDX) = X(KDX)
               IDX = IDX + 1
3000        CONTINUE
         ENDIF
C 
         IOFF = IOFF + NPIXC
         KOFF = KOFF + NPIXA
5000  CONTINUE
C 
      RETURN
      END

      SUBROUTINE EXPA2(FLAG,X,NPIXC,INB)
C 
      IMPLICIT NONE
C 
      REAL     X(*)
C 
      INTEGER  FLAG,NPIXC(*),INB
      INTEGER  N,NN,IDX,KDX,IOFF,KOFF,MX,MMX
C 
      MX = NPIXC(1)
C 
      IF (FLAG.EQ.1) THEN
         IOFF = 0
         KOFF = 2*INB * MX
         MMX = MX
      ELSE
         IOFF = 2*INB * MX
         KOFF = 0
         MMX = -MX
      ENDIF
C 
      DO 2000, NN=1,INB
         IDX = IOFF + 1
         KDX = KOFF + 1
         DO 1000, N=1,MX
            X(IDX) = X(KDX)
            IDX = IDX + 1
            KDX = KDX + 1
1000     CONTINUE
C 
         IOFF = IOFF + MMX
         KOFF = KOFF - MMX
2000  CONTINUE
C 
      RETURN
      END

      SUBROUTINE SORTIT(RA,N,NH,LL)
C 
C  N  = full dimension of RA
C  NH = index of median - we only sort fully until NH
C  then we insert
C 
C  this is the Heapsort algorithm from "Numerical Recipes", page 231
C 
      IMPLICIT NONE
C 
      INTEGER   N,NH,LL
      INTEGER   L,IR,J,I
      INTEGER   JL,JU,JM,M
C 
      REAL      RA(N),RRA
C 
      L = LL                      !that is L = NH/2 + 1  
      IR = NH
C 
10    IF (L.GT.1) THEN            !still hiring
         L = L - 1
         RRA = RA(L)
      ELSE                        !in retirement + promotion phase
         RRA = RA(IR)             !clear a space at end of array
         RA(IR) = RA(1)           !retire the top of the heap into it
         IR = IR - 1           !decrease the size of the corporation
         IF (IR.EQ.1) THEN        !done with last promotion
            RA(1) = RRA           !the least competent guy of all...
            GOTO 1000             !now move to insertion part
         ELSE
         ENDIF
      ENDIF
C 
      I = L                       !in hiring as well as in promotion phase
      J = L + L              !we here set up to sift down RRA to its right level
C 
20    IF (J.LE.IR) THEN           !"Do while J<=IR"
         IF (J.LT.IR) THEN
            IF (RA(J).LT.RA(J+1)) J = J + 1   !compare to the better underling
         ENDIF
         IF (RRA.LT.RA(J)) THEN   !demote RRA
            RA(I) = RA(J)
            I = J
            J = J + J
         ELSE
            J = IR + 1
         ENDIF
         GOTO 20
      ENDIF
C 
      RA(I) = RRA                 !put RRA into its slot
      GOTO 10
C 
C 
1000  I = NH + 1
      DO 2000 L=I,N
         IF (RA(L).LT.RA(NH)) THEN
            RRA = RA(L)
            JL = 0
            JU = I
C 
1010        JM = (JU+JL) / 2
            IF (RRA.LT.RA(JM)) THEN
               JU = JM
            ELSE
               JL = JM
            ENDIF
            IF ((JU-JL).GT.1) GOTO 1010               !loop
C 
CC            DO 1100 M=L,JL+2,-1          !RA kept intact
            DO 1100 M=NH,JL+2,-1           !ojo - RA not kept intact that way!
               RA(M) = RA(M-1)             !but faster
1100        CONTINUE
            RA(JL+1) = RRA
         ENDIF
2000  CONTINUE
C 
      RETURN
      END

      SUBROUTINE XSORT1(RA,NDIM,NH,RINDX)
C 
C  RA  = real array of which the elements were used for sorting INDARR
C  NDIM  = full dimension of RA
C  NH = index of median - we only sort fully until NH
C  then we insert
C 
C  this is the Heapsort algorithm from "Numerical Recipes", page 233
C 
      IMPLICIT NONE
C 
      INTEGER   NDIM,NH,RINDX
      INTEGER   L,IR,J,I
      INTEGER   INDXT
C 
      INTEGER   INDARR(2500),BACKAR(2500),SAVARR(2500)
      COMMON    /TMPARR/  INDARR,BACKAR,SAVARR
C 
      REAL      RA(NDIM),RRA
C 
      DO 5, J=1,NDIM
         INDARR(J) = J
5     CONTINUE
C 
      L = NDIM/2 + 1
      IR = NDIM
C 
10    IF (L.GT.1) THEN            !still hiring
         L = L - 1
         INDXT = INDARR(L)
         RRA = RA(INDXT)
      ELSE                        !in retirement + promotion phase
         INDXT = INDARR(IR)
         RRA = RA(INDXT)          !clear a space at end of array
         INDARR(IR) = INDARR(1)   !retire the top of the heap into it
         IR = IR - 1              !decrease the size of the corporation
         IF (IR.EQ.1) THEN        !done with last promotion
            INDARR(1) = INDXT     !the least competent guy of all...
            GOTO 1000
         ENDIF
      ENDIF
C 
      I = L                 !in hiring as well as in promotion phase
      J = L + L             !we here set up to sift down RRA to its right level
C 
20    IF (J.LE.IR) THEN                    !"Do while J<=IR"
         IF (J.LT.IR) THEN
            IF (RA(INDARR(J)).LT.RA(INDARR(J+1)))
     +         J = J + 1                   !compare to the better underling
         ENDIF
         IF (RRA.LT.RA(INDARR(J))) THEN    !demote RRA
            INDARR(I) = INDARR(J)
            I = J
            J = J + J
         ELSE
            J = IR + 1
         ENDIF
         GOTO 20
      ENDIF
C 
      INDARR(I) = INDXT                 !put RRA into its slot
      GOTO 10
C
1000  DO 4000, J=1,NDIM                 !set up the back index array
         BACKAR(INDARR(J)) = J
4000  CONTINUE
C 
      RINDX = INDARR(NH)
C 
      RETURN
      END

      SUBROUTINE XSORT2(RA,NDIM,NH,JWOFF,NEWRA,XSIZE,NCOUNT,RINDX)
C 
C  RA  = real array of which the elements were used for sorting INDARR
C  NDIM  = full dimension of RA
C  NH = index of median - we only sorted fully until NH then we inserted
C  IOFF = index of RA element in question
C  INDARR = semi sorted index array from last sorting of RA
C 
C 
      IMPLICIT NONE
C 
      INTEGER   NDIM,NH,JWOFF,XSIZE,NCOUNT,RINDX
      INTEGER   IOFF,INDXT,JL,JU,JM,M,N
C 
      INTEGER   INDARR(2500),BACKAR(2500),SAVARR(2500)
      COMMON    /TMPARR/  INDARR,BACKAR,SAVARR
C 
      REAL      RA(NDIM),NEWRA(NCOUNT)
      REAL      NEWR,QR
C 
      IOFF = JWOFF
C 
      DO 5000, N=1,NCOUNT
C 
         NEWR = NEWRA(N) 
         RA(IOFF) = NEWR               !new value
         INDXT = BACKAR(IOFF)
C 
         IF (INDXT.EQ.1) THEN
            IF (NEWR.LE.RA(INDARR(2)))             !it stays very first
     +         GOTO 4900
         ELSE
            IF (NEWR.LE.RA(INDARR(1))) THEN        !it's very first
               DO 1000, M=INDXT,2,-1               !remove old entry
                  INDARR(M) = INDARR(M-1)          !by shifting up
                  BACKAR(INDARR(M)) = M
1000           CONTINUE
               INDARR(1) = IOFF 
               BACKAR(IOFF) = 1
               GOTO 4900
            ENDIF
         ENDIF
C 
         IF (INDXT.EQ.NDIM) THEN
            IF (NEWR.GE.RA(INDARR(NDIM-1)))        !it stays very last
     +         GOTO 4900
         ELSE
            IF (NEWR.GE.RA(INDARR(NDIM))) THEN     !it very last
               DO 2000, M=INDXT,NDIM-1             !remove old entry
                  INDARR(M) = INDARR(M+1)          !by shifting down
                  BACKAR(INDARR(M)) = M
2000           CONTINUE
               INDARR(NDIM) = IOFF
               BACKAR(IOFF) = NDIM
               GOTO 4900
            ENDIF
         ENDIF
C 
         IF (INDXT.EQ.1) THEN
            QR = NEWR - 1.0                        !so we satisy following test
         ELSE
            QR = RA(INDARR(INDXT-1))
         ENDIF
C 
         IF (NEWR.GT.QR) THEN                      !INDXT = 1 also done here...
            IF (NEWR.LE.RA(INDARR(INDXT+1))) GOTO 4900
            JL = INDXT
            JU = NDIM + 1
C
2010        JM = (JU+JL) / 2
            IF (NEWR.LT.RA(INDARR(JM))) THEN
               JU = JM
            ELSE
               JL = JM
            ENDIF
            IF ((JU-JL).GT.1) GOTO 2010          !loop
C
            DO 2100, M=INDXT,JL-1
               INDARR(M) = INDARR(M+1)
               BACKAR(INDARR(M)) = M
2100        CONTINUE
            INDARR(JL) = IOFF
            BACKAR(IOFF) = JL 
C            
         ELSE IF (NEWR.LT.QR) THEN
            JL = 0
            JU = INDXT
C
3010        JM = (JU+JL) / 2
            IF (NEWR.LT.RA(INDARR(JM))) THEN
               JU = JM
            ELSE
               JL = JM
            ENDIF
            IF ((JU-JL).GT.1) GOTO 3010          !loop
C
            DO 3100, M=INDXT,JL+2,-1
               INDARR(M) = INDARR(M-1)
               BACKAR(INDARR(M)) = M
3100        CONTINUE
            INDARR(JL+1) = IOFF
            BACKAR(IOFF) = JL + 1
         ENDIF
C 
4900     IOFF = IOFF + XSIZE
5000  CONTINUE
C 
      RINDX = INDARR(NH)
C
      RETURN
      END

      SUBROUTINE XSAVX(DIRFLG,NDIM)
C 
      INTEGER   DIRFLG, NDIM
      INTEGER   N,TEST
C
      INTEGER   INDARR(2500),BACKAR(2500),SAVARR(2500)
      COMMON    /TMPARR/  INDARR,BACKAR,SAVARR
C 
      IF (DIRFLG.EQ.1) THEN
         DO 100, N=1,NDIM
            TEST = INDARR(N)
            SAVARR(N) = TEST
100      CONTINUE
      ELSE
C 
         DO 500, N=1,NDIM
            TEST = SAVARR(N)
            INDARR(N) = TEST
500      CONTINUE
         DO 600, N=1,NDIM                       !rebuild back index array
            BACKAR(INDARR(N)) = N
600      CONTINUE
      ENDIF
C 
      RETURN
      END