      SUBROUTINE PACCRN(X,IX,IXPAK,NWDS,NBITS,XMIN,XMAX,KIND) 
C 
C     * JUL 07/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)
C     * MAR 22/93 - F.MAJAESS   - DECLARE SOME VARIABLES WITH REAL*8.
C     * JAN 06/93 - F.MAJAESS   - IMPLEMENT CHANGES FOR FIXED FIELD.
C     * MAY 22/92 - A.J. STACEY - RELEASED (PACCRN2 RENAMED TO PACCRN).
C     * MAY 15/92 - A.J. STACEY - BACK OUT INTERPOLATION BETWEEN XMIN
C     *                           AND XMAX.  OPTIMIZE LOOP TEST FOR XMAX
C     *                           CORRECTION.
C     * APR 30/92 - A.J. STACEY - MODIFY PACCRN TO USE "XMIN" AND "XMAX"
C     *                           INSTEAD OF "XMIN" AND "XSCAL" AS THE 
C     *                           CONSTANTS THAT ARE STORED WITH THE 
C     *                           PACKED DATA.
C     * APR 15/92 - A.J. STACEY - PACCRN RENAMED TO PACCRN2.
C     * MAR 23/92 - A.J. STACEY - EXTENSIVE CHANGES TO ALLOW PACKER TO 
C     *                           RUN ON BOTH 32-BIT AND 64-BIT
C     *                           COMPUTERS.
C     * JUN 14/89 - F. MAJAESS (CHANGE GBYTES/SBYTES CALL TO GBYTESB/
C     *                         SBYTESB) 
C     * FEB 06/84 - R.LAPRISE, J.F.FORTIN.
C
C     * FOR KIND.GE.0, PACK NWDS OF FLOATING POINT ARRAY X
C     *                INTO ARRAY IXPAK.
C     * FOR KIND.LT.0, UNPACK IXPAK INTO NWDS OF X. 
C     * PACKING STEPS, (FOR EACH ARRAY ELEMENT) 
C     *    1- FIND THE MINIMUM (XMIN) AND MAXIMUM (XMAX)
C     *       VALUE IN THE ARRAY X. 
C     *    2- REMOVE BIAS (SUBTRACT XMIN FROM X). 
C     *    3- SCALE BY FACTOR (2**NBITS-1)/RANGE SO THAT
C     *       THE VALUES ARE BETWEEN 0. AND 1.*2**NBITS-1.
C     *    4- TRANSFORM INTO INTEGER ARRAY. 
C     *    5- PACK THIS ARRAY WITH SBYTESB. 
C     *
C     * UNPACKING, (FOR EACH ARRAY ELEMENT) 
C     *    1- UNPACK WITH GBYTESC.
C     *    2- TRANSFORM TO REAL ARRAY.
C     *    3- RESCALE FIELD WITH XMIN AND XSCAL.
C     *
C     * PARAMETERS, 
C     *    X     = NON-PACKED ARRAY. (INPUT IF KIND .GE. 0, OUTPUT IF 
C     *                               KIND .LT. 0). 
C     *    IX    = NON-PACKED ARRAY SERVING AS INPUT FOR SBYTESB
C     *            AND OUTPUT FOR GBYTESC, SAME PLACE IN MEMORY THAN X. 
C     *    IXPAK = PACKED ARRAY. (OUTPUT IF KIND .GE. 0, INPUT IF 
C     *                           KIND .LT. 0)
C     *    NWDS  = NUMBER OF ELEMENTS IN ARRAY X (INPUT). 
C     *    NBITS = NUMBER OF BITS TO CODE THE ARRAY ELEMENTS IN MEMORY. 
C     *    XMIN  = MINIMUM OF X ARRAY. (KIND .GE. 0=OUTPUT,ELSE INPUT)
C     *    XSCAL = SCALE USED = (2**NBITS-1)/(XMAX-XMIN). (OUTPUT KIND
C     *                                            .GE. 0, INPUT ELSE)
C     *    KIND  = PACKING (KIND .GE. 0) OR UNPACKING (KIND .LT. 0).
C 
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      DIMENSION X(NWDS)
      DIMENSION IX(NWDS),IXPAK(1) 
      REAL*8  BIGGEST,XSCAL,TEST,XSCALI,BIGSTI
C     INTEGER ISMIN, ISMAX
      COMMON /MACHTYP/ MACHINE,INTSIZE
C
C     * DECLARE AN ARRAY TO IMPLEMENT TRUNCATION OF FLOATING POINT NUMBS.
C
      INTEGER IXSCALI(2),IXTMP(2)
      EQUIVALENCE (XSCALI,IXSCALI),(XTMP,IXTMP)
C 
CCRN  DATA MAXINT/47/ 
      DATA MAXINT/32/
CCRN  DATA SMALL,VLARGE/1.D-13,1.D+38/
      DATA SMALL,VLARGE/1.E-38,1.E+38/
C ------------------------------------------------------------------- 

C     * ARITHMETIC IS NOT GUARANTEED IF NBITS IS TOO LARGE. 
C 
      IF (NBITS .GT. MAXINT) CALL XIT ('PACCRN',-1) 
      IF (NBITS .LE. 1)      CALL XIT ('PACCRN',-2)
C
C     * CALCULATE LARGEST INTEGER THAT CAN BE REPRESENTED BY "NBITS".
C     * THE ORIGINAL CODE WON'T WORK ON 32-BIT MACHINE SINCE 2**NBITS CANNOT BE 
C     * COMPUTED IF NBITS=32 (WHICH IS A VERY LIKELY CHOICE OF 
C     * PACKING FACTOR).  NOTE ALSO THAT ON 32 BIT MACHINES, THE LARGEST
C     * REPRESENTABLE *POSITIVE* INTEGER CANNOT HAVE A ONE IN BIT POSITION 31.
C     * (THIS IS THE SIGN BIT!)
C     * THUS FOR A PACKING FACTOR OF NPACK64=2, ONLY 31 BITS CAN BE USED
C     * TO STORE PACKED DATA.  FOR ANY PACKING FACTOR GREATER THAN 2, ALL THE
C     * BITS DEFINED BY NBITS MAY BE USED TO REPRESENT FLOATING POINT DATA.
C     *
C     * NOTE THAT FOR NPACK=2 "BIGGEST" WILL BE SLIGHTLY DIFFERENT DEPENDING WHETHER
C     * FLOATING POINT NUMBERS ARE 32-BITS OR 64-BITS BECAUSE THE "FLOAT"
C     * OPERATION WILL TRUNCATE SOME OF THE LEAST SIGNIFICANT BITS IN THE 
C     * LARGEST 31-BIT INTEGER.  THIS WILL LEAD TO A SYSTEMATICALLY SMALLER
C     * SET OF DATA IF DATA IS PACKED BY A 64-BIT MACHINE AND ANALYZED IN
C     * 32-BITS.  THIS WILL CAUSE AN ERROR IN (APPROX) THE 7-TH SIGNIFICANT
C     * DIGIT FOR NPACK=2.  OTHER PACKING FACTORS (>2) ARE STILL OK.
C 
CCRN  LARGEST=2**NBITS -1 
      LARGEST=IBITS(-1,0,MIN(31,NBITS))
      BIGGEST=DFLOAT(LARGEST)
C 
C     * TEST FOR PACKING OR UNPACKING CASE. 
C 
      IF (KIND .LT. 0) GO TO 250
C ------------------------------------------------------------------- 
C     * PACKING CASE. 
C 
C     * FIND EXTREMA. 
C 
C     XMIN=X(ISMIN(NWDS,X,1)) 
C     XMAX=X(ISMAX(NWDS,X,1)) 
C
C     * FOLLOWING CODE IS MORE PORTABLE AND SHOULD HAVE SAME PERFORMANCE
C     * IF COMPILER RECOGNIZES THE CONSTRUCT AND SUBSTITUTES THE ABOVE CODE.
C
C      
      XMIN = VLARGE
      XMAX = -VLARGE
      DO 1 I=1,NWDS
        XMIN = MIN(X(I),XMIN)
        XMAX = MAX(X(I),XMAX)
    1 CONTINUE

      RANGE = XMAX - XMIN 
CCRN  IF(ABS(RANGE).LE.1.E-2000) RANGE=1. 
      IF ( ABS(RANGE). LE. SMALL) THEN
        DO 100 I = 1,NWDS*INTSIZE                                     
           IX(I) = 0                                                 
  100   CONTINUE 
        IF ( XMAX.GE.0 .AND. XMIN.GE.0 ) THEN
           XMAX=XMIN
        ELSE
           IF ( XMAX.LE.0 .AND. XMIN.LE.0 ) THEN
             XMIN=XMAX
           ELSE
             IF ( ABS(XMIN).LE.ABS(XMAX) ) THEN
              XMAX=XMIN
             ELSE
              XMIN=XMAX
             ENDIF
           ENDIF
        ENDIF
        CALL SBYTESB(IXPAK,IX,NBITS,NWDS)
        DO 170 I = 1,NWDS
          X(I) = XMIN
  170   CONTINUE 
        RETURN
      ENDIF
C 
C     * SCALE FIELD AND REPRESENT AS INTEGER. 
C     *
C     * A REAL*8 ARRAY "X" AND AN INTEGER ARRAY "IX" SHARE THE SAME MEMORY
C     * SPACE.  SINCE INTEGERS MAY BE 32 BITS LONG AND THE FLOATING POINT
C     * VARIABLES ARE ALWAY REAL*8 HERE, WE ADOPT THE CONVENTION THAT THE
C     * INTEGERS ARE STORED IN THE LOW-ORDER 32-BIT PARCEL OF A 64-BIT WORD
C     * ON BOTH 32 AND 64-BIT MACHINES.
C
      XSCAL = BIGGEST / RANGE 
      IF (INTSIZE .EQ. 1) THEN
        DO 200 I = 1,NWDS 
  200   IX(I)   = INT (XSCAL * (X(I) - XMIN) + 0.5E0) 
      ELSE
        DO 201 I = 1,NWDS 
  201   IX(I*INTSIZE)   = INT (XSCAL * (X(I) - XMIN) + 0.5E0) 
        DO 202 I = 1,2*NWDS,2 
  202   IX(I) = 0
      END IF
C 
C     * PACK INTEGER VALUES INTO NBITS. 
C 
  210 CONTINUE
      CALL SBYTESB (IXPAK,IX,NBITS,NWDS)
C 
C     * RESTORE (APPROXIMATE) FIELD OF X BEFORE LEAVING.
C 
      GO TO 290 
C ------------------------------------------------------------------- 
C     * UNPACKING CASE. 
C 
C     * EXPAND EVERY NBITS OF IXPAK INTO WORD OF IX.
C 
  250 CALL GBYTESB (IXPAK,IX,NBITS,NWDS)
C 
C     * RESCALE X FIELD. (NOTE THAT "IX" MAY BE A 32-BIT INTEGER ARRAY.)
C 
      RANGE = XMAX - XMIN 
CCRN  IF(ABS(RANGE).LE.1.E-2000) RANGE=1. 
      IF ( ABS(RANGE). LE. SMALL) THEN
        IF ( XMAX.GE.0 .AND. XMIN.GE.0 ) THEN
           XMAX=XMIN
        ELSE
           IF ( XMAX.LE.0 .AND. XMIN.LE.0 ) THEN
             XMIN=XMAX
           ELSE
             IF ( ABS(XMIN).LE.ABS(XMAX) ) THEN
              XMAX=XMIN
             ELSE
              XMIN=XMAX
             ENDIF
           ENDIF
        ENDIF
        DO 270 I = 1,NWDS
          X(I) = XMIN
  270   CONTINUE 
        RETURN
      ENDIF
                  
  290 CONTINUE

      TEST  = RANGE/BIGGEST/2.0E0
      XSCALI = RANGE / BIGGEST
CTAPE CALL TRUNC(IXSCALI,1,5)
      BIGSTI = 1.E0/BIGGEST
      IF (INTSIZE .EQ. 1) THEN
        DO 300 I = 1,NWDS 
CBUG
CBUG      * USE DFLOAT INSTEAD OF FLOAT UNTIL MIPS BUG
CBUG      * IS FIXED.
CBUG
          X(I) = DFLOAT (IX(I))         * XSCALI + XMIN
CTAPE     XTMP = DFLOAT(IX(I))*XSCALI
CTAPE     CALL TRUNC(IXTMP,1,5)
CTAPE     XTMP = XMIN + XTMP
CTAPE     CALL TRUNC(IXTMP,1,5)
CTAPE     X(I) = XTMP
C         X(I) = FLOAT (IX(I*INTSIZE)) * XSCALI + XMIN
C         X(I) = XMAX*(FLOAT(IX(I*INTSIZE))*BIGSTI) +
C    1           XMIN*(1.0 - (FLOAT(IX(I*INTSIZE))*BIGSTI)) 
C         R    = FLOAT(IX(I))*BIGSTI
C         X(I) = XMIN - XMIN*R + XMAX*R
C
C         * THE FOLLOWING CODE TWEAKS THE VALUE OF X(I) TO EXACTLY EQUAL
C         * XMAX WHEN X(I)=XMAX WITHIN ROUNDOFF ERROR.  THIS CORRECTS A
C         * SINGLE BIT ERROR ON THE SX-3 DUE TO ROUNDOFF ERROR.
C         * IT IS NOT NECESSARY TO PERFORM THE SAME TWEAKING FOR "XMIN"
C         * IF "X(I) = FLOAT(IX(I))*XSCALI + XMIN" IS USED.
C
          IF (ABS(X(I)-XMAX).LT.TEST) X(I) = XMAX
  300   CONTINUE
      ELSE
        DO 310 I = 1,NWDS 
CBUG
CBUG      * USE DFLOAT INSTEAD OF FLOAT UNTIL MIPS BUG
CBUG      * IS FIXED.
CBUG
          X(I) = DFLOAT (IX(I*INTSIZE)) * XSCALI + XMIN
CTAPE     XTMP = DFLOAT(IX(I*INTSIZE))*XSCALI
CTAPE     CALL TRUNC(IXTMP,1,5)
CTAPE     XTMP = XMIN + XTMP
CTAPE     CALL TRUNC(IXTMP,1,5)
CTAPE     X(I) = XTMP
C         X(I) = XMAX*(FLOAT(IX(I*INTSIZE))*BIGSTI) +
C    1           XMIN*(1.0 - (FLOAT(IX(I*INTSIZE))*BIGSTI)) 
C         R    = FLOAT(IX(I*2))*BIGSTI
C         X(I) = XMIN - XMIN*R + XMAX*R
          IF (ABS(X(I)-XMAX).LT.TEST) X(I) = XMAX
  310   CONTINUE
      END IF
      RETURN
C ------------------------------------------------------------------- 
      END 
