SUBROUTINE W3FI75 (IBITL,ITYPE,ITOSS,FLD,IFLD,IBMAP,IBDSFL,
& NPTS,BDS11,IPFLD,PFLD,LEN,LENBDS,IBERR,PDS,IGDS)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . . .
C SUBPROGRAM: W3FI75 GRIB PACK DATA AND FORM BDS OCTETS(1-11)
C PRGMMR: FARLEY ORG: NMC421 DATE:94-11-22
C
C ABSTRACT: THIS ROUTINE PACKS A GRIB FIELD AND FORMS OCTETS(1-11)
C OF THE BINARY DATA SECTION (BDS).
C
C PROGRAM HISTORY LOG:
C 92-07-10 M. FARLEY ORIGINAL AUTHOR
C 92-10-01 R.E.JONES CORRECTION FOR FIELD OF CONSTANT DATA
C 92-10-16 R.E.JONES GET RID OF ARRAYS FP AND INT
C 93-08-06 CAVANAUGH ADDED ROUTINES FI7501, FI7502, FI7503
C TO ALLOW SECOND ORDER PACKING IN PDS.
C 93-07-21 STACKPOLE ASSORTED REPAIRS TO GET 2ND DIFF PACK IN
C 93-10-28 CAVANAUGH COMMENTED OUT NONOPERATIONAL PRINTS AND
C WRITE STATEMENTS
C 93-12-15 CAVANAUGH CORRECTED LOCATION OF START OF FIRST ORDER
C VALUES AND START OF SECOND ORDER VALUES TO
C REFLECT A BYTE LOCATION IN THE BDS INSTEAD
C OF AN OFFSET IN SUBROUTINE FI7501.
C 94-01-27 CAVANAUGH ADDED IGDS AS INPUT ARGUMENT TO THIS ROUTINE
C AND ADDED PDS AND IGDS ARRAYS TO THE CALL TO
C W3FI82 TO PROVIDE INFORMATION NEEDED FOR
C BOUSTROPHEDONIC PROCESSING.
C 94-05-25 CAVANAUGH SUBROUTINE FI7503 HAS BEEN ADDED TO PROVIDE
C FOR ROW BY ROW OR COLUMN BY COLUMN SECOND
C ORDER PACKING. THIS FEATURE CAN BE ACTIVATED
C BY SETTING IBDSFL(7) TO ZERO.
C 94-07-08 CAVANAUGH COMMENTED OUT PRINT STATEMENTS USED FOR DEBUG
C 94-11-22 FARLEY ENLARGED WORK ARRAYS TO HANDLE .5DEGREE GRIDS
C 95-06-01 R.E.JONES CORRECTION FOR NUMBER OF UNUSED BITS AT END
C OF SECTION 4, IN BDS BYTE 4, BITS 5-8.
C 95-10-31 IREDELL REMOVED SAVES AND PRINTS
C 2001-06-06 GILBERT CHanged gbyte/sbyte calls to refer to
C Wesley Ebisuzaki's endian independent
C versions gbytec/sbytec.
C Use f90 standard routine bit_size to get
C number of bits in an integer instead of w3fi01.
C
C USAGE: CALL W3FI75 (IBITL,ITYPE,ITOSS,FLD,IFLD,IBMAP,IBDSFL,
C & NPTS,BDS11,IPFLD,PFLD,LEN,LENBDS,IBERR,PDS,IGDS)
C INPUT ARGUMENT LIST:
C IBITL - 0, COMPUTER COMPUTES PACKING LENGTH FROM POWER
C OF 2 THAT BEST FITS THE DATA.
C 8, 12, ETC. COMPUTER RESCALES DATA TO FIT INTO
C SET NUMBER OF BITS.
C ITYPE - 0 = IF INPUT DATA IS FLOATING POINT (FLD)
C 1 = IF INPUT DATA IS INTEGER (IFLD)
C ITOSS - 0 = NO BIT MAP IS INCLUDED (DON'T TOSS DATA)
C 1 = TOSS NULL DATA ACCORDING TO IBMAP
C FLD - REAL ARRAY OF DATA TO BE PACKED IF ITYPE=0
C IFLD - INTEGER ARRAY TO BE PACKED IF ITYPE=1
C IBMAP - BIT MAP SUPPLIED FROM USER
C IBDSFL - INTEGER ARRAY CONTAINING TABLE 11 FLAG INFO
C BDS OCTET 4:
C (1) 0 = GRID POINT DATA
C 1 = SPHERICAL HARMONIC COEFFICIENTS
C (2) 0 = SIMPLE PACKING
C 1 = SECOND ORDER PACKING
C (3) 0 = ORIGINAL DATA WERE FLOATING POINT VALUES
C 1 = ORIGINAL DATA WERE INTEGER VALUES
C (4) 0 = NO ADDITIONAL FLAGS AT OCTET 14
C 1 = OCTET 14 CONTAINS FLAG BITS 5-12
C (5) 0 = RESERVED - ALWAYS SET TO 0
C (6) 0 = SINGLE DATUM AT EACH GRID POINT
C 1 = MATRIX OF VALUES AT EACH GRID POINT
C (7) 0 = NO SECONDARY BIT MAPS
C 1 = SECONDARY BIT MAPS PRESENT
C (8) 0 = SECOND ORDER VALUES HAVE CONSTANT WIDTH
C 1 = SECOND ORDER VALUES HAVE DIFFERENT WIDTHS
C NPTS - NUMBER OF GRIDPOINTS IN ARRAY TO BE PACKED
C IGDS - ARRAY OF GDS INFORMATION
C
C OUTPUT ARGUMENT LIST:
C BDS11 - FIRST 11 OCTETS OF BDS
C PFLD - PACKED GRIB FIELD
C LEN - LENGTH OF PFLD
C LENBDS - LENGTH OF BDS
C IBERR - 1, ERROR CONVERTING IEEE F.P. NUMBER TO IBM370 F.P.
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C LANGUAGE: IBM VS FORTRAN 77, CRAY CFT77 FORTRAN
C MACHINE: HDS, CRAY C916/256, Y-MP8/64, Y-MP EL92/256
C
C$$$
C
REAL FLD(*)
C REAL FWORK(260000)
C
C FWORK CAN USE DYNAMIC ALLOCATION OF MEMORY ON CRAY
C
REAL FWORK(NPTS)
REAL RMIN,REFNCE
C
character(len=1) IPFLD(*)
INTEGER IBDSFL(*)
INTEGER IBMAP(*)
INTEGER IFLD(*),IGDS(*)
C INTEGER IWORK(260000)
C
C IWORK CAN USE DYNAMIC ALLOCATION OF MEMORY ON CRAY
C
INTEGER IWORK(NPTS)
C
LOGICAL CONST
C
CHARACTER * 1 BDS11(11),PDS(*)
CHARACTER * 1 PFLD(*)
C
C 1.0 PACK THE FIELD.
C
C 1.1 TOSS DATA IF BITMAP BEING USED,
C MOVING 'DATA' TO WORK AREA...
C
CONST = .FALSE.
IBERR = 0
IW = 0
C
IF (ITOSS .EQ. 1) THEN
IF (ITYPE .EQ. 0) THEN
DO 110 IT=1,NPTS
IF (IBMAP(IT) .EQ. 1) THEN
IW = IW + 1
FWORK(IW) = FLD(IT)
ENDIF
110 CONTINUE
NPTS = IW
ELSE IF (ITYPE .EQ. 1) THEN
DO 111 IT=1,NPTS
IF (IBMAP(IT) .EQ. 1) THEN
IW = IW + 1
IWORK(IW) = IFLD(IT)
ENDIF
111 CONTINUE
NPTS = IW
ENDIF
C
C ELSE, JUST MOVE DATA TO WORK ARRAY
C
ELSE IF (ITOSS .EQ. 0) THEN
IF (ITYPE .EQ. 0) THEN
DO 112 IT=1,NPTS
FWORK(IT) = FLD(IT)
112 CONTINUE
ELSE IF (ITYPE .EQ. 1) THEN
DO 113 IT=1,NPTS
IWORK(IT) = IFLD(IT)
113 CONTINUE
ENDIF
ENDIF
C
C 1.2 CONVERT DATA IF NEEDED PRIOR TO PACKING.
C (INTEGER TO F.P. OR F.P. TO INTEGER)
C ITYPE = 0...FLOATING POINT DATA
C IBITL = 0...PACK IN LEAST # BITS...CONVERT TO INTEGER
C ITYPE = 1...INTEGER DATA
C IBITL > 0...PACK IN FIXED # BITS...CONVERT TO FLOATING POINT
C
IF (ITYPE .EQ. 0 .AND. IBITL .EQ. 0) THEN
DO 120 IF=1,NPTS
IWORK(IF) = NINT(FWORK(IF))
120 CONTINUE
ELSE IF (ITYPE .EQ. 1 .AND. IBITL .NE. 0) THEN
DO 123 IF=1,NPTS
FWORK(IF) = FLOAT(IWORK(IF))
123 CONTINUE
ENDIF
C
C 1.3 PACK THE DATA.
C
IF (IBDSFL(2).NE.0) THEN
C SECOND ORDER PACKING
C
C PRINT*,' DOING SECOND ORDER PACKING...'
IF (IBITL.EQ.0) THEN
C
C PRINT*,' AND VARIABLE BIT PACKING'
C
C WORKING WITH INTEGER VALUES
C SINCE DOING VARIABLE BIT PACKING
C
MAX = IWORK(1)
MIN = IWORK(1)
DO 300 I = 2, NPTS
IF (IWORK(I).LT.MIN) THEN
MIN = IWORK(I)
ELSE IF (IWORK(I).GT.MAX) THEN
MAX = IWORK(I)
END IF
300 CONTINUE
C EXTRACT MINIMA
DO 400 I = 1, NPTS
C IF (IWORK(I).LT.0) THEN
C PRINT *,'MINIMA 400',I,IWORK(I),NPTS
C END IF
IWORK(I) = IWORK(I) - MIN
400 CONTINUE
REFNCE = MIN
IDIFF = MAX - MIN
C PRINT *,'REFERENCE VALUE',REFNCE
C
C WRITE (6,FMT='('' MINIMA REMOVED = '',/,
C & 10(3X,10I10,/))') (IWORK(I),I=1,6)
C WRITE (6,FMT='('' END OF ARRAY = '',/,
C & 10(3X,10I10,/))') (IWORK(I),I=NPTS-5,NPTS)
C
C FIND BIT WIDTH OF IDIFF
C
CALL FI7505 (IDIFF,KWIDE)
C PRINT*,' BIT WIDTH FOR ORIGINAL DATA', KWIDE
ISCAL2 = 0
C
C MULTIPLICATIVE SCALE FACTOR SET TO 1
C IN ANTICIPATION OF POSSIBLE USE IN GLAHN 2DN DIFF
C
SCAL2 = 1.
C
ELSE
C
C PRINT*,' AND FIXED BIT PACKING, IBITL = ', IBITL
C FIXED BIT PACKING
C - LENGTH OF FIELD IN IBITL
C - MUST BE REAL DATA
C FLOATING POINT INPUT
C
RMAX = FWORK(1)
RMIN = FWORK(1)
DO 100 I = 2, NPTS
IF (FWORK(I).LT.RMIN) THEN
RMIN = FWORK(I)
ELSE IF (FWORK(I).GT.RMAX) THEN
RMAX = FWORK(I)
END IF
100 CONTINUE
REFNCE = RMIN
C PRINT *,'100 REFERENCE',REFNCE
C EXTRACT MINIMA
DO 200 I = 1, NPTS
FWORK(I) = FWORK(I) - RMIN
200 CONTINUE
C PRINT *,'REFERENCE VALUE',REFNCE
C WRITE (6,FMT='('' MINIMA REMOVED = '',/,
C & 10(3X,10F8.2,/))') (FWORK(I),I=1,6)
C WRITE (6,FMT='('' END OF ARRAY = '',/,
C & 10(3X,10F8.2,/))') (FWORK(I),I=NPTS-5,NPTS)
C FIND LARGEST DELTA
IDELT = NINT(RMAX - RMIN)
C DO BINARY SCALING
C FIND OUT WHAT BINARY SCALE FACTOR
C PERMITS CONTAINMENT OF
C LARGEST DELTA
CALL FI7505 (IDELT,IWIDE)
C
C BINARY SCALING
C
ISCAL2 = IWIDE - IBITL
C PRINT *,'SCALING NEEDED TO FIT =',ISCAL2
C PRINT*,' RANGE OF = ',IDELT
C
C EXPAND DATA WITH BINARY SCALING
C CONVERT TO INTEGER
SCAL2 = 2.0**ISCAL2
SCAL2 = 1./ SCAL2
DO 600 I = 1, NPTS
IWORK(I) = NINT(FWORK(I) * SCAL2)
600 CONTINUE
KWIDE = IBITL
END IF
C
C *****************************************************************
C
C FOLLOWING IS FOR GLAHN SECOND DIFFERENCING
C NOT STANDARD GRIB
C
C TEST FOR SECOND DIFFERENCE PACKING
C BASED OF SIZE OF PDS - SIZE IN FIRST 3 BYTES
C
CALL GBYTEC(PDS,IPDSIZ,0,24)
IF (IPDSIZ.EQ.50) THEN
C PRINT*,' DO SECOND DIFFERENCE PACKING '
C
C GLAHN PACKING TO 2ND DIFFS
C
C WRITE (6,FMT='('' CALL TO W3FI82 WITH = '',/,
C & 10(3X,10I6,/))') (IWORK(I),I=1,NPTS)
C
CALL W3FI82 (IWORK,FVAL1,FDIFF1,NPTS,PDS,IGDS)
C
C PRINT *,'GLAHN',FVAL1,FDIFF1
C WRITE (6,FMT='('' OUT FROM W3FI82 WITH = '',/,
C & 10(3X,10I6,/))') (IWORK(I),I=1,NPTS)
C
C MUST NOW RE-REMOVE THE MINIMUM VALUE
C OF THE SECOND DIFFERENCES TO ASSURE
C ALL POSITIVE NUMBERS FOR SECOND ORDER GRIB PACKING
C
C ORIGINAL REFERENCE VALUE ADDED TO FIRST POINT
C VALUE FROM THE 2ND DIFF PACKER TO BE ADDED
C BACK IN WHEN THE 2ND DIFF VALUES ARE
C RECONSTRUCTED BACK TO THE BASIC VALUES
C
C ALSO, THE REFERENCE VALUE IS
C POWER-OF-TWO SCALED TO MATCH
C FVAL1. ALL OF THIS SCALING
C WILL BE REMOVED AFTER THE
C GLAHN SECOND DIFFERENCING IS UNDONE.
C THE SCALING FACTOR NEEDED TO DO THAT
C IS SAVED IN THE PDS AS A SIGNED POSITIVE
C TWO BYTE INTEGER
C
C THE SCALING FOR THE 2ND DIF PACKED
C VALUES IS PROPERLY SET TO ZERO
C
FVAL1 = FVAL1 + REFNCE*SCAL2
C FIRST TEST TO SEE IF
C ON 32 OR 64 BIT COMPUTER
C CALL W3FI01(LW)
IF (bit_size(LW).EQ.32) THEN
CALL W3FI76 (FVAL1,IEXP,IMANT,32)
ELSE
CALL W3FI76 (FVAL1,IEXP,IMANT,64)
END IF
CALL SBYTEC(PDS,IEXP,320,8)
CALL SBYTEC(PDS,IMANT,328,24)
C
IF (bit_size(LW).EQ.32) THEN
CALL W3FI76 (FDIFF1,IEXP,IMANT,32)
ELSE
CALL W3FI76 (FDIFF1,IEXP,IMANT,64)
END IF
CALL SBYTEC(PDS,IEXP,352,8)
CALL SBYTEC(PDS,IMANT,360,24)
C
C TURN ISCAL2 INTO SIGNED POSITIVE INTEGER
C AND STORE IN TWO BYTES
C
IF(ISCAL2.GE.0) THEN
CALL SBYTEC(PDS,ISCAL2,384,16)
ELSE
CALL SBYTEC(PDS,1,384,1)
ISCAL2 = - ISCAL2
CALL SBYTEC( PDS,ISCAL2,385,15)
ENDIF
C
MAX = IWORK(1)
MIN = IWORK(1)
DO 700 I = 2, NPTS
IF (IWORK(I).LT.MIN) THEN
MIN = IWORK(I)
ELSE IF (IWORK(I).GT.MAX) THEN
MAX = IWORK(I)
END IF
700 CONTINUE
C EXTRACT MINIMA
DO 710 I = 1, NPTS
IWORK(I) = IWORK(I) - MIN
710 CONTINUE
REFNCE = MIN
C PRINT *,'710 REFERENCE',REFNCE
ISCAL2 = 0
C
C AND RESET VALUE OF KWIDE - THE BIT WIDTH
C FOR THE RANGE OF THE VALUES
C
IDIFF = MAX - MIN
CALL FI7505 (IDIFF,KWIDE)
C
C PRINT*,'BIT WIDTH (KWIDE) OF 2ND DIFFS', KWIDE
C
C **************************** END OF GLAHN PACKING ************
ELSE IF (IBDSFL(2).EQ.1.AND.IBDSFL(7).EQ.0) THEN
C HAVE SECOND ORDER PACKING WITH NO SECOND ORDER
C BIT MAP. ERGO ROW BY ROW - COL BY COL
CALL FI7503 (IWORK,IPFLD,NPTS,IBDSFL,BDS11,
* LEN,LENBDS,PDS,REFNCE,ISCAL2,KWIDE,IGDS)
RETURN
END IF
C WRITE (6,FMT='('' CALL TO FI7501 WITH = '',/,
C & 10(3X,10I6,/))') (IWORK(I),I=1,NPTS)
C WRITE (6,FMT='('' END OF ARRAY = '',/,
C & 10(3X,10I6,/))') (IWORK(I),I=NPTS-5,NPTS)
C PRINT*,' REFNCE,ISCAL2, KWIDE AT CALL TO FI7501',
C & REFNCE, ISCAL2,KWIDE
C
C SECOND ORDER PACKING
C
CALL FI7501 (IWORK,IPFLD,NPTS,IBDSFL,BDS11,
* LEN,LENBDS,PDS,REFNCE,ISCAL2,KWIDE)
C
C BDS COMPLETELY ASSEMBLED IN FI7501 FOR SECOND ORDER
C PACKING.
C
ELSE
C SIMPLE PACKING
C
C PRINT*,' SIMPLE FIRST ORDER PACKING...'
IF (IBITL.EQ.0) THEN
C PRINT*,' WITH VARIABLE BIT LENGTH'
C
C WITH VARIABLE BIT LENGTH, ADJUSTED
C TO ACCOMMODATE LARGEST VALUE
C BINARY SCALING ALWAYS = 0
C
CALL W3FI58(IWORK,NPTS,IWORK,PFLD,NBITS,LEN,KMIN)
RMIN = KMIN
REFNCE = RMIN
ISCALE = 0
C PRINT*,' BIT LENGTH CAME OUT AT ...',NBITS
C
C SET CONST .TRUE. IF ALL VALUES ARE THE SAME
C
IF (LEN.EQ.0.AND.NBITS.EQ.0) CONST = .TRUE.
C
ELSE
C PRINT*,' FIXED BIT LENGTH, IBITL = ', IBITL
C
C FIXED BIT LENGTH PACKING (VARIABLE PRECISION)
C VALUES SCALED BY POWER OF 2 (ISCALE) TO
C FIT LARGEST VALUE INTO GIVEN BIT LENGTH (IBITL)
C
CALL W3FI59(FWORK,NPTS,IBITL,IWORK,PFLD,ISCALE,LEN,RMIN)
REFNCE = RMIN
C PRINT *,' SCALING NEEDED TO FIT IS ...', ISCALE
NBITS = IBITL
C
C SET CONST .TRUE. IF ALL VALUES ARE THE SAME
C
IF (LEN.EQ.0) THEN
CONST = .TRUE.
NBITS = 0
END IF
END IF
C
C$ COMPUTE LENGTH OF BDS IN OCTETS
C
INUM = NPTS * NBITS + 88
C PRINT *,'NUMBER OF BITS BEFORE FILL ADDED',INUM
C
C NUMBER OF FILL BITS
NFILL = 0
NLEFT = MOD(INUM,16)
IF (NLEFT.NE.0) THEN
INUM = INUM + 16 - NLEFT
NFILL = 16 - NLEFT
END IF
C PRINT *,'NUMBER OF BITS AFTER FILL ADDED',INUM
C LENGTH OF BDS IN BYTES
LENBDS = INUM / 8
C
C 2.0 FORM THE BINARY DATA SECTION (BDS).
C
C CONCANTENATE ALL FIELDS FOR BDS
C
C BYTES 1-3
CALL SBYTEC (BDS11,LENBDS,0,24)
C
C BYTE 4
C FLAGS
CALL SBYTEC (BDS11,IBDSFL(1),24,1)
CALL SBYTEC (BDS11,IBDSFL(2),25,1)
CALL SBYTEC (BDS11,IBDSFL(3),26,1)
CALL SBYTEC (BDS11,IBDSFL(4),27,1)
C NR OF FILL BITS
CALL SBYTEC (BDS11,NFILL,28,4)
C
C$ FILL OCTETS 5-6 WITH THE SCALE FACTOR.
C
C BYTE 5-6
IF (ISCALE.LT.0) THEN
CALL SBYTEC (BDS11,1,32,1)
ISCALE = - ISCALE
CALL SBYTEC (BDS11,ISCALE,33,15)
ELSE
CALL SBYTEC (BDS11,ISCALE,32,16)
END IF
C
C$ FILL OCTET 7-10 WITH THE REFERENCE VALUE
C CONVERT THE FLOATING POINT OF YOUR MACHINE TO IBM370 32 BIT
C FLOATING POINT NUMBER
C
C BYTE 7-10
C REFERENCE VALUE
C FIRST TEST TO SEE IF
C ON 32 OR 64 BIT COMPUTER
C CALL W3FI01(LW)
IF (bit_size(LW).EQ.32) THEN
CALL W3FI76 (REFNCE,IEXP,IMANT,32)
ELSE
CALL W3FI76 (REFNCE,IEXP,IMANT,64)
END IF
CALL SBYTEC (BDS11,IEXP,48,8)
CALL SBYTEC (BDS11,IMANT,56,24)
C
C
C$ FILL OCTET 11 WITH THE NUMBER OF BITS.
C
C BYTE 11
CALL SBYTEC (BDS11,NBITS,80,8)
END IF
C
RETURN
END
SUBROUTINE FI7501 (IWORK,IPFLD,NPTS,IBDSFL,BDS11,
* LEN,LENBDS,PDS,REFNCE,ISCAL2,KWIDE)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . . .
C SUBPROGRAM: FI7501 BDS SECOND ORDER PACKING
C PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 93-08-06
C
C ABSTRACT: PERFORM SECONDARY PACKING ON GRID POINT DATA,
C GENERATING ALL BDS INFORMATION.
C
C PROGRAM HISTORY LOG:
C 93-08-06 CAVANAUGH
C 93-12-15 CAVANAUGH CORRECTED LOCATION OF START OF FIRST ORDER
C VALUES AND START OF SECOND ORDER VALUES TO
C REFLECT A BYTE LOCATION IN THE BDS INSTEAD
C OF AN OFFSET.
C 95-10-31 IREDELL REMOVED SAVES AND PRINTS
C
C USAGE: CALL FI7501 (IWORK,IPFLD,NPTS,IBDSFL,BDS11,
C * LEN,LENBDS,PDS,REFNCE,ISCAL2,KWIDE)
C INPUT ARGUMENT LIST:
C IWORK - INTEGER SOURCE ARRAY
C NPTS - NUMBER OF POINTS IN IWORK
C IBDSFL - FLAGS
C
C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS)
C IPFLD - CONTAINS BDS FROM BYTE 12 ON
C BDS11 - CONTAINS FIRST 11 BYTES FOR BDS
C LEN - NUMBER OF BYTES FROM 12 ON
C LENBDS - TOTAL LENGTH OF BDS
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C LANGUAGE: IBM VS FORTRAN 77, CRAY CFT77 FORTRAN
C MACHINE: HDS, CRAY C916/256, Y-MP8/64, Y-MP EL92/256
C
C$$$
CHARACTER*1 BDS11(*),PDS(*)
C
REAL REFNCE
C
INTEGER ISCAL2,KWIDE
INTEGER LENBDS
CHARACTER(len=1) IPFLD(*)
INTEGER LEN,KBDS(22)
INTEGER IWORK(*)
C OCTET NUMBER IN SECTION, FIRST ORDER PACKING
C INTEGER KBDS(12)
C FLAGS
INTEGER IBDSFL(*)
C EXTENDED FLAGS
C INTEGER KBDS(14)
C OCTET NUMBER FOR SECOND ORDER PACKING
C INTEGER KBDS(15)
C NUMBER OF FIRST ORDER VALUES
C INTEGER KBDS(17)
C NUMBER OF SECOND ORDER PACKED VALUES
C INTEGER KBDS(19)
C WIDTH OF SECOND ORDER PACKING
character(len=1) ISOWID(400000)
C SECONDARY BIT MAP
character(len=1) ISOBMP(65600)
C FIRST ORDER PACKED VALUES
character(len=1) IFOVAL(400000)
C SECOND ORDER PACKED VALUES
character(len=1) ISOVAL(800000)
C
C INTEGER KBDS(11)
C BIT WIDTH TABLE
INTEGER IBITS(31)
C
DATA IBITS/1,3,7,15,31,63,127,255,511,1023,
* 2047,4095,8191,16383,32767,65535,131072,
* 262143,524287,1048575,2097151,4194303,
* 8388607,16777215,33554431,67108863,
* 134217727,268435455,536870911,
* 1073741823,2147483647/
C ----------------------------------
C INITIALIZE ARRAYS
DO I = 1, 400000
IFOVAL(I) = char(0)
ISOWID(I) = char(0)
ENDDO
C
DO 101 I = 1, 65600
ISOBMP(I) = char(0)
101 CONTINUE
DO 102 I = 1, 800000
ISOVAL(I) = char(0)
102 CONTINUE
C INITIALIZE POINTERS
C SECONDARY BIT WIDTH POINTER
IWDPTR = 0
C SECONDARY BIT MAP POINTER
IBMP2P = 0
C FIRST ORDER VALUE POINTER
IFOPTR = 0
C BYTE POINTER TO START OF 1ST ORDER VALUES
KBDS(12) = 0
C BYTE POINTER TO START OF 2ND ORDER VALUES
KBDS(15) = 0
C TO CONTAIN NUMBER OF FIRST ORDER VALUES
KBDS(17) = 0
C TO CONTAIN NUMBER OF SECOND ORDER VALUES
KBDS(19) = 0
C SECOND ORDER PACKED VALUE POINTER
ISOPTR = 0
C =======================================================
C
C DATA IS IN IWORK
C
KBDS(11) = KWIDE
C
C DATA PACKING
C
ITER = 0
INEXT = 1
ISTART = 1
C -----------------------------------------------------------
KOUNT = 0
C DO 1 I = 1, NPTS, 10
C PRINT *,I,(IWORK(K),K=I, I+9)
C 1 CONTINUE
2000 CONTINUE
ITER = ITER + 1
C PRINT *,'NEXT ITERATION STARTS AT',ISTART
IF (ISTART.GT.NPTS) THEN
GO TO 4000
ELSE IF (ISTART.EQ.NPTS) THEN
KPTS = 1
MXDIFF = 0
GO TO 2200
END IF
C
C LOOK FOR REPITITIONS OF A SINGLE VALUE
CALL FI7502 (IWORK,ISTART,NPTS,ISAME)
IF (ISAME.GE.15) THEN
KOUNT = KOUNT + 1
C PRINT *,'FI7501 - FOUND IDENTICAL SET OF ',ISAME
MXDIFF = 0
KPTS = ISAME
ELSE
C
C LOOK FOR SETS OF VALUES IN TREND SELECTED RANGE
CALL FI7513 (IWORK,ISTART,NPTS,NMAX,NMIN,INRNGE)
C PRINT *,'ISTART ',ISTART,' INRNGE',INRNGE,NMAX,NMIN
IEND = ISTART + INRNGE - 1
C DO 2199 NM = ISTART, IEND, 10
C PRINT *,' ',(IWORK(NM+JK),JK=0,9)
C2199 CONTINUE
MXDIFF = NMAX - NMIN
KPTS = INRNGE
END IF
2200 CONTINUE
C PRINT *,' RANGE ',MXDIFF,' MAX',NMAX,' MIN',NMIN
C INCREMENT NUMBER OF FIRST ORDER VALUES
KBDS(17) = KBDS(17) + 1
C ENTER FIRST ORDER VALUE
IF (MXDIFF.GT.0) THEN
DO 2220 LK = 0, KPTS-1
IWORK(ISTART+LK) = IWORK(ISTART+LK) - NMIN
2220 CONTINUE
CALL SBYTEC (IFOVAL,NMIN,IFOPTR,KBDS(11))
ELSE
CALL SBYTEC (IFOVAL,IWORK(ISTART),IFOPTR,KBDS(11))
END IF
IFOPTR = IFOPTR + KBDS(11)
C PROCESS SECOND ORDER BIT WIDTH
IF (MXDIFF.GT.0) THEN
DO 2330 KWIDE = 1, 31
IF (MXDIFF.LE.IBITS(KWIDE)) THEN
GO TO 2331
END IF
2330 CONTINUE
2331 CONTINUE
ELSE
KWIDE = 0
END IF
CALL SBYTEC (ISOWID,KWIDE,IWDPTR,8)
IWDPTR = IWDPTR + 8
C PRINT *,KWIDE,' IFOVAL=',NMIN,IWORK(ISTART),KPTS
C IF KWIDE NE 0, SAVE SECOND ORDER VALUE
IF (KWIDE.GT.0) THEN
CALL SBYTESC (ISOVAL,IWORK(ISTART),ISOPTR,KWIDE,0,KPTS)
ISOPTR = ISOPTR + KPTS * KWIDE
KBDS(19) = KBDS(19) + KPTS
C PRINT *,' SECOND ORDER VALUES'
C PRINT *,(IWORK(ISTART+I),I=0,KPTS-1)
END IF
C ADD TO SECOND ORDER BITMAP
CALL SBYTEC (ISOBMP,1,IBMP2P,1)
IBMP2P = IBMP2P + KPTS
ISTART = ISTART + KPTS
GO TO 2000
C --------------------------------------------------------------
4000 CONTINUE
C PRINT *,'THERE WERE ',ITER,' SECOND ORDER GROUPS'
C PRINT *,'THERE WERE ',KOUNT,' STRINGS OF CONSTANTS'
C CONCANTENATE ALL FIELDS FOR BDS
C
C REMAINDER GOES INTO IPFLD
IPTR = 0
C BYTES 12-13
C VALUE FOR N1
C LEAVE SPACE FOR THIS
IPTR = IPTR + 16
C BYTE 14
C EXTENDED FLAGS
CALL SBYTEC (IPFLD,IBDSFL(5),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(6),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(7),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(8),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(9),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(10),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(11),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(12),IPTR,1)
IPTR = IPTR + 1
C BYTES 15-16
C SKIP OVER VALUE FOR N2
IPTR = IPTR + 16
C BYTES 17-18
C P1
CALL SBYTEC (IPFLD,KBDS(17),IPTR,16)
IPTR = IPTR + 16
C BYTES 19-20
C P2
CALL SBYTEC (IPFLD,KBDS(19),IPTR,16)
IPTR = IPTR + 16
C BYTE 21 - RESERVED LOCATION
CALL SBYTEC (IPFLD,0,IPTR,8)
IPTR = IPTR + 8
C BYTES 22 - ?
C WIDTHS OF SECOND ORDER PACKING
IX = (IWDPTR + 32) / 32
C CALL SBYTESC (IPFLD,ISOWID,IPTR,32,0,IX)
ijk=IWDPTR/8
jst=(iptr/8)+1
ipfld(jst:jst+ijk)=ISOWID(1:ijk)
IPTR = IPTR + IWDPTR
C SECONDARY BIT MAP
IJ = (IBMP2P + 32) / 32
C CALL SBYTESC (IPFLD,ISOBMP,IPTR,32,0,IJ)
ijk=(IBMP2P/8)+1
jst=(iptr/8)+1
ipfld(jst:jst+ijk)=ISOBMP(1:ijk)
IPTR = IPTR + IBMP2P
IF (MOD(IPTR,8).NE.0) THEN
IPTR = IPTR + 8 - MOD(IPTR,8)
END IF
C DETERMINE LOCATION FOR START
C OF FIRST ORDER PACKED VALUES
KBDS(12) = IPTR / 8 + 12
C STORE LOCATION
CALL SBYTEC (IPFLD,KBDS(12),0,16)
C MOVE IN FIRST ORDER PACKED VALUES
IPASS = (IFOPTR + 32) / 32
C CALL SBYTESC (IPFLD,IFOVAL,IPTR,32,0,IPASS)
ijk=(IFOPTR/8)+1
jst=(iptr/8)+1
ipfld(jst:jst+ijk)=ifoval(1:ijk)
IPTR = IPTR + IFOPTR
IF (MOD(IPTR,8).NE.0) THEN
IPTR = IPTR + 8 - MOD(IPTR,8)
END IF
C PRINT *,'IFOPTR =',IFOPTR,' ISOPTR =',ISOPTR
C DETERMINE LOCATION FOR START
C OF SECOND ORDER VALUES
KBDS(15) = IPTR / 8 + 12
C SAVE LOCATION OF SECOND ORDER VALUES
CALL SBYTEC (IPFLD,KBDS(15),24,16)
C MOVE IN SECOND ORDER PACKED VALUES
IX = (ISOPTR + 32) / 32
c CALL SBYTESC (IPFLD,ISOVAL,IPTR,32,0,IX)
ijk=(ISOPTR/8)+1
jst=(iptr/8)+1
ipfld(jst:jst+ijk)=isoval(1:ijk)
IPTR = IPTR + ISOPTR
NLEFT = MOD(IPTR+88,16)
IF (NLEFT.NE.0) THEN
NLEFT = 16 - NLEFT
IPTR = IPTR + NLEFT
END IF
C COMPUTE LENGTH OF DATA PORTION
LEN = IPTR / 8
C COMPUTE LENGTH OF BDS
LENBDS = LEN + 11
C -----------------------------------
C BYTES 1-3
C THIS FUNCTION COMPLETED BELOW
C WHEN LENGTH OF BDS IS KNOWN
CALL SBYTEC (BDS11,LENBDS,0,24)
C BYTE 4
CALL SBYTEC (BDS11,IBDSFL(1),24,1)
CALL SBYTEC (BDS11,IBDSFL(2),25,1)
CALL SBYTEC (BDS11,IBDSFL(3),26,1)
CALL SBYTEC (BDS11,IBDSFL(4),27,1)
C ENTER NUMBER OF FILL BITS
CALL SBYTEC (BDS11,NLEFT,28,4)
C BYTE 5-6
IF (ISCAL2.LT.0) THEN
CALL SBYTEC (BDS11,1,32,1)
ISCAL2 = - ISCAL2
ELSE
CALL SBYTEC (BDS11,0,32,1)
END IF
CALL SBYTEC (BDS11,ISCAL2,33,15)
C
C$ FILL OCTET 7-10 WITH THE REFERENCE VALUE
C CONVERT THE FLOATING POINT OF YOUR MACHINE TO IBM370 32 BIT
C FLOATING POINT NUMBER
C REFERENCE VALUE
C FIRST TEST TO SEE IF
C ON 32 OR 64 BIT COMPUTER
C CALL W3FI01(LW)
IF (bit_size(LW).EQ.32) THEN
CALL W3FI76 (REFNCE,IEXP,IMANT,32)
ELSE
CALL W3FI76 (REFNCE,IEXP,IMANT,64)
END IF
CALL SBYTEC (BDS11,IEXP,48,8)
CALL SBYTEC (BDS11,IMANT,56,24)
C
C BYTE 11
C
CALL SBYTEC (BDS11,KBDS(11),80,8)
C
RETURN
END
SUBROUTINE FI7502 (IWORK,ISTART,NPTS,ISAME)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . . .
C SUBPROGRAM: FI7502 SECOND ORDER SAME VALUE COLLECTION
C PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 93-06-23
C
C ABSTRACT: COLLECT SEQUENTIAL SAME VALUES FOR PROCESSING
C AS SECOND ORDER VALUE FOR GRIB MESSAGES.
C
C PROGRAM HISTORY LOG:
C 93-06-23 CAVANAUGH
C 95-10-31 IREDELL REMOVED SAVES AND PRINTS
C
C USAGE: CALL FI7502 (IWORK,ISTART,NPTS,ISAME)
C INPUT ARGUMENT LIST:
C IWORK - ARRAY CONTAINING SOURCE DATA
C ISTART - STARTING LOCATION FOR THIS TEST
C NPTS - NUMBER OF POINTS IN IWORK
C
C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS)
C ISAME - NUMBER OF SEQUENTIAL POINTS HAVING THE SAME VALUE
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C LANGUAGE: IBM VS FORTRAN 77, CRAY CFT77 FORTRAN
C MACHINE: HDS, CRAY C916/256, Y-MP8/64, Y-MP EL92/256
C
C$$$
INTEGER IWORK(*)
INTEGER ISTART
INTEGER ISAME
INTEGER K
INTEGER NPTS
C -------------------------------------------------------------
ISAME = 0
DO 100 K = ISTART, NPTS
IF (IWORK(K).NE.IWORK(ISTART)) THEN
RETURN
END IF
ISAME = ISAME + 1
100 CONTINUE
RETURN
END
SUBROUTINE FI7503 (IWORK,IPFLD,NPTS,IBDSFL,BDS11,
* LEN,LENBDS,PDS,REFNCE,ISCAL2,KWIDE,IGDS)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . . .
C SUBPROGRAM: FI7501 ROW BY ROW, COL BY COL PACKING
C PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 94-05-20
C
C ABSTRACT: PERFORM ROW BY ROW OR COLUMN BY COLUMN PACKING
C GENERATING ALL BDS INFORMATION.
C
C PROGRAM HISTORY LOG:
C 93-08-06 CAVANAUGH
C 95-10-31 IREDELL REMOVED SAVES AND PRINTS
C
C USAGE: CALL FI7503 (IWORK,IPFLD,NPTS,IBDSFL,BDS11,
C * LEN,LENBDS,PDS,REFNCE,ISCAL2,KWIDE,IGDS)
C INPUT ARGUMENT LIST:
C IWORK - INTEGER SOURCE ARRAY
C NPTS - NUMBER OF POINTS IN IWORK
C IBDSFL - FLAGS
C
C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS)
C IPFLD - CONTAINS BDS FROM BYTE 12 ON
C BDS11 - CONTAINS FIRST 11 BYTES FOR BDS
C LEN - NUMBER OF BYTES FROM 12 ON
C LENBDS - TOTAL LENGTH OF BDS
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C LANGUAGE: IBM VS FORTRAN 77, CRAY CFT77 FORTRAN
C MACHINE: HDS, CRAY C916/256, Y-MP8/64, Y-MP EL92/256
C
C$$$
CHARACTER*1 BDS11(*),PDS(*),IPFLD(*)
C
REAL REFNCE
C
INTEGER ISCAL2,KWIDE
INTEGER LENBDS
INTEGER IGDS(*)
INTEGER LEN,KBDS(22)
INTEGER IWORK(*)
C OCTET NUMBER IN SECTION, FIRST ORDER PACKING
C INTEGER KBDS(12)
C FLAGS
INTEGER IBDSFL(*)
C EXTENDED FLAGS
C INTEGER KBDS(14)
C OCTET NUMBER FOR SECOND ORDER PACKING
C INTEGER KBDS(15)
C NUMBER OF FIRST ORDER VALUES
C INTEGER KBDS(17)
C NUMBER OF SECOND ORDER PACKED VALUES
C INTEGER KBDS(19)
C WIDTH OF SECOND ORDER PACKING
character(len=1) ISOWID(400000)
C SECONDARY BIT MAP
character(len=1) ISOBMP(65600)
C FIRST ORDER PACKED VALUES
character(len=1) IFOVAL(400000)
C SECOND ORDER PACKED VALUES
character(len=1) ISOVAL(800000)
C
C INTEGER KBDS(11)
C ----------------------------------
C INITIALIZE ARRAYS
C
DO I = 1, 400000
IFOVAL(I) = char(0)
ISOWID(I) = char(0)
ENDDO
C
DO 101 I = 1, 65600
ISOBMP(I) = char(0)
101 CONTINUE
DO 102 I = 1, 800000
ISOVAL(I) = char(0)
102 CONTINUE
C INITIALIZE POINTERS
C SECONDARY BIT WIDTH POINTER
IWDPTR = 0
C SECONDARY BIT MAP POINTER
IBMP2P = 0
C FIRST ORDER VALUE POINTER
IFOPTR = 0
C BYTE POINTER TO START OF 1ST ORDER VALUES
KBDS(12) = 0
C BYTE POINTER TO START OF 2ND ORDER VALUES
KBDS(15) = 0
C TO CONTAIN NUMBER OF FIRST ORDER VALUES
KBDS(17) = 0
C TO CONTAIN NUMBER OF SECOND ORDER VALUES
KBDS(19) = 0
C SECOND ORDER PACKED VALUE POINTER
ISOPTR = 0
C =======================================================
C BUILD SECOND ORDER BIT MAP IN EITHER
C ROW BY ROW OR COL BY COL FORMAT
IF (IAND(IGDS(13),32).NE.0) THEN
C COLUMN BY COLUMN
KOUT = IGDS(4)
KIN = IGDS(5)
C PRINT *,'COLUMN BY COLUMN',KOUT,KIN
ELSE
C ROW BY ROW
KOUT = IGDS(5)
KIN = IGDS(4)
C PRINT *,'ROW BY ROW',KOUT,KIN
END IF
KBDS(17) = KOUT
KBDS(19) = NPTS
C
C DO 4100 J = 1, NPTS, 53
C WRITE (6,4101) (IWORK(K),K=J,J+52)
4101 FORMAT (1X,25I4)
C PRINT *,' '
C4100 CONTINUE
C
C INITIALIZE BIT MAP POINTER
IBMP2P = 0
C CONSTRUCT WORKING BIT MAP
DO 2000 I = 1, KOUT
DO 1000 J = 1, KIN
IF (J.EQ.1) THEN
CALL SBYTEC (ISOBMP,1,IBMP2P,1)
ELSE
CALL SBYTEC (ISOBMP,0,IBMP2P,1)
END IF
IBMP2P = IBMP2P + 1
1000 CONTINUE
2000 CONTINUE
LEN = IBMP2P / 32 + 1
C CALL BINARY(ISOBMP,LEN)
C
C PROCESS OUTER LOOP OF ROW BY ROW OR COL BY COL
C
KPTR = 1
KBDS(11) = KWIDE
DO 6000 I = 1, KOUT
C IN CURRENT ROW OR COL
C FIND FIRST ORDER VALUE
JPTR = KPTR
LOWEST = IWORK(JPTR)
DO 4000 J = 1, KIN
IF (IWORK(JPTR).LT.LOWEST) THEN
LOWEST = IWORK(JPTR)
END IF
JPTR = JPTR + 1
4000 CONTINUE
C SAVE FIRST ORDER VALUE
CALL SBYTEC (IFOVAL,LOWEST,IFOPTR,KWIDE)
IFOPTR = IFOPTR + KWIDE
C PRINT *,'FOVAL',I,LOWEST,KWIDE
C SUBTRACT FIRST ORDER VALUE FROM OTHER VALS
C GETTING SECOND ORDER VALUES
JPTR = KPTR
IBIG = IWORK(JPTR) - LOWEST
DO 4200 J = 1, KIN
IWORK(JPTR) = IWORK(JPTR) - LOWEST
IF (IWORK(JPTR).GT.IBIG) THEN
IBIG = IWORK(JPTR)
END IF
JPTR = JPTR + 1
4200 CONTINUE
C HOW MANY BITS TO CONTAIN LARGEST SECOND
C ORDER VALUE IN SEGMENT
CALL FI7505 (IBIG,NWIDE)
C SAVE BIT WIDTH
CALL SBYTEC (ISOWID,NWIDE,IWDPTR,8)
IWDPTR = IWDPTR + 8
C PRINT *,I,'SOVAL',IBIG,' IN',NWIDE,' BITS'
C WRITE (6,4101) (IWORK(K),K=KPTR,KPTR+52)
C SAVE SECOND ORDER VALUES OF THIS SEGMENT
DO 5000 J = 0, KIN-1
CALL SBYTEC (ISOVAL,IWORK(KPTR+J),ISOPTR,NWIDE)
ISOPTR = ISOPTR + NWIDE
5000 CONTINUE
KPTR = KPTR + KIN
6000 CONTINUE
C =======================================================
C CONCANTENATE ALL FIELDS FOR BDS
C
C REMAINDER GOES INTO IPFLD
IPTR = 0
C BYTES 12-13
C VALUE FOR N1
C LEAVE SPACE FOR THIS
IPTR = IPTR + 16
C BYTE 14
C EXTENDED FLAGS
CALL SBYTEC (IPFLD,IBDSFL(5),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(6),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(7),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(8),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(9),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(10),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(11),IPTR,1)
IPTR = IPTR + 1
CALL SBYTEC (IPFLD,IBDSFL(12),IPTR,1)
IPTR = IPTR + 1
C BYTES 15-16
C SKIP OVER VALUE FOR N2
IPTR = IPTR + 16
C BYTES 17-18
C P1
CALL SBYTEC (IPFLD,KBDS(17),IPTR,16)
IPTR = IPTR + 16
C BYTES 19-20
C P2
CALL SBYTEC (IPFLD,KBDS(19),IPTR,16)
IPTR = IPTR + 16
C BYTE 21 - RESERVED LOCATION
CALL SBYTEC (IPFLD,0,IPTR,8)
IPTR = IPTR + 8
C BYTES 22 - ?
C WIDTHS OF SECOND ORDER PACKING
IX = (IWDPTR + 32) / 32
C CALL SBYTESC (IPFLD,ISOWID,IPTR,32,0,IX)
ijk=IWDPTR/8
jst=(iptr/8)+1
ipfld(jst:jst+ijk)=ISOWID(1:ijk)
IPTR = IPTR + IWDPTR
C PRINT *,'ISOWID',IWDPTR,IX
C CALL BINARY (ISOWID,IX)
C
C NO SECONDARY BIT MAP
C DETERMINE LOCATION FOR START
C OF FIRST ORDER PACKED VALUES
KBDS(12) = IPTR / 8 + 12
C STORE LOCATION
CALL SBYTEC (IPFLD,KBDS(12),0,16)
C MOVE IN FIRST ORDER PACKED VALUES
IPASS = (IFOPTR + 32) / 32
c CALL SBYTESC (IPFLD,IFOVAL,IPTR,32,0,IPASS)
ijk=(IFOPTR/8)+1
jst=(iptr/8)+1
ipfld(jst:jst+ijk)=ifoval(1:ijk)
IPTR = IPTR + IFOPTR
C PRINT *,'IFOVAL',IFOPTR,IPASS,KWIDE
C CALL BINARY (IFOVAL,IPASS)
IF (MOD(IPTR,8).NE.0) THEN
IPTR = IPTR + 8 - MOD(IPTR,8)
END IF
C PRINT *,'IFOPTR =',IFOPTR,' ISOPTR =',ISOPTR
C DETERMINE LOCATION FOR START
C OF SECOND ORDER VALUES
KBDS(15) = IPTR / 8 + 12
C SAVE LOCATION OF SECOND ORDER VALUES
CALL SBYTEC (IPFLD,KBDS(15),24,16)
C MOVE IN SECOND ORDER PACKED VALUES
IX = (ISOPTR + 32) / 32
C CALL SBYTESC (IPFLD,ISOVAL,IPTR,32,0,IX)
ijk=(ISOPTR/8)+1
jst=(iptr/8)+1
ipfld(jst:jst+ijk)=isoval(1:ijk)
IPTR = IPTR + ISOPTR
C PRINT *,'ISOVAL',ISOPTR,IX
C CALL BINARY (ISOVAL,IX)
NLEFT = MOD(IPTR+88,16)
IF (NLEFT.NE.0) THEN
NLEFT = 16 - NLEFT
IPTR = IPTR + NLEFT
END IF
C COMPUTE LENGTH OF DATA PORTION
LEN = IPTR / 8
C COMPUTE LENGTH OF BDS
LENBDS = LEN + 11
C -----------------------------------
C BYTES 1-3
C THIS FUNCTION COMPLETED BELOW
C WHEN LENGTH OF BDS IS KNOWN
CALL SBYTEC (BDS11,LENBDS,0,24)
C BYTE 4
CALL SBYTEC (BDS11,IBDSFL(1),24,1)
CALL SBYTEC (BDS11,IBDSFL(2),25,1)
CALL SBYTEC (BDS11,IBDSFL(3),26,1)
CALL SBYTEC (BDS11,IBDSFL(4),27,1)
C ENTER NUMBER OF FILL BITS
CALL SBYTEC (BDS11,NLEFT,28,4)
C BYTE 5-6
IF (ISCAL2.LT.0) THEN
CALL SBYTEC (BDS11,1,32,1)
ISCAL2 = - ISCAL2
ELSE
CALL SBYTEC (BDS11,0,32,1)
END IF
CALL SBYTEC (BDS11,ISCAL2,33,15)
C
C$ FILL OCTET 7-10 WITH THE REFERENCE VALUE
C CONVERT THE FLOATING POINT OF YOUR MACHINE TO IBM370 32 BIT
C FLOATING POINT NUMBER
C REFERENCE VALUE
C FIRST TEST TO SEE IF
C ON 32 OR 64 BIT COMPUTER
C CALL W3FI01(LW)
IF (bit_size(LW).EQ.32) THEN
CALL W3FI76 (REFNCE,IEXP,IMANT,32)
ELSE
CALL W3FI76 (REFNCE,IEXP,IMANT,64)
END IF
CALL SBYTEC (BDS11,IEXP,48,8)
CALL SBYTEC (BDS11,IMANT,56,24)
C
C BYTE 11
C
CALL SBYTEC (BDS11,KBDS(11),80,8)
C
KLEN = LENBDS / 4 + 1
C PRINT *,'BDS11 LISTING',4,LENBDS
C CALL BINARY (BDS11,4)
C PRINT *,'IPFLD LISTING'
C CALL BINARY (IPFLD,KLEN)
RETURN
END
SUBROUTINE FI7505 (N,NBITS)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . . .
C SUBPROGRAM: FI7505 DETERMINE NUMBER OF BITS TO CONTAIN VALUE
C PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 93-06-23
C
C ABSTRACT: CALCULATE NUMBER OF BITS TO CONTAIN VALUE N, WITH A
C MAXIMUM OF 32 BITS.
C
C PROGRAM HISTORY LOG:
C 93-06-23 CAVANAUGH
C 95-10-31 IREDELL REMOVED SAVES AND PRINTS
C
C USAGE: CALL FI7505 (N,NBITS)
C INPUT ARGUMENT LIST:
C N - INTEGER VALUE
C
C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS)
C NBITS - NUMBER OF BITS TO CONTAIN N
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C LANGUAGE: IBM VS FORTRAN 77, CRAY CFT77 FORTRAN
C MACHINE: HDS, CRAY C916/256, Y-MP8/64, Y-MP EL92/256
C
C$$$
INTEGER N,NBITS
INTEGER IBITS(31)
C
DATA IBITS/1,3,7,15,31,63,127,255,511,1023,2047,
* 4095,8191,16383,32767,65535,131071,262143,
* 524287,1048575,2097151,4194303,8388607,
* 16777215,33554431,67108863,134217727,268435455,
* 536870911,1073741823,2147483647/
C ----------------------------------------------------------------
C
DO 1000 NBITS = 1, 31
IF (N.LE.IBITS(NBITS)) THEN
RETURN
END IF
1000 CONTINUE
RETURN
END
SUBROUTINE FI7513 (IWORK,ISTART,NPTS,MAX,MIN,INRNGE)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . . .
C SUBPROGRAM: FI7513 SELECT BLOCK OF DATA FOR PACKING
C PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 94-01-21
C
C ABSTRACT: SELECT A BLOCK OF DATA FOR PACKING
C
C PROGRAM HISTORY LOG:
C 94-01-21 CAVANAUGH
C 95-10-31 IREDELL REMOVED SAVES AND PRINTS
C
C USAGE: CALL FI7513 (IWORK,ISTART,NPTS,MAX,MIN,INRNGE)
C INPUT ARGUMENT LIST:
C * - RETURN ADDRESS IF ENCOUNTER SET OF SAME VALUES
C IWORK -
C ISTART -
C NPTS -
C
C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS)
C MAX -
C MIN -
C INRNGE -
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C LANGUAGE: IBM VS FORTRAN 77, CRAY CFT77 FORTRAN
C MACHINE: HDS, CRAY C916/256, Y-MP8/64, Y-MP EL92/256
C
C$$$
INTEGER IWORK(*),NPTS,ISTART,INRNGE,INRNGA,INRNGB
INTEGER MAX,MIN,MXVAL,MAXB,MINB,MXVALB
INTEGER IBITS(31)
C
DATA IBITS/1,3,7,15,31,63,127,255,511,1023,2047,
* 4095,8191,16383,32767,65535,131071,262143,
* 524287,1048575,2097151,4194303,8388607,
* 16777215,33554431,67108863,134217727,268435455,
* 536870911,1073741823,2147483647/
C ----------------------------------------------------------------
C IDENTIFY NEXT BLOCK OF DATA FOR PACKING AND
C RETURN TO CALLER
C ********************************************************************
ISTRTA = ISTART
C
C GET BLOCK A
CALL FI7516 (IWORK,NPTS,INRNGA,ISTRTA,
* MAX,MIN,MXVAL,LWIDE)
C ********************************************************************
C
ISTRTB = ISTRTA + INRNGA
2000 CONTINUE
C IF HAVE PROCESSED ALL DATA, RETURN
IF (ISTRTB.GT.NPTS) THEN
C NO MORE DATA TO LOOK AT
INRNGE = INRNGA
RETURN
END IF
C GET BLOCK B
CALL FI7502 (IWORK,ISTRTB,NPTS,ISAME)
IF (ISAME.GE.15) THEN
C PRINT *,'BLOCK B HAS ALL IDENTICAL VALUES'
C PRINT *,'BLOCK A HAS INRNGE =',INRNGA
C BLOCK B CONTAINS ALL IDENTICAL VALUES
INRNGE = INRNGA
C EXIT WITH BLOCK A
RETURN
END IF
C GET BLOCK B
C
ISTRTB = ISTRTA + INRNGA
CALL FI7516 (IWORK,NPTS,INRNGB,ISTRTB,
* MAXB,MINB,MXVALB,LWIDEB)
C PRINT *,'BLOCK A',INRNGA,' BLOCK B',INRNGB
C ********************************************************************
C PERFORM TREND ANALYSIS TO DETERMINE
C IF DATA COLLECTION CAN BE IMPROVED
C
KTRND = LWIDE - LWIDEB
C PRINT *,'TREND',LWIDE,LWIDEB
IF (KTRND.LE.0) THEN
C PRINT *,'BLOCK A - SMALLER, SHOULD EXTEND INTO BLOCK B'
MXVAL = IBITS(LWIDE)
C
C IF BLOCK A REQUIRES THE SAME OR FEWER BITS
C LOOK AHEAD
C AND GATHER THOSE DATA POINTS THAT CAN
C BE RETAINED IN BLOCK A
C BECAUSE THIS BLOCK OF DATA
C USES FEWER BITS
C
CALL FI7518 (IRET,IWORK,NPTS,ISTRTA,INRNGA,INRNGB,
* MAX,MIN,LWIDE,MXVAL)
IF(IRET.EQ.1) GO TO 8000
C PRINT *,'18 INRNGA IS NOW ',INRNGA
IF (INRNGB.LT.20) THEN
RETURN
ELSE
GO TO 2000
END IF
ELSE
C PRINT *,'BLOCK A - LARGER, B SHOULD EXTEND BACK INTO A'
MXVALB = IBITS(LWIDEB)
C
C IF BLOCK B REQUIRES FEWER BITS
C LOOK BACK
C SHORTEN BLOCK A BECAUSE NEXT BLOCK OF DATA
C USES FEWER BITS
C
CALL FI7517 (IRET,IWORK,NPTS,ISTRTB,INRNGA,
* MAXB,MINB,LWIDEB,MXVALB)
IF(IRET.EQ.1) GO TO 8000
C PRINT *,'17 INRNGA IS NOW ',INRNGA
END IF
C
C PACK UP BLOCK A
C UPDATA POINTERS
8000 CONTINUE
INRNGE = INRNGA
C GET NEXT BLOCK A
9000 CONTINUE
RETURN
END
SUBROUTINE FI7516 (IWORK,NPTS,INRNG,ISTART,MAX,MIN,MXVAL,LWIDTH)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . . .
C SUBPROGRAM: FI7516 SCAN NUMBER OF POINTS
C PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 94-01-21
C
C ABSTRACT: SCAN FORWARD FROM CURRENT POSITION. COLLECT POINTS AND
C DETERMINE MAXIMUM AND MINIMUM VALUES AND THE NUMBER
C OF POINTS THAT ARE INCLUDED. FORWARD SEARCH IS TERMINATED
C BY ENCOUNTERING A SET OF IDENTICAL VALUES, BY REACHING
C THE NUMBER OF POINTS SELECTED OR BY REACHING THE END
C OF DATA.
C
C PROGRAM HISTORY LOG:
C 94-01-21 CAVANAUGH
C 95-10-31 IREDELL REMOVED SAVES AND PRINTS
C
C USAGE: CALL FI7516 (IWORK,NPTS,INRNG,ISTART,MAX,MIN,MXVAL,LWIDTH)
C INPUT ARGUMENT LIST:
C * - RETURN ADDRESS IF ENCOUNTER SET OF SAME VALUES
C IWORK - DATA ARRAY
C NPTS - NUMBER OF POINTS IN DATA ARRAY
C ISTART - STARTING LOCATION IN DATA
C
C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS)
C INRNG - NUMBER OF POINTS SELECTED
C MAX - MAXIMUM VALUE OF POINTS
C MIN - MINIMUM VALUE OF POINTS
C MXVAL - MAXIMUM VALUE THAT CAN BE CONTAINED IN LWIDTH BITS
C LWIDTH - NUMBER OF BITS TO CONTAIN MAX DIFF
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C LANGUAGE: IBM VS FORTRAN 77, CRAY CFT77 FORTRAN
C MACHINE: HDS, CRAY C916/256, Y-MP8/64, Y-MP EL92/256
C
C$$$
INTEGER IWORK(*),NPTS,ISTART,INRNG,MAX,MIN,LWIDTH,MXVAL
INTEGER IBITS(31)
C
DATA IBITS/1,3,7,15,31,63,127,255,511,1023,2047,
* 4095,8191,16383,32767,65535,131071,262143,
* 524287,1048575,2097151,4194303,8388607,
* 16777215,33554431,67108863,134217727,268435455,
* 536870911,1073741823,2147483647/
C ----------------------------------------------------------------
C
INRNG = 1
JQ = ISTART + 19
MAX = IWORK(ISTART)
MIN = IWORK(ISTART)
DO 1000 I = ISTART+1, JQ
CALL FI7502 (IWORK,I,NPTS,ISAME)
IF (ISAME.GE.15) THEN
GO TO 5000
END IF
INRNG = INRNG + 1
IF (IWORK(I).GT.MAX) THEN
MAX = IWORK(I)
ELSE IF (IWORK(I).LT.MIN) THEN
MIN = IWORK(I)
END IF
1000 CONTINUE
5000 CONTINUE
KRNG = MAX - MIN
C
DO 9000 LWIDTH = 1, 31
IF (KRNG.LE.IBITS(LWIDTH)) THEN
C PRINT *,'RETURNED',INRNG,' VALUES'
RETURN
END IF
9000 CONTINUE
RETURN
END
SUBROUTINE FI7517 (IRET,IWORK,NPTS,ISTRTB,INRNGA,
* MAXB,MINB,MXVALB,LWIDEB)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . . .
C SUBPROGRAM: FI7517 SCAN BACKWARD
C PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 94-01-21
C
C ABSTRACT: SCAN BACKWARDS UNTIL A VALUE EXCEEDS RANGE OF GROUP B
C THIS MAY SHORTEN GROUP A
C
C PROGRAM HISTORY LOG:
C 94-01-21 CAVANAUGH
C 95-10-31 IREDELL REMOVED SAVES AND PRINTS
C 98-06-17 IREDELL REMOVED ALTERNATE RETURN
C
C USAGE: CALL FI7517 (IRET,IWORK,NPTS,ISTRTB,INRNGA,
C * MAXB,MINB,MXVALB,LWIDEB)
C INPUT ARGUMENT LIST:
C IWORK -
C ISTRTB -
C NPTS -
C INRNGA -
C
C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS)
C IRET -
C JLAST -
C MAXB -
C MINB -
C LWIDTH - NUMBER OF BITS TO CONTAIN MAX DIFF
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C LANGUAGE: IBM VS FORTRAN 77, CRAY CFT77 FORTRAN
C MACHINE: HDS, CRAY C916/256, Y-MP8/64, Y-MP EL92/256
C
C$$$
INTEGER IWORK(*),NPTS,ISTRTB,INRNGA
INTEGER MAXB,MINB,LWIDEB,MXVALB
INTEGER IBITS(31)
C
DATA IBITS/1,3,7,15,31,63,127,255,511,1023,2047,
* 4095,8191,16383,32767,65535,131071,262143,
* 524287,1048575,2097151,4194303,8388607,
* 16777215,33554431,67108863,134217727,268435455,
* 536870911,1073741823,2147483647/
C ----------------------------------------------------------------
IRET=0
C PRINT *,' FI7517'
NPOS = ISTRTB - 1
ITST = 0
KSET = INRNGA
C
1000 CONTINUE
C PRINT *,'TRY NPOS',NPOS,IWORK(NPOS),MAXB,MINB
ITST = ITST + 1
IF (ITST.LE.KSET) THEN
IF (IWORK(NPOS).GT.MAXB) THEN
IF ((IWORK(NPOS)-MINB).GT.MXVALB) THEN
C PRINT *,'WENT OUT OF RANGE AT',NPOS
IRET=1
RETURN
ELSE
MAXB = IWORK(NPOS)
END IF
ELSE IF (IWORK(NPOS).LT.MINB) THEN
IF ((MAXB-IWORK(NPOS)).GT.MXVALB) THEN
C PRINT *,'WENT OUT OF RANGE AT',NPOS
IRET=1
RETURN
ELSE
MINB = IWORK(NPOS)
END IF
END IF
INRNGA = INRNGA - 1
NPOS = NPOS - 1
GO TO 1000
END IF
C ----------------------------------------------------------------
C
9000 CONTINUE
RETURN
END
SUBROUTINE FI7518 (IRET,IWORK,NPTS,ISTRTA,INRNGA,INRNGB,
* MAXA,MINA,LWIDEA,MXVALA)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . . .
C SUBPROGRAM: FI7518 SCAN FORWARD
C PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 94-01-21
C
C ABSTRACT: SCAN FORWARD FROM START OF BLOCK B TOWARDS END OF BLOCK B
C IF NEXT POINT UNDER TEST FORCES A LARGER MAXVALA THEN
C TERMINATE INDICATING LAST POINT TESTED FOR INCLUSION
C INTO BLOCK A.
C
C PROGRAM HISTORY LOG:
C 94-01-21 CAVANAUGH
C 95-10-31 IREDELL REMOVED SAVES AND PRINTS
C 98-06-17 IREDELL REMOVED ALTERNATE RETURN
C
C USAGE: CALL FI7518 (IRET,IWORK,NPTS,ISTRTA,INRNGA,INRNGB,
C * MAXA,MINA,LWIDEA,MXVALA)
C INPUT ARGUMENT LIST:
C IFLD -
C JSTART -
C NPTS -
C
C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS)
C IRET -
C JLAST -
C MAX -
C MIN -
C LWIDTH - NUMBER OF BITS TO CONTAIN MAX DIFF
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C LANGUAGE: IBM VS FORTRAN 77, CRAY CFT77 FORTRAN
C MACHINE: HDS, CRAY C916/256, Y-MP8/64, Y-MP EL92/256
C
C$$$
INTEGER IWORK(*),NPTS,ISTRTA,INRNGA
INTEGER MAXA,MINA,LWIDEA,MXVALA
INTEGER IBITS(31)
C
DATA IBITS/1,3,7,15,31,63,127,255,511,1023,2047,
* 4095,8191,16383,32767,65535,131071,262143,
* 524287,1048575,2097151,4194303,8388607,
* 16777215,33554431,67108863,134217727,268435455,
* 536870911,1073741823,2147483647/
C ----------------------------------------------------------------
IRET=0
C PRINT *,' FI7518'
NPOS = ISTRTA + INRNGA
ITST = 0
C
1000 CONTINUE
ITST = ITST + 1
IF (ITST.LE.INRNGB) THEN
C PRINT *,'TRY NPOS',NPOS,IWORK(NPOS),MAXA,MINA
IF (IWORK(NPOS).GT.MAXA) THEN
IF ((IWORK(NPOS)-MINA).GT.MXVALA) THEN
C PRINT *,'FI7518A -',ITST,' RANGE EXCEEDS MAX'
IRET=1
RETURN
ELSE
MAXA = IWORK(NPOS)
END IF
ELSE IF (IWORK(NPOS).LT.MINA) THEN
IF ((MAXA-IWORK(NPOS)).GT.MXVALA) THEN
C PRINT *,'FI7518B -',ITST,' RANGE EXCEEDS MAX'
IRET=1
RETURN
ELSE
MINA = IWORK(NPOS)
END IF
END IF
INRNGA = INRNGA + 1
C PRINT *,' ',ITST,INRNGA
NPOS = NPOS +1
GO TO 1000
END IF
C ----------------------------------------------------------------
9000 CONTINUE
RETURN
END