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