! !################################################################## !################################################################## !###### ###### !###### NMCDECODE ###### !###### ###### !################################################################## !################################################################## ! ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Decode the NMC GRIB files. ! !----------------------------------------------------------------------- ! ! AUTHOR: NWS/NMC ! ! MODIFICATION HISTORY: ! 01/16/1996 (Yuhe Liu) ! Changed subroutine name of SBYTE(S) and GBYTE(S) to GRBSBYTE(S) ! and GRBGBYTE(S) to avoid the conflict with the subroutines of ! same names in NCAR graphics package. ! !----------------------------------------------------------------------- ! ! Original comments: ! ! FILE UPDATED ON 11-22-94 BY M. FARLEY. ! CONTAINS SUBROUTINES NEEDED TO DECODE GRIB MESSAGES: ! W3FI63 AND W3FI83. ! SUBROUTINE w3fi63(msga,kpds,kgds,kbms,DATA,kptr,kret) 1,24 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: W3FI63 UNPK GRIB FIELD TO GRIB GRID ! PRGMMR: FARLEY ORG: NMC421 DATE:94-11-22 ! ! ABSTRACT: UNPACK A GRIB (EDITION 1) FIELD TO THE EXACT GRID ! SPECIFIED IN THE GRIB MESSAGE, ISOLATE THE BIT MAP, AND MAKE ! THE VALUES OF THE PRODUCT DESCRIPTON SECTION (PDS) AND THE ! GRID DESCRIPTION SECTION (GDS) AVAILABLE IN RETURN ARRAYS. ! ! WHEN DECODING IS COMPLETED, DATA AT EACH GRID POINT HAS BEEN ! RETURNED IN THE UNITS SPECIFIED IN THE GRIB MANUAL. ! ! PROGRAM HISTORY LOG: ! 91-09-13 CAVANAUGH ! 91-11-12 CAVANAUGH MODIFIED SIZE OF ECMWF GRIDS 5-8 ! 91-12-22 CAVANAUGH CORRECTED PROCESSING OF MERCATOR PROJECTIONS ! IN GRID DEFINITION SECTION (GDS) IN ! ROUTINE FI633 ! 92-08-05 CAVANAUGH CORRECTED MAXIMUM GRID SIZE TO ALLOW FOR ! ONE DEGREE BY ONE DEGREE GLOBAL GRIDS ! 92-08-27 CAVANAUGH CORRECTED TYPO ERROR, ADDED CODE TO COMPARE ! TOTAL BYTE SIZE FROM SECTION 0 WITH SUM OF ! SECTION SIZES. ! 92-10-21 CAVANAUGH CORRECTIONS WERE MADE (IN FI634) TO REDUCE ! PROCESSING TIME FOR INTERNATIONAL GRIDS. ! REMOVED A TYPOGRAPHICAL ERROR IN FI635. ! 93-01-07 CAVANAUGH CORRECTIONS WERE MADE (IN FI635) TO ! FACILITATE USE OF THESE ROUTINES ON A PC. ! A TYPOGRAPHICAL ERROR WAS ALSO CORRECTED ! 93-01-13 CAVANAUGH CORRECTIONS WERE MADE (IN FI632) TO ! PROPERLY HANDLE CONDITION WHEN ! TIME RANGE INDICATOR = 10. ! ADDED U.S.GRID 87. ! 93-02-04 CAVANAUGH ADDED U.S.GRIDS 85 AND 86 ! 93-02-26 CAVANAUGH ADDED GRIDS 2, 3, 37 THRU 44,AND ! GRIDS 55, 56, 90, 91, 92, AND 93 TO ! LIST OF U.S. GRIDS. ! 93-04-07 CAVANAUGH ADDED GRIDS 67 THRU 77 TO ! LIST OF U.S. GRIDS. ! 93-04-20 CAVANAUGH INCREASED MAX SIZE TO ACCOMODATE ! GAUSSIAN GRIDS. ! 93-05-26 CAVANAUGH CORRECTED GRID RANGE SELECTION IN FI634 ! FOR RANGES 67-71 & 75-77 ! 93-06-08 CAVANAUGH CORRECTED FI635 TO ACCEPT GRIB MESSAGES ! WITH SECOND ORDER PACKING. ADDED ROUTINE FI636 ! TO PROCESS MESSAGES WITH SECOND ORDER PACKING. ! 93-09-22 CAVANAUGH MODIFIED TO EXTRACT SUB-CENTER NUMBER FROM ! PDS BYTE 26 ! 93-10-13 CAVANAUGH MODIFIED FI634 TO CORRECT GRID SIZES FOR ! GRIDS 204 AND 208 ! 93-10-14 CAVANAUGH INCREASED SIZE OF KGDS TO INCLUDE ENTRIES FOR ! NUMBER OF POINTS IN GRID AND NUMBER OF WORDS ! IN EACH ROW ! 93-12-08 CAVANAUGH CORRECTED TEST FOR EDITION NUMBER INSTEAD ! OF VERSION NUMBER ! 93-12-15 CAVANAUGH MODIFIED SECOND ORDER POINTERS TO FIRST ORDER ! VALUES AND SECOND ORDER VALUES CORRECTLY ! IN ROUTINE FI636 ! 94-03-02 CAVANAUGH ADDED CALL TO W3FI83 WITHIN DECODER. USER ! NO LONGER NEEDS TO MAKE CALL TO THIS ROUTINE ! 94-04-22 CAVANAUGH MODIFIED FI635, FI636 TO PROCESS ROW BY ROW ! SECOND ORDER PACKING, ADDED SCALING CORRECTION ! TO FI635, AND CORRECTED TYPOGRAPHICAL ERRORS ! IN COMMENT FIELDS IN FI634 ! 94-05-17 CAVANAUGH CORRECTED ERROR IN FI633 TO EXTRACT RESOLUTION ! FOR LAMBERT-CONFORMAL GRIDS. ADDED CLARIFYING ! INFORMATION TO DOCBLOCK ENTRIES ! 94-05-25 CAVANAUGH ADDED CODE TO PROCESS COLUMN BY COLUMN AS WELL ! AS ROW BY ROW ORDERING OF SECOND ORDER DATA ! 94-06-27 CAVANAUGH ADDED PROCESSING FOR GRIDS 45, 94 AND 95. ! INCLUDES CONSTRUCTION OF SECOND ORDER BIT MAPS ! FOR THINNED GRIDS IN FI636. ! 94-07-08 CAVANAUGH COMMENTED OUT PRINT OUTS USED FOR DEBUGGING ! 94-09-08 CAVANAUGH ADDED GRIDS 220, 221, 223 FOR FNOC ! 94-11-10 FARLEY INCREASED MXSIZE FROM 72960 TO 260000 ! FOR .5 DEGREE SST ANALYSIS FIELDS ! ! ! USAGE: CALL W3FI63(MSGA,KPDS,KGDS,KBMS,DATA,KPTR,KRET) ! INPUT ARGUMENT LIST: ! MSGA - GRIB FIELD - "GRIB" THRU "7777" CHAR*1 ! (MESSAGE CAN BE PRECEDED BY JUNK CHARS) ! ! OUTPUT ARGUMENT LIST: ! DATA - ARRAY CONTAINING DATA ELEMENTS ! KPDS - ARRAY CONTAINING PDS ELEMENTS. (EDITION 1) ! (1) - ID OF CENTER ! (2) - GENERATING PROCESS ID NUMBER ! (3) - GRID DEFINITION ! (4) - GDS/BMS FLAG (RIGHT ADJ COPY OF OCTET 8) ! (5) - INDICATOR OF PARAMETER ! (6) - TYPE OF LEVEL ! (7) - HEIGHT/PRESSURE , ETC OF LEVEL ! (8) - YEAR INCLUDING (CENTURY-1) ! (9) - MONTH OF YEAR ! (10) - DAY OF MONTH ! (11) - HOUR OF DAY ! (12) - MINUTE OF HOUR ! (13) - INDICATOR OF FORECAST TIME UNIT ! (14) - TIME RANGE 1 ! (15) - TIME RANGE 2 ! (16) - TIME RANGE FLAG ! (17) - NUMBER INCLUDED IN AVERAGE ! (18) - VERSION NR OF GRIB SPECIFICATION ! (19) - VERSION NR OF PARAMETER TABLE ! (20) - NR MISSING FROM AVERAGE/ACCUMULATION ! (21) - CENTURY OF REFERENCE TIME OF DATA ! (22) - UNITS DECIMAL SCALE FACTOR ! (23) - SUBCENTER NUMBER ! (26-35) - RESERVED ! (24-25) - RESERVED FOR FUTURE USE ! (36-N) - CONSECUTIVE BYTES EXTRACTED FROM PROGRAM ! DEFINITION SECTION (PDS) OF GRIB MESSAGE ! KGDS - ARRAY CONTAINING GDS ELEMENTS. ! (1) - DATA REPRESENTATION TYPE ! (19) - NUMBER OF VERTICAL COORDINATE PARAMETERS ! (20) - OCTET NUMBER OF THE LIST OF VERTICAL COORDINATE ! PARAMETERS ! OR ! OCTET NUMBER OF THE LIST OF NUMBERS OF POINTS ! IN EACH ROW ! OR ! 255 IF NEITHER ARE PRESENT ! (21) - FOR GRIDS WITH PL, NUMBER OF POINTS IN GRID ! (22) - NUMBER OF WORDS IN EACH ROW ! LATITUDE/LONGITUDE GRIDS ! (2) - N(I) NR POINTS ON LATITUDE CIRCLE ! (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN ! (4) - LA(1) LATITUDE OF ORIGIN ! (5) - LO(1) LONGITUDE OF ORIGIN ! (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17) ! (7) - LA(2) LATITUDE OF EXTREME POINT ! (8) - LO(2) LONGITUDE OF EXTREME POINT ! (9) - DI LATITUDINAL DIRECTION OF INCREMENT ! (10) - DJ LONGITUDINAL DIRECTION INCREMENT ! (11) - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28) ! GAUSSIAN GRIDS ! (2) - N(I) NR POINTS ON LATITUDE CIRCLE ! (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN ! (4) - LA(1) LATITUDE OF ORIGIN ! (5) - LO(1) LONGITUDE OF ORIGIN ! (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17) ! (7) - LA(2) LATITUDE OF EXTREME POINT ! (8) - LO(2) LONGITUDE OF EXTREME POINT ! (9) - DI LATITUDINAL DIRECTION OF INCREMENT ! (10) - N - NR OF CIRCLES POLE TO EQUATOR ! (11) - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28) ! (12) - NV - NR OF VERT COORD PARAMETERS ! (13) - PV - OCTET NR OF LIST OF VERT COORD PARAMETERS ! OR ! PL - LOCATION OF THE LIST OF NUMBERS OF POINTS IN ! EACH ROW (IF NO VERT COORD PARAMETERS ! ARE PRESENT ! OR ! 255 IF NEITHER ARE PRESENT ! POLAR STEREOGRAPHIC GRIDS ! (2) - N(I) NR POINTS ALONG LAT CIRCLE ! (3) - N(J) NR POINTS ALONG LON CIRCLE ! (4) - LA(1) LATITUDE OF ORIGIN ! (5) - LO(1) LONGITUDE OF ORIGIN ! (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17) ! (7) - LOV GRID ORIENTATION ! (8) - DX - X DIRECTION INCREMENT ! (9) - DY - Y DIRECTION INCREMENT ! (10) - PROJECTION CENTER FLAG ! (11) - SCANNING MODE (RIGHT ADJ COPY OF OCTET 28) ! SPHERICAL HARMONIC COEFFICIENTS ! (2) - J PENTAGONAL RESOLUTION PARAMETER ! (3) - K " " " ! (4) - M " " " ! (5) - REPRESENTATION TYPE ! (6) - COEFFICIENT STORAGE MODE ! MERCATOR GRIDS ! (2) - N(I) NR POINTS ON LATITUDE CIRCLE ! (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN ! (4) - LA(1) LATITUDE OF ORIGIN ! (5) - LO(1) LONGITUDE OF ORIGIN ! (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17) ! (7) - LA(2) LATITUDE OF LAST GRID POINT ! (8) - LO(2) LONGITUDE OF LAST GRID POINT ! (9) - LATIT - LATITUDE OF PROJECTION INTERSECTION ! (10) - RESERVED ! (11) - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28) ! (12) - LONGITUDINAL DIR GRID LENGTH ! (13) - LATITUDINAL DIR GRID LENGTH ! LAMBERT CONFORMAL GRIDS ! (2) - NX NR POINTS ALONG X-AXIS ! (3) - NY NR POINTS ALONG Y-AXIS ! (4) - LA1 LAT OF ORIGIN (LOWER LEFT) ! (5) - LO1 LON OF ORIGIN (LOWER LEFT) ! (6) - RESOLUTION (RIGHT ADJ COPY OF OCTET 17) ! (7) - LOV - ORIENTATION OF GRID ! (8) - DX - X-DIR INCREMENT ! (9) - DY - Y-DIR INCREMENT ! (10) - PROJECTION CENTER FLAG ! (11) - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28) ! (12) - LATIN 1 - FIRST LAT FROM POLE OF SECANT CONE INTER ! (13) - LATIN 2 - SECOND LAT FROM POLE OF SECANT CONE INTER ! KBMS - BITMAP DESCRIBING LOCATION OF OUTPUT ELEMENTS. ! (ALWAYS CONSTRUCTED) ! KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS ! (1) - TOTAL LENGTH OF GRIB MESSAGE ! (2) - LENGTH OF INDICATOR (SECTION 0) ! (3) - LENGTH OF PDS (SECTION 1) ! (4) - LENGTH OF GDS (SECTION 2) ! (5) - LENGTH OF BMS (SECTION 3) ! (6) - LENGTH OF BDS (SECTION 4) ! (7) - VALUE OF CURRENT BYTE ! (8) - BIT POINTER ! (9) - GRIB START BIT NR ! (10) - GRIB/GRID ELEMENT COUNT ! (11) - NR UNUSED BITS AT END OF SECTION 3 ! (12) - BIT MAP FLAG (COPY OF BMS OCTETS 5,6) ! (13) - NR UNUSED BITS AT END OF SECTION 2 ! (14) - BDS FLAGS (RIGHT ADJ COPY OF OCTET 4) ! (15) - NR UNUSED BITS AT END OF SECTION 4 ! KRET - FLAG INDICATING QUALITY OF COMPLETION ! ! REMARKS: WHEN DECODING IS COMPLETED, DATA AT EACH GRID POINT HAS BEEN ! RETURNED IN THE UNITS SPECIFIED IN THE GRIB MANUAL. ! ! VALUES FOR RETURN FLAG (KRET) ! KRET = 0 - NORMAL RETURN, NO ERRORS ! = 1 - 'GRIB' NOT FOUND IN FIRST 100 CHARS ! = 2 - '7777' NOT IN CORRECT LOCATION ! = 3 - UNPACKED FIELD IS LARGER THAN 65160 ! = 4 - GDS/ GRID NOT ONE OF CURRENTLY ACCEPTED VALUES ! = 5 - GRID NOT CURRENTLY AVAIL FOR CENTER INDICATED ! = 8 - TEMP GDS INDICATED, BUT GDS FLAG IS OFF ! = 9 - GDS INDICATES SIZE MISMATCH WITH STD GRID ! =10 - INCORRECT CENTER INDICATOR ! =11 - BINARY DATA SECTION (BDS) NOT COMPLETELY PROCESSED. ! PROGRAM IS NOT SET TO PROCESS FLAG COMBINATIONS ! SHOWN IN OCTETS 4 AND 14. ! =12 - BINARY DATA SECTION (BDS) NOT COMPLETELY PROCESSED. ! PROGRAM IS NOT SET TO PROCESS FLAG COMBINATIONS ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 77 ! MACHINE: NAS9000 ! !$$$ ! 4 AUG 1988 ! W3FI63 ! ! ! GRIB UNPACKING ROUTINE ! ! ! THIS ROUTINE WILL UNPACK A 'GRIB' FIELD TO THE EXACT GRID ! TYPE SPECIFIED IN THE MESSAGE, RETURN A BIT MAP AND MAKE THE ! VALUES OF THE PRODUCT DEFINITION SEC (PDS) AND THE GRID ! DESCRIPTION SEC (GDS) AVAILABLE IN RETURN ARRAYS. ! SEE "GRIB - THE WMO FORMAT FOR THE STORAGE OF WEATHER PRODUCT ! INFORMATION AND THE EXCHANGE OF WEATHER PRODUCT MESSAGES IN ! GRIDDED BINARY FORM" DATED JULY 1, 1988 BY JOHN D. STACKPOLE ! DOC, NOAA, NWS, NATIONAL METEOROLOGICAL CENTER. ! ! THE CALL TO THE GRIB UNPACKING ROUTINE IS AS FOLLOWS: ! ! CALL W3FI63(MSGA,KPDS,KGDS,LBMS,DATA,KPTR,KRET) ! ! INPUT: ! ! MSGA = CONTAINS THE GRIB MESSAGE TO BE UNPACKED. CHARACTERS ! "GRIB" MAY BEGIN ANYWHERE WITHIN FIRST 100 BYTES. ! ! OUTPUT: ! ! KPDS(100) INTEGER*4 ! ARRAY TO CONTAIN THE ELEMENTS OF THE PRODUCT ! DEFINITION SEC . ! (VERSION 1) ! KPDS(1) - ID OF CENTER ! KPDS(2) - MODEL IDENTIFICATION (SEE "GRIB" TABLE 1) ! KPDS(3) - GRID IDENTIFICATION (SEE "GRIB" TABLE 2) ! KPDS(4) - GDS/BMS FLAG ! BIT DEFINITION ! 25 0 - GDS OMITTED ! 1 - GDS INCLUDED ! 26 0 - BMS OMITTED ! 1 - BMS INCLUDED ! NOTE:- LEFTMOST BIT = 1, ! RIGHTMOST BIT = 32 ! KPDS(5) - INDICATOR OF PARAMETER (SEE "GRIB" TABLE 5) ! KPDS(6) - TYPE OF LEVEL (SEE "GRIB" TABLES 6 & 7) ! KPDS(7) - HEIGHT,PRESSURE,ETC OF LEVEL ! KPDS(8) - YEAR INCLUDING CENTURY ! KPDS(9) - MONTH OF YEAR ! KPDS(10) - DAY OF MONTH ! KPDS(11) - HOUR OF DAY ! KPDS(12) - MINUTE OF HOUR ! KPDS(13) - INDICATOR OF FORECAST TIME UNIT (SEE "GRIB" ! TABLE 8) ! KPDS(14) - TIME 1 (SEE "GRIB" TABLE 8A) ! KPDS(15) - TIME 2 (SEE "GRIB" TABLE 8A) ! KPDS(16) - TIME RANGE INDICATOR (SEE "GRIB" TABLE 8A) ! KPDS(17) - NUMBER INCLUDED IN AVERAGE ! KPDS(18) - EDITION NR OF GRIB SPECIFICATION ! KPDS(19) - VERSION NR OF PARAMETER TABLE ! ! KGDS(13) INTEGER*4 ! ARRAY CONTAINING GDS ELEMENTS. ! ! KGDS(1) - DATA REPRESENTATION TYPE ! ! LATITUDE/LONGITUDE GRIDS (SEE "GRIB" TABLE 10) ! KGDS(2) - N(I) NUMBER OF POINTS ON LATITUDE ! CIRCLE ! KGDS(3) - N(J) NUMBER OF POINTS ON LONGITUDE ! CIRCLE ! KGDS(4) - LA(1) LATITUDE OF ORIGIN ! KGDS(5) - LO(1) LONGITUDE OF ORIGIN ! KGDS(6) - RESOLUTION FLAG ! BIT MEANING ! 25 0 - DIRECTION INCREMENTS NOT ! GIVEN ! 1 - DIRECTION INCREMENTS GIVEN ! KGDS(7) - LA(2) LATITUDE OF EXTREME POINT ! KGDS(8) - LO(2) LONGITUDE OF EXTREME POINT ! KGDS(9) - DI LATITUDINAL DIRECTION INCREMENT ! KGDS(10) - REGULAR LAT/LON GRID ! DJ - LONGITUDINAL DIRECTION ! INCREMENT ! GAUSSIAN GRID ! N - NUMBER OF LATITUDE CIRCLES ! BETWEEN A POLE AND THE EQUATOR ! KGDS(11) - SCANNING MODE FLAG ! BIT MEANING ! 25 0 - POINTS ALONG A LATITUDE ! SCAN FROM WEST TO EAST ! 1 - POINTS ALONG A LATITUDE ! SCAN FROM EAST TO WEST ! 26 0 - POINTS ALONG A MERIDIAN ! SCAN FROM NORTH TO SOUTH ! 1 - POINTS ALONG A MERIDIAN ! SCAN FROM SOUTH TO NORTH ! 27 0 - POINTS SCAN FIRST ALONG ! CIRCLES OF LATITUDE, THEN ! ALONG MERIDIANS ! (FORTRAN: (I,J)) ! 1 - POINTS SCAN FIRST ALONG ! MERIDIANS THEN ALONG ! CIRCLES OF LATITUDE ! (FORTRAN: (J,I)) ! ! POLAR STEREOGRAPHIC GRIDS (SEE GRIB TABLE 12) ! KGDS(2) - N(I) NR POINTS ALONG LAT CIRCLE ! KGDS(3) - N(J) NR POINTS ALONG LON CIRCLE ! KGDS(4) - LA(1) LATITUDE OF ORIGIN ! KGDS(5) - LO(1) LONGITUDE OF ORIGIN ! KGDS(6) - RESERVED ! KGDS(7) - LOV GRID ORIENTATION ! KGDS(8) - DX - X DIRECTION INCREMENT ! KGDS(9) - DY - Y DIRECTION INCREMENT ! KGDS(10) - PROJECTION CENTER FLAG ! KGDS(11) - SCANNING MODE ! ! SPHERICAL HARMONIC COEFFICIENTS (SEE "GRIB" TABLE 14) ! KGDS(2) - J PENTAGONAL RESOLUTION PARAMETER ! KGDS(3) - K PENTAGONAL RESOLUTION PARAMETER ! KGDS(4) - M PENTAGONAL RESOLUTION PARAMETER ! KGDS(5) - REPRESENTATION TYPE ! KGDS(6) - COEFFICIENT STORAGE MODE ! ! MERCATOR GRIDS ! KGDS(2) - N(I) NR POINTS ON LATITUDE CIRCLE ! KGDS(3) - N(J) NR POINTS ON LONGITUDE MERIDIAN ! KGDS(4) - LA(1) LATITUDE OF ORIGIN ! KGDS(5) - LO(1) LONGITUDE OF ORIGIN ! KGDS(6) - RESOLUTION FLAG ! KGDS(7) - LA(2) LATITUDE OF LAST GRID POINT ! KGDS(8) - LO(2) LONGITUDE OF LAST GRID POINT ! KGDS(9) - LATIN - LATITUDE OF PROJECTION INTERSECTION ! KGDS(10) - RESERVED ! KGDS(11) - SCANNING MODE FLAG ! KGDS(12) - LONGITUDINAL DIR GRID LENGTH ! KGDS(13) - LATITUDINAL DIR GRID LENGTH ! LAMBERT CONFORMAL GRIDS ! KGDS(2) - NX NR POINTS ALONG X-AXIS ! KGDS(3) - NY NR POINTS ALONG Y-AXIS ! KGDS(4) - LA1 LAT OF ORIGIN (LOWER LEFT) ! KGDS(5) - LO1 LON OF ORIGIN (LOWER LEFT) ! KGDS(6) - RESOLUTION (RIGHT ADJ COPY OF OCTET 17) ! KGDS(7) - LOV - ORIENTATION OF GRID ! KGDS(8) - DX - X-DIR INCREMENT ! KGDS(9) - DY - Y-DIR INCREMENT ! KGDS(10) - PROJECTION CENTER FLAG ! KGDS(11) - SCANNING MODE FLAG ! KGDS(12) - LATIN 1 - FIRST LAT FROM POLE OF ! SECANT CONE INTERSECTION ! KGDS(13) - LATIN 2 - SECOND LAT FROM POLE OF ! SECANT CONE INTERSECTION ! ! LBMS(65160) LOGICAL ! ARRAY TO CONTAIN THE BIT MAP DESCRIBING THE ! PLACEMENT OF DATA IN THE OUTPUT ARRAY. IF A ! BIT MAP IS NOT INCLUDED IN THE SOURCE MESSAGE, ! ONE WILL BE GENERATED AUTOMATICALLY BY THE ! UNPACKING ROUTINE. ! ! ! DATA(65160) REAL*4 ! THIS ARRAY WILL CONTAIN THE UNPACKED DATA POINTS. ! ! NOTE:- 65160 IS MAXIMUN FIELD SIZE ALLOWABLE ! ! KPTR(10) INTEGER*4 ! ARRAY CONTAINING STORAGE FOR THE FOLLOWING ! PARAMETERS. ! ! (1) - UNUSED ! (2) - UNUSED ! (3) - LENGTH OF PDS (IN BYTES) ! (4) - LENGTH OF GDS (IN BYTES) ! (5) - LENGTH OF BMS (IN BYTES) ! (6) - LENGTH OF BDS (IN BYTES) ! (7) - USED BY UNPACKING ROUTINE ! (8) - NUMBER OF DATA POINTS FOR GRID ! (9) - "GRIB" CHARACTERS START IN BYTE NUMBER ! (10) - USED BY UNPACKING ROUTINE ! ! ! KRET INTEGER*4 ! THIS VARIABLE WILL CONTAIN THE RETURN INDICATOR. ! ! 0 - NO ERRORS DETECTED. ! ! 1 - 'GRIB' NOT FOUND IN FIRST 100 ! CHARACTERS. ! ! 2 - '7777' NOT FOUND, EITHER MISSING OR ! TOTAL OF SEC COUNTS OF INDIVIDUAL ! SECTIONS IS INCORRECT. ! ! 3 - UNPACKED FIELD IS LARGER THAN 65160. ! ! 4 - IN GDS, DATA REPRESENTATION TYPE ! NOT ONE OF THE CURRENTLY ACCEPTABLE ! VALUES. SEE "GRIB" TABLE 9. VALUE ! OF INCORRECT TYPE RETURNED IN KGDS(1). ! ! 5 - GRID INDICATED IN KPDS(3) IS NOT ! AVAILABLE FOR THE CENTER INDICATED IN ! KPDS(1) AND NO GDS SENT. ! ! 7 - EDITION INDICATED IN KPDS(18) HAS NOT ! YET BEEN INCLUDED IN THE DECODER. ! ! 8 - GRID IDENTIFICATION = 255 (NOT STANDARD ! GRID) BUT FLAG INDICATING PRESENCE OF ! GDS IS TURNED OFF. NO METHOD OF ! GENERATING PROPER GRID. ! ! 9 - PRODUCT OF KGDS(2) AND KGDS(3) DOES NOT ! MATCH STANDARD NUMBER OF POINTS FOR THIS ! GRID (FOR OTHER THAN SPECTRALS). THIS ! WILL OCCUR ONLY IF THE GRID. ! IDENTIFICATION, KPDS(3), AND A ! TRANSMITTED GDS ARE INCONSISTENT. ! ! 10 - CENTER INDICATOR WAS NOT ONE INDICATED ! IN "GRIB" TABLE 1. PLEASE CONTACT AD ! PRODUCTION MANAGEMENT BRANCH (W/NMC42) ! IF THIS ERROR IS ENCOUNTERED. ! ! 11 - BINARY DATA SECTION (BDS) NOT COMPLETELY ! PROCESSED. PROGRAM IS NOT SET TO PROCESS ! FLAG COMBINATIONS AS SHOWN IN ! OCTETS 4 AND 14. ! ! ! LIST OF TEXT MESSAGES FROM CODE ! ! ! W3FI63/FI632 ! ! 'HAVE ENCOUNTERED A NEW GRID FOR NMC, PLEASE NOTIFY ! AUTOMATION DIVISION, PRODUCTION MANAGEMENT BRANCH ! (W/NMC42)' ! ! 'HAVE ENCOUNTERED A NEW GRID FOR ECMWF, PLEASE NOTIFY ! AUTOMATION DIVISION, PRODUCTION MANAGEMENT BRANCH ! (W/NMC42)' ! ! 'HAVE ENCOUNTERED A NEW GRID FOR U.K. METEOROLOGICAL ! OFFICE, BRACKNELL. PLEASE NOTIFY AUTOMATION DIVISION, ! PRODUCTION MANAGEMENT BRANCH (W/NMC42)' ! ! 'HAVE ENCOUNTERED A NEW GRID FOR FNOC, PLEASE NOTIFY ! AUTOMATION DIVISION, PRODUCTION MANAGEMENT BRANCH ! (W/NMC42)' ! ! ! W3FI63/FI633 ! ! 'POLAR STEREO PROCESSING NOT AVAILABLE' * ! ! W3FI63/FI634 ! ! 'WARNING - BIT MAP MAY NOT BE ASSOCIATED WITH SPHERICAL ! COEFFICIENTS' ! ! ! W3FI63/FI637 ! ! 'NO CURRENT LISTING OF FNOC GRIDS' * ! ! ! * WILL BE AVAILABLE IN NEXT UPDATE ! *************************************************************** ! ! INCOMING MESSAGE HOLDER CHARACTER :: msga(*) ! BIT MAP LOGICAL*1 :: kbms(*) ! ! ELEMENTS OF PRODUCT DESCRIPTION SEC (PDS) INTEGER :: kpds(*) ! ELEMENTS OF GRID DESCRIPTION SEC (PDS) INTEGER :: kgds(*) ! ! CONTAINER FOR GRIB GRID REAL :: DATA(*) ! ! ARRAY OF POINTERS AND COUNTERS INTEGER :: kptr(*) ! ! ***************************************************************** INTEGER :: kkk,jsgn,jexp,ifr,npts CHARACTER (LEN=1) :: kk(8) REAL :: realkk,fval1,fdiff1 EQUIVALENCE (kk(1),kkk) ! ***************************************************************** ! 1.0 LOCATE BEGINNING OF 'GRIB' MESSAGE ! FIND 'GRIB' CHARACTERS ! 2.0 USE COUNTS IN EACH DESCRIPTION SEC TO DETERMINE ! IF '7777' IS IN PROPER PLACE. ! 3.0 PARSE PRODUCT DEFINITION SECTION. ! 4.0 PARSE GRID DESCRIPTION SEC (IF INCLUDED) ! 5.0 PARSE BIT MAP SEC (IF INCLUDED) ! 6.0 USING INFORMATION FROM PRODUCT DEFINITION, GRID ! DESCRIPTION, AND BIT MAP SECTIONS.. EXTRACT ! DATA AND PLACE INTO PROPER ARRAY. ! ******************************************************************* ! ! MAIN DRIVER ! ! ******************************************************************* SAVE kptr(10) = 0 ! SEE IF PROPER 'GRIB' KEY EXISTS, THEN ! USING SEC COUNTS, DETERMINE IF '7777' ! IS IN THE PROPER LOCATION ! CALL fi631(msga,kptr,kpds,kret) IF(kret /= 0) THEN GO TO 900 END IF ! PRINT *,'FI631 KPTR',(KPTR(I),I=1,16) ! ! PARSE PARAMETERS FROM PRODUCT DESCRIPTION SECTION ! CALL fi632(msga,kptr,kpds,kret) IF(kret /= 0) THEN GO TO 900 END IF ! PRINT *,'FI632 KPTR',(KPTR(I),I=1,16) ! ! IF AVAILABLE, EXTRACT NEW GRID DESCRIPTION ! IF (IAND(kpds(4),128) /= 0) THEN CALL fi633(msga,kptr,kgds,kret) IF(kret /= 0) THEN GO TO 900 END IF ! PRINT *,'FI633 KPTR',(KPTR(I),I=1,16) END IF ! ! EXTRACT OR GENERATE BIT MAP ! CALL fi634(msga,kptr,kpds,kgds,kbms,kret) IF(kret /= 0) THEN GO TO 900 END IF ! PRINT *,'FI634 KPTR',(KPTR(I),I=1,16) ! ! USING INFORMATION FROM PDS, BMS AND BIT DATA SEC , ! EXTRACT AND SAVE IN GRIB GRID, ALL DATA ENTRIES. ! IF (kpds(18) == 1) THEN CALL fi635(msga,kptr,kpds,kgds,kbms,DATA,kret) IF (kptr(3) >= 50) THEN ! ! PDS GREATER THAN 28 BYTES ! THEREFORE SOMETHING SPECIAL IS GOING ON ! ! IN THIS CASE 2ND DIFFERENCE PACKING ! NEEDS TO BE UNDONE. ! ! EXTRACT FIRST VALUE FROM BYTE 41-44 PDS ! KPTR(9) CONTAINS OFFSET TO START OF ! GRIB MESSAGE. ! EXTRACT FIRST FIRST-DIFFERENCE FROM BYTES 45-48 PDS ! ! AND EXTRACT SCALE FACTOR (E) TO UNDO 2**E ! THAT WAS APPLIED PRIOR TO 2ND ORDER PACKING ! AND PLACED IN PDS BYTES 49-51 ! FACTOR IS A SIGNED TWO BYTE INTEGER ! ! ALSO NEED THE DECIMAL SCALING FROM PDS(27-28) ! (AVAILABLE IN KPDS(22) FROM UNPACKER) ! TO UNDO THE DECIMAL SCALING APPLIED TO THE ! SECOND DIFFERENCES DURING UNPACKING. ! SECOND DIFFS ALWAYS PACKED WITH 0 DECIMAL SCALE ! BUT UNPACKER DOESNT KNOW THAT. ! ! CALL GRBGBYTE (MSGA,FVAL1,KPTR(9)+384,32) ! ! NOTE INTEGERS, CHARACTERS AND EQUIVALENCES ! DEFINED ABOVE TO MAKE THIS KKK EXTRACTION ! WORK AND LINE UP ON WORD BOUNDARIES ! CALL grbgbyte (msga,kkk,kptr(9)+384,32) ! ! THE NEXT CODE WILL CONVERT THE IBM370 FOATING POINT ! TO THE FLOATING POINT USED ON YOUR MACHINE. ! ! 1ST TEST TO SEE IN ON 32 OR 64 BIT WORD MACHINE ! LW = 4 OR 8; IF 8 MAY BE A CRAY ! CALL w3fi01(lw) IF (lw == 4) THEN CALL grbgbyte (kk,jsgn,0,1) CALL grbgbyte (kk,jexp,1,7) CALL grbgbyte (kk,ifr,8,24) ELSE CALL grbgbyte (kk,jsgn,32,1) CALL grbgbyte (kk,jexp,33,7) CALL grbgbyte (kk,ifr,40,24) END IF ! IF (ifr == 0) THEN realkk = 0.0 ELSE IF (jexp == 0.AND.ifr == 0) THEN realkk = 0.0 ELSE realkk = FLOAT(ifr) * 16.0 ** (jexp - 64 - 6) IF (jsgn /= 0) realkk = -realkk END IF fval1 = realkk ! ! CALL GRBGBYTE (MSGA,FDIFF1,KPTR(9)+416,32) ! (REPLACED BY FOLLOWING EXTRACTION) ! CALL grbgbyte (msga,kkk,kptr(9)+416,32) ! ! THE NEXT CODE WILL CONVERT THE IBM370 FOATING POINT ! TO THE FLOATING POINT USED ON YOUR MACHINE. ! ! 1ST TEST TO SEE IN ON 32 OR 64 BIT WORD MACHINE ! LW = 4 OR 8; IF 8 MAY BE A CRAY ! CALL w3fi01(lw) IF (lw == 4) THEN CALL grbgbyte (kk,jsgn,0,1) CALL grbgbyte (kk,jexp,1,7) CALL grbgbyte (kk,ifr,8,24) ELSE CALL grbgbyte (kk,jsgn,32,1) CALL grbgbyte (kk,jexp,33,7) CALL grbgbyte (kk,ifr,40,24) END IF ! IF (ifr == 0) THEN realkk = 0.0 ELSE IF (jexp == 0.AND.ifr == 0) THEN realkk = 0.0 ELSE realkk = FLOAT(ifr) * 16.0 ** (jexp - 64 - 6) IF (jsgn /= 0) realkk = -realkk END IF fdiff1 = realkk ! CALL grbgbyte (msga,ISIGN,kptr(9)+448,1) CALL grbgbyte (msga,iscal2,kptr(9)+449,15) IF(ISIGN > 0) THEN iscal2 = - iscal2 END IF ! PRINT *,'DELTA POINT 1-',FVAL1 ! PRINT *,'DELTA POINT 2-',FDIFF1 ! PRINT *,'DELTA POINT 3-',ISCAL2 npts = kptr(10) ! WRITE (6,FMT='('' 2ND DIFF POINTS IN FIELD = '',/, ! & 10(3X,10F12.2,/))') (DATA(I),I=1,NPTS) ! PRINT *,'DELTA POINT 4-',KPDS(22) CALL w3fi83 (DATA,npts,fval1,fdiff1, & iscal2,kpds(22),kpds,kgds) ! WRITE (6,FMT='('' 2ND DIFF EXPANDED POINTS IN FIELD = '', ! & /,10(3X,10F12.2,/))') (DATA(I),I=1,NPTS) ! WRITE (6,FMT='('' END OF ARRAY IN FIELD = '',/, ! & 10(3X,10F12.2,/))') (DATA(I),I=NPTS-5,NPTS) END IF ELSE PRINT *,'FI635 NOT PROGRAMMED FOR EDITION NR',kpds(18) kret = 7 END IF ! 900 RETURN END SUBROUTINE w3fi63 SUBROUTINE fi631(msga,kptr,kpds,kret) 1,10 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: FI631 FIND 'GRIB' CHARS & RESET POINTERS ! PRGMMR: BILL CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 ! ! ABSTRACT: FIND 'GRIB; CHARACTERS AND SET POINTERS TO THE NEXT ! BYTE FOLLOWING 'GRIB'. IF THEY EXIST EXTRACT COUNTS FROM GDS AND ! BMS. EXTRACT COUNT FROM BDS. DETERMINE IF SUM OF COUNTS ACTUALLY ! PLACES TERMINATOR '7777' AT THE CORRECT LOCATION. ! ! PROGRAM HISTORY LOG: ! 91-09-13 CAVANAUGH ! ! USAGE: CALL FI631(MSGA,KPTR,KPDS,KRET) ! INPUT ARGUMENT LIST: ! MSGA - GRIB FIELD - "GRIB" THRU "7777" ! KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS ! (1) - TOTAL LENGTH OF GRIB MESSAGE ! (2) - LENGTH OF INDICATOR (SECTION 0) ! (3) - LENGTH OF PDS (SECTION 1) ! (4) - LENGTH OF GDS (SECTION 2) ! (5) - LENGTH OF BMS (SECTION 3) ! (6) - LENGTH OF BDS (SECTION 4) ! (7) - VALUE OF CURRENT BYTE ! (8) - BIT POINTER ! (9) - GRIB START BIT NR ! (10) - GRIB/GRID ELEMENT COUNT ! (11) - NR UNUSED BITS AT END OF SECTION 3 ! (12) - BIT MAP FLAG ! (13) - NR UNUSED BITS AT END OF SECTION 2 ! (14) - BDS FLAGS ! (15) - NR UNUSED BITS AT END OF SECTION 4 ! ! OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) ! KPDS - ARRAY CONTAINING PDS ELEMENTS. ! (1) - ID OF CENTER ! (2) - MODEL IDENTIFICATION ! (3) - GRID IDENTIFICATION ! (4) - GDS/BMS FLAG ! (5) - INDICATOR OF PARAMETER ! (6) - TYPE OF LEVEL ! (7) - HEIGHT/PRESSURE , ETC OF LEVEL ! (8) - YEAR OF CENTURY ! (9) - MONTH OF YEAR ! (10) - DAY OF MONTH ! (11) - HOUR OF DAY ! (12) - MINUTE OF HOUR ! (13) - INDICATOR OF FORECAST TIME UNIT ! (14) - TIME RANGE 1 ! (15) - TIME RANGE 2 ! (16) - TIME RANGE FLAG ! (17) - NUMBER INCLUDED IN AVERAGE ! KPTR - SEE INPUT LIST ! KRET - ERROR RETURN ! ! REMARKS: ! ERROR RETURNS ! KRET = 1 - NO 'GRIB' ! 2 - NO '7777' OR MISLOCATED (BY COUNTS) ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 77 ! MACHINE: NAS9000 ! !$$$ ! ! INCOMING MESSAGE HOLDER CHARACTER :: msga(*) ! ARRAY OF POINTERS AND COUNTERS INTEGER :: kptr(*) ! PRODUCT DESCRIPTION SECTION DATA. INTEGER :: kpds(*) ! INTEGER :: kret ! ! ****************************************************************** SAVE kret = 0 ! ------------------- FIND 'GRIB' KEY DO i = 0, 839, 8 CALL grbgbyte (msga,mgrib,i,32) IF (mgrib == 1196575042) THEN kptr(9) = i GO TO 60 END IF END DO kret = 1 RETURN 60 CONTINUE ! -------------FOUND 'GRIB' ! SKIP GRIB CHARACTERS ! PRINT *,'FI631 GRIB AT',I kptr(8) = kptr(9) + 32 CALL grbgbyte (msga,itotal,kptr(8),24) ! HAVE LIFTED WHAT MAY BE A MSG TOTAL BYTE COUNT ipoint = kptr(9) + itotal * 8 - 32 CALL grbgbyte (msga,i7777,ipoint,32) IF (i7777 == 926365495) THEN ! HAVE FOUND END OF MESSAGE '7777' IN PROPER LOCATION ! MARK AND PROCESS AS GRIB VERSION 1 OR HIGHER ! PRINT *,'FI631 7777 AT',IPOINT kptr(8) = kptr(8) + 24 kptr(1) = itotal kptr(2) = 8 CALL grbgbyte (msga,kpds(18),kptr(8),8) kptr(8) = kptr(8) + 8 ELSE ! CANNOT FIND END OF GRIB EDITION 1 MESSAGE kret = 2 RETURN END IF ! ------------------- PROCESS SECTION 1 ! EXTRACT COUNT FROM PDS ! PRINT *,'START OF PDS',KPTR(8) CALL grbgbyte (msga,kptr(3),kptr(8),24) look = kptr(8) + 56 ! EXTRACT GDS/BMS FLAG CALL grbgbyte (msga,kpds(4),look,8) kptr(8) = kptr(8) + kptr(3) * 8 ! PRINT *,'START OF GDS',KPTR(8) IF (IAND(kpds(4),128) /= 0) THEN ! EXTRACT COUNT FROM GDS CALL grbgbyte (msga,kptr(4),kptr(8),24) kptr(8) = kptr(8) + kptr(4) * 8 ELSE kptr(4) = 0 END IF ! PRINT *,'START OF BMS',KPTR(8) IF (IAND(kpds(4),64) /= 0) THEN ! EXTRACT COUNT FROM BMS CALL grbgbyte (msga,kptr(5),kptr(8),24) ELSE kptr(5) = 0 END IF kptr(8) = kptr(8) + kptr(5) * 8 ! PRINT *,'START OF BDS',KPTR(8) ! EXTRACT COUNT FROM BDS CALL grbgbyte (msga,kptr(6),kptr(8),24) ! --------------- TEST FOR '7777' ! PRINT *,(KPTR(KJ),KJ=1,10) kptr(8) = kptr(8) + kptr(6) * 8 ! EXTRACT FOUR BYTES FROM THIS LOCATION ! PRINT *,'FI631 LOOKING FOR 7777 AT',KPTR(8) CALL grbgbyte (msga,k7777,kptr(8),32) match = kptr(2) + kptr(3) + kptr(4) + kptr(5) + kptr(6) + 4 IF (k7777 /= 926365495.OR.match /= kptr(1)) THEN kret = 2 ELSE ! PRINT *,'FI631 7777 AT',KPTR(8) IF (kpds(18) == 0) THEN kptr(1) = kptr(2) + kptr(3) + kptr(4) + kptr(5) + & kptr(6) + 4 END IF END IF ! PRINT *,'KPTR',(KPTR(I),I=1,16) RETURN END SUBROUTINE fi631 SUBROUTINE fi632(msga,kptr,kpds,kret) 1,24 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: FI632 GATHER INFO FROM PRODUCT DEFINITION SEC ! PRGMMR: BILL CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 ! ! ABSTRACT: EXTRACT INFORMATION FROM THE PRODUCT DESCRIPTION ! SEC , AND GENERATE LABEL INFORMATION TO PERMIT STORAGE ! IN OFFICE NOTE 84 FORMAT. ! ! PROGRAM HISTORY LOG: ! 91-09-13 CAVANAUGH ! 93-12-08 CAVANAUGH CORRECTED TEST FOR EDITION NUMBER INSTEAD ! OF VERSION NUMBER ! ! USAGE: CALL FI632(MSGA,KPTR,KPDS,KRET) ! INPUT ARGUMENT LIST: ! MSGA - ARRAY CONTAINING GRIB MESSAGE ! KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS ! (1) - TOTAL LENGTH OF GRIB MESSAGE ! (2) - LENGTH OF INDICATOR (SECTION 0) ! (3) - LENGTH OF PDS (SECTION 1) ! (4) - LENGTH OF GDS (SECTION 2) ! (5) - LENGTH OF BMS (SECTION 3) ! (6) - LENGTH OF BDS (SECTION 4) ! (7) - VALUE OF CURRENT BYTE ! (8) - BIT POINTER ! (9) - GRIB START BIT NR ! (10) - GRIB/GRID ELEMENT COUNT ! (11) - NR UNUSED BITS AT END OF SECTION 3 ! (12) - BIT MAP FLAG ! (13) - NR UNUSED BITS AT END OF SECTION 2 ! (14) - BDS FLAGS ! (15) - NR UNUSED BITS AT END OF SECTION 4 ! ! OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) ! KPDS - ARRAY CONTAINING PDS ELEMENTS. ! (1) - ID OF CENTER ! (2) - MODEL IDENTIFICATION ! (3) - GRID IDENTIFICATION ! (4) - GDS/BMS FLAG ! (5) - INDICATOR OF PARAMETER ! (6) - TYPE OF LEVEL ! (7) - HEIGHT/PRESSURE , ETC OF LEVEL ! (8) - YEAR OF CENTURY ! (9) - MONTH OF YEAR ! (10) - DAY OF MONTH ! (11) - HOUR OF DAY ! (12) - MINUTE OF HOUR ! (13) - INDICATOR OF FORECAST TIME UNIT ! (14) - TIME RANGE 1 ! (15) - TIME RANGE 2 ! (16) - TIME RANGE FLAG ! (17) - NUMBER INCLUDED IN AVERAGE ! (18) - ! (19) - ! (20) - NUMBER MISSING FROM AVGS/ACCUMULATIONS ! (21) - CENTURY ! (22) - UNITS DECIMAL SCALE FACTOR ! (23) - SUBCENTER ! KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS ! SEE INPUT LIST ! KRET - ERROR RETURN ! ! REMARKS: ! ERROR RETURN = 0 - NO ERRORS ! = 8 - TEMP GDS INDICATED, BUT NO GDS ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 77 ! MACHINE: NAS9000 ! !$$$ ! ! INCOMING MESSAGE HOLDER CHARACTER :: msga(*) ! ! ARRAY OF POINTERS AND COUNTERS INTEGER :: kptr(*) ! PRODUCT DESCRIPTION SECTION ENTRIES INTEGER :: kpds(*) ! INTEGER :: kret ! ------------------- PROCESS SECTION 1 SAVE kptr(8) = kptr(9) + kptr(2) * 8 + 24 ! BYTE 4 ! PARAMETER TABLE VERSION NR CALL grbgbyte (msga,kpds(19),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 5 IDENTIFICATION OF CENTER CALL grbgbyte (msga,kpds(1),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 6 ! GET GENERATING PROCESS ID NR CALL grbgbyte (msga,kpds(2),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 7 ! GRID DEFINITION CALL grbgbyte (msga,kpds(3),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 8 ! GDS/BMS FLAGS ! CALL GRBGBYTE (MSGA,KPDS(4),KPTR(8),8) kptr(8) = kptr(8) + 8 ! BYTE 9 ! INDICATOR OF PARAMETER CALL grbgbyte (msga,kpds(5),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 10 ! TYPE OF LEVEL CALL grbgbyte (msga,kpds(6),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 11,12 ! HEIGHT/PRESSURE CALL grbgbyte (msga,kpds(7),kptr(8),16) kptr(8) = kptr(8) + 16 ! BYTE 13 ! YEAR OF CENTURY CALL grbgbyte (msga,kpds(8),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 14 ! MONTH OF YEAR CALL grbgbyte (msga,kpds(9),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 15 ! DAY OF MONTH CALL grbgbyte (msga,kpds(10),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 16 ! HOUR OF DAY CALL grbgbyte (msga,kpds(11),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 17 ! MINUTE CALL grbgbyte (msga,kpds(12),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 18 ! INDICATOR TIME UNIT RANGE CALL grbgbyte (msga,kpds(13),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 19 ! P1 - PERIOD OF TIME CALL grbgbyte (msga,kpds(14),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 20 ! P2 - PERIOD OF TIME CALL grbgbyte (msga,kpds(15),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 21 ! TIME RANGE INDICATOR CALL grbgbyte (msga,kpds(16),kptr(8),8) kptr(8) = kptr(8) + 8 ! ! IF TIME RANGE INDICATOR IS 10, P1 IS PACKED IN ! PDS BYTES 19-20 ! IF (kpds(16) == 10) THEN kpds(14) = kpds(14) * 256 + kpds(15) kpds(15) = 0 END IF ! BYTE 22,23 ! NUMBER INCLUDED IN AVERAGE CALL grbgbyte (msga,kpds(17),kptr(8),16) kptr(8) = kptr(8) + 16 ! BYTE 24 ! NUMBER MISSING FROM AVERAGES/ACCUMULATIONS CALL grbgbyte (msga,kpds(20),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 25 ! IDENTIFICATION OF CENTURY CALL grbgbyte (msga,kpds(21),kptr(8),8) kptr(8) = kptr(8) + 8 IF (kptr(3) > 25) THEN ! BYTE 26 SUB CENTER NUMBER CALL grbgbyte (msga,kpds(23),kptr(8),8) kptr(8) = kptr(8) + 8 IF (kptr(3) >= 28) THEN ! BYTE 27-28 ! UNITS DECIMAL SCALE FACTOR CALL grbgbyte (msga,ISIGN,kptr(8),1) kptr(8) = kptr(8) + 1 CALL grbgbyte (msga,idec,kptr(8),15) kptr(8) = kptr(8) + 15 IF (ISIGN > 0) THEN kpds(22) = - idec ELSE kpds(22) = idec END IF isiz = kptr(3) - 28 IF (isiz <= 12) THEN ! BYTES 29-40 CURRENTLY RESERVED FOR FUTURE USE kptr(8) = kptr(8) + isiz * 8 ELSE ! BYTES 29-40 CURRENTLY RESERVED FOR FUTURE USE kptr(8) = kptr(8) + 12 * 8 ! BYTES 40 - N LOCAL USE DATA CALL w3fi01(lw) mwdbit = lw * 8 isiz = kptr(3) - 40 iter = isiz / lw CALL grbgbytes (msga,kpds(36),kptr(8),mwdbit,0,iter) kptr(8) = kptr(8) + isiz * 8 END IF END IF END IF ! ----------- TEST FOR NEW GRID IF (IAND(kpds(4),128) /= 0) THEN IF (IAND(kpds(4),64) /= 0) THEN IF (kpds(3) /= 255) THEN IF (kpds(3) >= 21.AND.kpds(3) <= 26)THEN RETURN ELSE IF (kpds(3) >= 37.AND.kpds(3) <= 44)THEN RETURN ELSE IF (kpds(3) >= 61.AND.kpds(3) <= 64) THEN RETURN END IF IF (kpds(1) == 7) THEN IF (kpds(3) >= 2.AND.kpds(3) <= 3)THEN ELSE IF (kpds(3) >= 5.AND.kpds(3) <= 6)THEN ELSE IF (kpds(3) >= 27.AND.kpds(3) <= 34)THEN ELSE IF (kpds(3) == 50) THEN ELSE IF (kpds(3) >= 70.AND.kpds(3) <= 77) THEN ELSE IF (kpds(3) >= 100.AND.kpds(3) <= 105) THEN ELSE IF (kpds(3) >= 201.AND.kpds(3) <= 214) THEN ELSE PRINT *,' HAVE ENCOUNTERED A NEW GRID FOR', & ' NMC WITHOUT A GRID DESCRIPTION SECTION' PRINT *,' PLEASE NOTIFY AUTOMATION DIVISION' PRINT *,' PRODUCTION MANAGEMENT BRANCH' PRINT *,' W/NMC42)' END IF ELSE IF (kpds(1) == 98) THEN IF (kpds(3) >= 1.AND.kpds(3) <= 16) THEN ELSE PRINT *,' HAVE ENCOUNTERED A NEW GRID FOR', & ' ECMWF WITHOUT A GRID DESCRIPTION SECTION' PRINT *,' PLEASE NOTIFY AUTOMATION DIVISION' PRINT *,' PRODUCTION MANAGEMENT BRANCH' PRINT *,' W/NMC42)' END IF ELSE IF (kpds(1) == 74) THEN IF (kpds(3) >= 1.AND.kpds(3) <= 12) THEN ELSE IF (kpds(3) >= 21.AND.kpds(3) <= 26)THEN ELSE IF (kpds(3) >= 61.AND.kpds(3) <= 64) THEN ELSE IF (kpds(3) >= 70.AND.kpds(3) <= 77) THEN ELSE PRINT *,' HAVE ENCOUNTERED A NEW GRID FOR', & ' U.K. MET OFFICE, BRACKNELL', & ' WITHOUT A GRID DESCRIPTION SECTION' PRINT *,' PLEASE NOTIFY AUTOMATION DIVISION' PRINT *,' PRODUCTION MANAGEMENT BRANCH' PRINT *,' W/NMC42)' END IF ELSE IF (kpds(1) == 58) THEN IF (kpds(3) >= 1.AND.kpds(3) <= 12) THEN ELSE PRINT *,' HAVE ENCOUNTERED A NEW GRID FOR', & ' FNOC WITHOUT A GRID DESCRIPTION SECTION' PRINT *,' PLEASE NOTIFY AUTOMATION DIVISION' PRINT *,' PRODUCTION MANAGEMENT BRANCH' PRINT *,' W/NMC42)' END IF END IF END IF END IF END IF RETURN END SUBROUTINE fi632 SUBROUTINE fi633(msga,kptr,kgds,kret) 1,58 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: FI633 EXTRACT INFO FROM GRIB-GDS ! PRGMMR: BILL CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 ! ! ABSTRACT: EXTRACT INFORMATION ON UNLISTED GRID TO ALLOW ! CONVERSION TO OFFICE NOTE 84 FORMAT. ! ! PROGRAM HISTORY LOG: ! 91-09-13 CAVANAUGH ! ! USAGE: CALL FI633(MSGA,KPTR,KGDS,KRET) ! INPUT ARGUMENT LIST: ! MSGA - ARRAY CONTAINING GRIB MESSAGE ! KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS ! (1) - TOTAL LENGTH OF GRIB MESSAGE ! (2) - LENGTH OF INDICATOR (SECTION 0) ! (3) - LENGTH OF PDS (SECTION 1) ! (4) - LENGTH OF GDS (SECTION 2) ! (5) - LENGTH OF BMS (SECTION 3) ! (6) - LENGTH OF BDS (SECTION 4) ! (7) - VALUE OF CURRENT BYTE ! (8) - BIT POINTER ! (9) - GRIB START BIT NR ! (10) - GRIB/GRID ELEMENT COUNT ! (11) - NR UNUSED BITS AT END OF SECTION 3 ! (12) - BIT MAP FLAG ! (13) - NR UNUSED BITS AT END OF SECTION 2 ! (14) - BDS FLAGS ! (15) - NR UNUSED BITS AT END OF SECTION 4 ! ! OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) ! KGDS - ARRAY CONTAINING GDS ELEMENTS. ! (1) - DATA REPRESENTATION TYPE ! (19) - NUMBER OF VERTICAL COORDINATE PARAMETERS ! (20) - OCTET NUMBER OF THE LIST OF VERTICAL COORDINATE ! PARAMETERS ! OR ! OCTET NUMBER OF THE LIST OF NUMBERS OF POINTS ! IN EACH ROW ! OR ! 255 IF NEITHER ARE PRESENT ! (21) - FOR GRIDS WITH PL, NUMBER OF POINTS IN GRID ! (22) - NUMBER OF WORDS IN EACH ROW ! LATITUDE/LONGITUDE GRIDS ! (2) - N(I) NR POINTS ON LATITUDE CIRCLE ! (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN ! (4) - LA(1) LATITUDE OF ORIGIN ! (5) - LO(1) LONGITUDE OF ORIGIN ! (6) - RESOLUTION FLAG ! (7) - LA(2) LATITUDE OF EXTREME POINT ! (8) - LO(2) LONGITUDE OF EXTREME POINT ! (9) - DI LATITUDINAL DIRECTION OF INCREMENT ! (10) - DJ LONGITUDINAL DIRECTION INCREMENT ! (11) - SCANNING MODE FLAG ! POLAR STEREOGRAPHIC GRIDS ! (2) - N(I) NR POINTS ALONG LAT CIRCLE ! (3) - N(J) NR POINTS ALONG LON CIRCLE ! (4) - LA(1) LATITUDE OF ORIGIN ! (5) - LO(1) LONGITUDE OF ORIGIN ! (6) - RESERVED ! (7) - LOV GRID ORIENTATION ! (8) - DX - X DIRECTION INCREMENT ! (9) - DY - Y DIRECTION INCREMENT ! (10) - PROJECTION CENTER FLAG ! (11) - SCANNING MODE ! SPHERICAL HARMONIC COEFFICIENTS ! (2) - J PENTAGONAL RESOLUTION PARAMETER ! (3) - K " " " ! (4) - M " " " ! (5) - REPRESENTATION TYPE ! (6) - COEFFICIENT STORAGE MODE ! MERCATOR GRIDS ! (2) - N(I) NR POINTS ON LATITUDE CIRCLE ! (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN ! (4) - LA(1) LATITUDE OF ORIGIN ! (5) - LO(1) LONGITUDE OF ORIGIN ! (6) - RESOLUTION FLAG ! (7) - LA(2) LATITUDE OF LAST GRID POINT ! (8) - LO(2) LONGITUDE OF LAST GRID POINT ! (9) - LATIN - LATITUDE OF PROJECTION INTERSECTION ! (10) - RESERVED ! (11) - SCANNING MODE FLAG ! (12) - LONGITUDINAL DIR GRID LENGTH ! (13) - LATITUDINAL DIR GRID LENGTH ! LAMBERT CONFORMAL GRIDS ! (2) - NX NR POINTS ALONG X-AXIS ! (3) - NY NR POINTS ALONG Y-AXIS ! (4) - LA1 LAT OF ORIGIN (LOWER LEFT) ! (5) - LO1 LON OF ORIGIN (LOWER LEFT) ! (6) - RESOLUTION (RIGHT ADJ COPY OF OCTET 17) ! (7) - LOV - ORIENTATION OF GRID ! (8) - DX - X-DIR INCREMENT ! (9) - DY - Y-DIR INCREMENT ! (10) - PROJECTION CENTER FLAG ! (11) - SCANNING MODE FLAG ! (12) - LATIN 1 - FIRST LAT FROM POLE OF SECANT CONE INTER ! (13) - LATIN 2 - SECOND LAT FROM POLE OF SECANT CONE INTER ! KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS ! SEE INPUT LIST ! KRET - ERROR RETURN ! ! REMARKS: ! KRET = 0 ! = 4 - DATA REPRESENTATION TYPE NOT CURRENTLY ACCEPTABLE ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 77 ! MACHINE: NAS9000 ! !$$$ ! ************************************************************ ! INCOMING MESSAGE HOLDER CHARACTER :: msga(*) ! ! ARRAY GDS ELEMENTS INTEGER :: kgds(*) ! ARRAY OF POINTERS AND COUNTERS INTEGER :: kptr(*) ! INTEGER :: kret ! --------------------------------------------------------------- SAVE kret = 0 ! PROCESS GRID DEFINITION SECTION (IF PRESENT) ! MAKE SURE BIT POINTER IS PROPERLY SET kptr(8) = kptr(9) + (kptr(2)*8) + (kptr(3)*8) + 24 nsave = kptr(8) - 24 ! BYTE 4 ! NV - NR OF VERT COORD PARAMETERS CALL grbgbyte (msga,kgds(19),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 5 ! PV - LOCATION - SEE FM92 MANUAL CALL grbgbyte (msga,kgds(20),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTE 6 ! DATA REPRESENTATION TYPE CALL grbgbyte (msga,kgds(1),kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTES 7-32 ARE GRID DEFINITION DEPENDING ON ! DATA REPRESENTATION TYPE IF (kgds(1) == 0) THEN GO TO 1000 ELSE IF (kgds(1) == 1) THEN GO TO 4000 ELSE IF (kgds(1) == 2.OR.kgds(1) == 5) THEN GO TO 2000 ELSE IF (kgds(1) == 3) THEN GO TO 5000 ELSE IF (kgds(1) == 4) THEN GO TO 1000 ! ELSE IF (KGDS(1).EQ.10) THEN ! ELSE IF (KGDS(1).EQ.14) THEN ! ELSE IF (KGDS(1).EQ.20) THEN ! ELSE IF (KGDS(1).EQ.24) THEN ! ELSE IF (KGDS(1).EQ.30) THEN ! ELSE IF (KGDS(1).EQ.34) THEN ELSE IF (kgds(1) == 50) THEN GO TO 3000 ! ELSE IF (KGDS(1).EQ.60) THEN ! ELSE IF (KGDS(1).EQ.70) THEN ! ELSE IF (KGDS(1).EQ.80) THEN ELSE ! MARK AS GDS/ UNKNOWN DATA REPRESENTATION TYPE kret = 4 RETURN END IF ! BYTE 33-N VERTICAL COORDINATE PARAMETERS ! ----------- ! BYTES 33-42 EXTENSIONS OF GRID DEFINITION FOR ROTATION ! OR STRETCHING OF THE COORDINATE SYSTEM OR ! LAMBERT CONFORMAL PROJECTION. ! BYTE 43-N VERTICAL COORDINATE PARAMETERS ! ----------- ! BYTES 33-52 EXTENSIONS OF GRID DEFINITION FOR STRETCHED ! AND ROTATED COORDINATE SYSTEM ! BYTE 53-N VERTICAL COORDINATE PARAMETERS ! ----------- ! ************************************************************ ! ------------------- LATITUDE/LONGITUDE GRIDS ! ! ------------------- BYTE 7-8 NR OF POINTS ALONG LATITUDE CIRCLE 1000 CONTINUE CALL grbgbyte (msga,kgds(2),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 9-10 NR OF POINTS ALONG LONG MERIDIAN CALL grbgbyte (msga,kgds(3),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 11-13 LATITUDE OF ORIGIN CALL grbgbyte (msga,kgds(4),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(4),8388608) /= 0) THEN kgds(4) = IAND(kgds(4),8388607) * (-1) END IF ! ------------------- BYTE 14-16 LONGITUDE OF ORIGIN CALL grbgbyte (msga,kgds(5),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(5),8388608) /= 0) THEN kgds(5) = - IAND(kgds(5),8388607) END IF ! ------------------- BYTE 17 RESOLUTION FLAG CALL grbgbyte (msga,kgds(6),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 18-20 LATITUDE OF LAST GRID POINT CALL grbgbyte (msga,kgds(7),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(7),8388608) /= 0) THEN kgds(7) = - IAND(kgds(7),8388607) END IF ! ------------------- BYTE 21-23 LONGITUDE OF LAST GRID POINT CALL grbgbyte (msga,kgds(8),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(8),8388608) /= 0) THEN kgds(8) = - IAND(kgds(8),8388607) END IF ! ------------------- BYTE 24-25 LATITUDINAL DIR INCREMENT CALL grbgbyte (msga,kgds(9),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 26-27 IF REGULAR LAT/LON GRID ! HAVE LONGIT DIR INCREMENT ! ELSE IF GAUSSIAN GRID ! HAVE NR OF LAT CIRCLES ! BETWEEN POLE AND EQUATOR CALL grbgbyte (msga,kgds(10),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 28 SCANNING MODE FLAGS CALL grbgbyte (msga,kgds(11),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 29-32 RESERVED ! SKIP TO START OF BYTE 33 CALL grbgbyte (msga,kgds(12),kptr(8),32) kptr(8) = kptr(8) + 32 ! ------------------- GO TO 900 ! ****************************************************************** ! ' POLAR STEREO PROCESSING ' ! ! ------------------- BYTE 7-8 NR OF POINTS ALONG X=AXIS 2000 CONTINUE CALL grbgbyte (msga,kgds(2),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 9-10 NR OF POINTS ALONG Y-AXIS CALL grbgbyte (msga,kgds(3),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 11-13 LATITUDE OF ORIGIN CALL grbgbyte (msga,kgds(4),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(4),8388608) /= 0) THEN kgds(4) = - IAND(kgds(4),8388607) END IF ! ------------------- BYTE 14-16 LONGITUDE OF ORIGIN CALL grbgbyte (msga,kgds(5),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(5),8388608) /= 0) THEN kgds(5) = - IAND(kgds(5),8388607) END IF ! ------------------- BYTE 17 RESERVED CALL grbgbyte (msga,kgds(6),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 18-20 LOV ORIENTATION OF THE GRID CALL grbgbyte (msga,kgds(7),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(7),8388608) /= 0) THEN kgds(7) = - IAND(kgds(7),8388607) END IF ! ------------------- BYTE 21-23 DX - THE X DIRECTION INCREMENT CALL grbgbyte (msga,kgds(8),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(8),8388608) /= 0) THEN kgds(8) = - IAND(kgds(8),8388607) END IF ! ------------------- BYTE 24-26 DY - THE Y DIRECTION INCREMENT CALL grbgbyte (msga,kgds(9),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(9),8388608) /= 0) THEN kgds(9) = - IAND(kgds(9),8388607) END IF ! ------------------- BYTE 27 PROJECTION CENTER FLAG CALL grbgbyte (msga,kgds(10),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 28 SCANNING MODE CALL grbgbyte (msga,kgds(11),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 29-32 RESERVED ! SKIP TO START OF BYTE 33 CALL grbgbyte (msga,kgds(12),kptr(8),32) kptr(8) = kptr(8) + 32 ! ! ------------------- GO TO 900 ! ! ****************************************************************** ! ------------------- GRID DESCRIPTION FOR SPHERICAL HARMONIC COEFF. ! ! ------------------- BYTE 7-8 J PENTAGONAL RESOLUTION PARAMETER 3000 CONTINUE CALL grbgbyte (msga,kgds(2),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 9-10 K PENTAGONAL RESOLUTION PARAMETER CALL grbgbyte (msga,kgds(3),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 11-12 M PENTAGONAL RESOLUTION PARAMETER CALL grbgbyte (msga,kgds(4),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 13 REPRESENTATION TYPE CALL grbgbyte (msga,kgds(5),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 14 COEFFICIENT STORAGE MODE CALL grbgbyte (msga,kgds(6),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- EMPTY FIELDS - BYTES 15 - 32 ! SET TO START OF BYTE 33 kptr(8) = kptr(8) + 18 * 8 GO TO 900 ! ****************************************************************** ! PROCESS MERCATOR GRIDS ! ! ------------------- BYTE 7-8 NR OF POINTS ALONG LATITUDE CIRCLE 4000 CONTINUE CALL grbgbyte (msga,kgds(2),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 9-10 NR OF POINTS ALONG LONG MERIDIAN CALL grbgbyte (msga,kgds(3),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 11-13 LATITUE OF ORIGIN CALL grbgbyte (msga,kgds(4),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(4),8388608) /= 0) THEN kgds(4) = - IAND(kgds(4),8388607) END IF ! ------------------- BYTE 14-16 LONGITUDE OF ORIGIN CALL grbgbyte (msga,kgds(5),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(5),8388608) /= 0) THEN kgds(5) = - IAND(kgds(5),8388607) END IF ! ------------------- BYTE 17 RESOLUTION FLAG CALL grbgbyte (msga,kgds(6),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 18-20 LATITUDE OF EXTREME POINT CALL grbgbyte (msga,kgds(7),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(7),8388608) /= 0) THEN kgds(7) = - IAND(kgds(7),8388607) END IF ! ------------------- BYTE 21-23 LONGITUDE OF EXTREME POINT CALL grbgbyte (msga,kgds(8),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(8),8388608) /= 0) THEN kgds(8) = - IAND(kgds(8),8388607) END IF ! ------------------- BYTE 24-26 LATITUDE OF PROJECTION INTERSECTION CALL grbgbyte (msga,kgds(9),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(9),8388608) /= 0) THEN kgds(9) = - IAND(kgds(9),8388607) END IF ! ------------------- BYTE 27 RESERVED CALL grbgbyte (msga,kgds(10),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 28 SCANNING MODE CALL grbgbyte (msga,kgds(11),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 29-31 LONGITUDINAL DIR INCREMENT CALL grbgbyte (msga,kgds(12),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(12),8388608) /= 0) THEN kgds(12) = - IAND(kgds(12),8388607) END IF ! ------------------- BYTE 32-34 LATITUDINAL DIR INCREMENT CALL grbgbyte (msga,kgds(13),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(13),8388608) /= 0) THEN kgds(13) = - IAND(kgds(13),8388607) END IF ! ------------------- BYTE 35-42 RESERVED ! SKIP TO START OF BYTE 43 kptr(8) = kptr(8) + 8 * 8 ! ------------------- GO TO 900 ! ****************************************************************** ! PROCESS LAMBERT CONFORMAL ! ! ------------------- BYTE 7-8 NR OF POINTS ALONG X-AXIS 5000 CONTINUE CALL grbgbyte (msga,kgds(2),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 9-10 NR OF POINTS ALONG Y-AXIS CALL grbgbyte (msga,kgds(3),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- BYTE 11-13 LATITUDE OF ORIGIN CALL grbgbyte (msga,kgds(4),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(4),8388608) /= 0) THEN kgds(4) = - IAND(kgds(4),8388607) END IF ! ------------------- BYTE 14-16 LONGITUDE OF ORIGIN (LOWER LEFT) CALL grbgbyte (msga,kgds(5),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(5),8388608) /= 0) THEN kgds(5) = - IAND(kgds(5),8388607) END IF ! ------------------- BYTE 17 RESOLUTION CALL grbgbyte (msga,kgds(6),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 18-20 LOV -ORIENTATION OF GRID CALL grbgbyte (msga,kgds(7),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(7),8388608) /= 0) THEN kgds(7) = - IAND(kgds(7),8388607) END IF ! ------------------- BYTE 21-23 DX - X-DIR INCREMENT CALL grbgbyte (msga,kgds(8),kptr(8),24) kptr(8) = kptr(8) + 24 ! ------------------- BYTE 24-26 DY - Y-DIR INCREMENT CALL grbgbyte (msga,kgds(9),kptr(8),24) kptr(8) = kptr(8) + 24 ! ------------------- BYTE 27 PROJECTION CENTER FLAG CALL grbgbyte (msga,kgds(10),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 28 SCANNING MODE CALL grbgbyte (msga,kgds(11),kptr(8),8) kptr(8) = kptr(8) + 8 ! ------------------- BYTE 29-31 LATIN1 - 1ST LAT FROM POLE CALL grbgbyte (msga,kgds(12),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(12),8388608) /= 0) THEN kgds(12) = - IAND(kgds(12),8388607) END IF ! ------------------- BYTE 32-34 LATIN2 - 2ND LAT FROM POLE CALL grbgbyte (msga,kgds(13),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(13),8388608) /= 0) THEN kgds(13) = - IAND(kgds(13),8388607) END IF ! ------------------- BYTE 35-37 LATITUDE OF SOUTHERN POLE CALL grbgbyte (msga,kgds(14),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(14),8388608) /= 0) THEN kgds(14) = - IAND(kgds(14),8388607) END IF ! ------------------- BYTE 38-40 LONGITUDE OF SOUTHERN POLE CALL grbgbyte (msga,kgds(15),kptr(8),24) kptr(8) = kptr(8) + 24 IF (IAND(kgds(15),8388608) /= 0) THEN kgds(15) = - IAND(kgds(15),8388607) END IF ! ------------------- BYTE 41-42 RESERVED CALL grbgbyte (msga,kgds(16),kptr(8),16) kptr(8) = kptr(8) + 16 ! ------------------- 900 CONTINUE ! ! MORE CODE FOR GRIDS WITH PL ! IF (kgds(19) == 0.AND.kgds(20) /= 255) THEN isum = 0 kptr(8) = nsave + (kgds(20) - 1) * 8 CALL grbgbytes (msga,kgds(22),kptr(8),16,0,kgds(3)) DO j = 1, kgds(3) isum = isum + kgds(21+j) END DO kgds(21) = isum END IF RETURN END SUBROUTINE fi633 SUBROUTINE fi634(msga,kptr,kpds,kgds,kbms,kret) 1,17 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: FI634 EXTRACT OR GENERATE BIT MAP FOR OUTPUT ! PRGMMR: BILL CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 ! ! ABSTRACT: IF BIT MAP SEC IS AVAILABLE IN GRIB MESSAGE, EXTRACT ! FOR PROGRAM USE, OTHERWISE GENERATE AN APPROPRIATE BIT MAP. ! ! PROGRAM HISTORY LOG: ! 91-09-13 CAVANAUGH ! 91-11-12 CAVANAUGH MODIFIED SIZE OF ECMWF GRIDS 5 - 8. ! ! USAGE: CALL FI634(MSGA,KPTR,KPDS,KGDS,KBMS,KRET) ! INPUT ARGUMENT LIST: ! MSGA - BUFR MESSAGE ! KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS ! (1) - TOTAL LENGTH OF GRIB MESSAGE ! (2) - LENGTH OF INDICATOR (SECTION 0) ! (3) - LENGTH OF PDS (SECTION 1) ! (4) - LENGTH OF GDS (SECTION 2) ! (5) - LENGTH OF BMS (SECTION 3) ! (6) - LENGTH OF BDS (SECTION 4) ! (7) - VALUE OF CURRENT BYTE ! (8) - BIT POINTER ! (9) - GRIB START BIT NR ! (10) - GRIB/GRID ELEMENT COUNT ! (11) - NR UNUSED BITS AT END OF SECTION 3 ! (12) - BIT MAP FLAG ! (13) - NR UNUSED BITS AT END OF SECTION 2 ! (14) - BDS FLAGS ! (15) - NR UNUSED BITS AT END OF SECTION 4 ! KPDS - ARRAY CONTAINING PDS ELEMENTS. ! (1) - ID OF CENTER ! (2) - MODEL IDENTIFICATION ! (3) - GRID IDENTIFICATION ! (4) - GDS/BMS FLAG ! (5) - INDICATOR OF PARAMETER ! (6) - TYPE OF LEVEL ! (7) - HEIGHT/PRESSURE , ETC OF LEVEL ! (8) - YEAR OF CENTURY ! (9) - MONTH OF YEAR ! (10) - DAY OF MONTH ! (11) - HOUR OF DAY ! (12) - MINUTE OF HOUR ! (13) - INDICATOR OF FORECAST TIME UNIT ! (14) - TIME RANGE 1 ! (15) - TIME RANGE 2 ! (16) - TIME RANGE FLAG ! (17) - NUMBER INCLUDED IN AVERAGE ! ! OUTPUT ARGUMENT LIST: ! KBMS - BITMAP DESCRIBING LOCATION OF OUTPUT ELEMENTS. ! KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS ! SEE INPUT LIST ! KRET - ERROR RETURN ! ! REMARKS: ! KRET = 0 - NO ERROR ! = 5 - GRID NOT AVAIL FOR CENTER INDICATED ! =10 - INCORRECT CENTER INDICATOR ! =12 - BYTES 5-6 ARE NOT ZERO IN BMS, PREDEFINED BIT MAP ! NOT PROVIDED BY THIS CENTER ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 77 ! MACHINE: NAS9000 ! !$$$ ! ! INCOMING MESSAGE HOLDER CHARACTER :: msga(*) ! ! BIT MAP LOGICAL*1 :: kbms(*) ! ! ARRAY OF POINTERS AND COUNTERS INTEGER :: kptr(*) ! ARRAY OF POINTERS AND COUNTERS INTEGER :: kpds(*) INTEGER :: kgds(*) ! INTEGER :: kret INTEGER :: mask(8) ! ----------------------GRID 21 AND GRID 22 ARE THE SAME LOGICAL*1 grd21( 1369) ! ----------------------GRID 23 AND GRID 24 ARE THE SAME LOGICAL*1 grd23( 1369) LOGICAL*1 grd25( 1368) LOGICAL*1 grd26( 1368) ! ----------------------GRID 27 AND GRID 28 ARE THE SAME ! ----------------------GRID 29 AND GRID 30 ARE THE SAME ! ----------------------GRID 33 AND GRID 34 ARE THE SAME LOGICAL*1 grd50( 1188) ! -----------------------GRID 61 AND GRID 62 ARE THE SAME LOGICAL*1 grd61( 4186) ! -----------------------GRID 63 AND GRID 64 ARE THE SAME LOGICAL*1 grd63( 4186) ! LOGICAL*1 GRD70(16380)/16380*.TRUE./ ! ------------------------------------------------------------- SAVE DATA grd21 /1333*.true.,36*.false./ DATA grd23 /.true.,36*.false.,1332*.true./ DATA grd25 /1297*.true.,71*.false./ DATA grd26 /.true.,71*.false.,1296*.true./ DATA grd50/ & ! LINE 1-4 7*.false.,22*.true.,14*.false.,22*.true., & 14*.false.,22*.true.,14*.false.,22*.true.,7*.false., & ! LINE 5-8 6*.false.,24*.true.,12*.false.,24*.true., & 12*.false.,24*.true.,12*.false.,24*.true.,6*.false., & ! LINE 9-12 5*.false.,26*.true.,10*.false.,26*.true., & 10*.false.,26*.true.,10*.false.,26*.true.,5*.false., & ! LINE 13-16 4*.false.,28*.true., 8*.false.,28*.true., & 8*.false.,28*.true., 8*.false.,28*.true.,4*.false., & ! LINE 17-20 3*.false.,30*.true., 6*.false.,30*.true., & 6*.false.,30*.true., 6*.false.,30*.true.,3*.false., & ! LINE 21-24 2*.false.,32*.true., 4*.false.,32*.true., & 4*.false.,32*.true., 4*.false.,32*.true.,2*.false., & ! LINE 25-28 .false.,34*.true., 2*.false.,34*.true., & 2*.false.,34*.true., 2*.false.,34*.true., .false., & ! LINE 29-33 180*.true./ DATA grd61 /4096*.true.,90*.false./ DATA grd63 /.true.,90*.false.,4095*.true./ DATA mask /128,64,32,16,8,4,2,1/ ! ! PRINT *,'FI634' IF (IAND(kpds(4),64) == 64) THEN ! ! SET UP BIT POINTER ! SECTION 0 SECTION 1 SECTION 2 kptr(8) = kptr(9) + (kptr(2)*8) + (kptr(3)*8) + (kptr(4)*8) + 24 ! ! BYTE 4 NUMBER OF UNUSED BITS AT END OF SECTION 3 ! CALL grbgbyte (msga,kptr(11),kptr(8),8) kptr(8) = kptr(8) + 8 ! ! BYTE 5,6 TABLE REFERENCE IF 0, BIT MAP FOLLOWS ! CALL grbgbyte (msga,kptr(12),kptr(8),16) kptr(8) = kptr(8) + 16 ! IF TABLE REFERENCE = 0, EXTRACT BIT MAP IF (kptr(12) == 0) THEN ! CALCULATE NR OF BITS IN BIT MAP ibits = (kptr(5) - 6) * 8 - kptr(11) kptr(10) = ibits IF (kpds(3) == 21.OR.kpds(3) == 22.OR.kpds(3) == 25 & .OR.kpds(3) == 61.OR.kpds(3) == 62) THEN ! NORTHERN HEMISPHERE 21, 22, 25, 61, 62 DO i = 1, ibits CALL grbgbyte (msga,ichk,kptr(8),1) kptr(8) = kptr(8) + 1 IF (ichk /= 0) THEN kbms(i) = .true. ELSE kbms(i) = .false. END IF END DO IF (kpds(3) == 25) THEN kadd = 71 ELSE IF (kpds(3) == 61.OR.kpds(3) == 62) THEN kadd = 90 ELSE kadd = 36 END IF DO i = 1, kadd kbms(i+ibits) = .false. END DO kptr(10) = kptr(10) + kadd RETURN ELSE IF (kpds(3) == 23.OR.kpds(3) == 24.OR.kpds(3) == 26 & .OR.kpds(3) == 63.OR.kpds(3) == 64) THEN ! SOUTHERN HEMISPHERE 23, 24, 26, 63, 64 IF (kpds(3) == 26) THEN kadd = 72 ELSE IF (kpds(3) == 63.OR.kpds(3) == 64) THEN kadd = 91 ELSE kadd = 37 END IF DO i = 1, kadd kbms(i+ibits) = .false. END DO DO i = 1, ibits CALL grbgbyte (msga,ichk,kptr(8),1) kptr(8) = kptr(8) + 1 IF (ichk /= 0) THEN kbms(i) = .true. ELSE kbms(i) = .false. END IF END DO kptr(10) = kptr(10) + kadd - 1 RETURN ELSE IF (kpds(3) == 50) THEN kpad = 7 kin = 22 kbits = 0 DO i = 1, 7 DO j = 1, 4 DO k = 1, kpad kbits = kbits + 1 kbms(kbits) = .false. END DO DO k = 1, kin CALL grbgbyte (msga,ichk,kptr(8),1) kptr(8) = kptr(8) + 1 kbits = kbits + 1 IF (ichk /= 0) THEN kbms(kbits) = .true. ELSE kbms(kbits) = .false. END IF END DO DO k = 1, kpad kbits = kbits + 1 kbms(kbits) = .false. END DO END DO kin = kin + 2 kpad = kpad - 1 END DO DO ii = 1, 5 DO j = 1, kin CALL grbgbyte (msga,ichk,kptr(8),1) kptr(8) = kptr(8) + 1 kbits = kbits + 1 IF (ichk /= 0) THEN kbms(kbits) = .true. ELSE kbms(kbits) = .false. END IF END DO END DO ELSE ! EXTRACT BIT MAP FROM BMS FOR OTHER GRIDS DO i = 1, ibits CALL grbgbyte (msga,ichk,kptr(8),1) kptr(8) = kptr(8) + 1 IF (ichk /= 0) THEN kbms(i) = .true. ELSE kbms(i) = .false. END IF END DO END IF RETURN ELSE PRINT *,'FI634-NO PREDEFINED BIT MAP PROVIDED BY THIS CENTER' kret = 12 RETURN END IF ! END IF kret = 0 ! ------------------------------------------------------- ! PROCESS NON-STANDARD GRID ! ------------------------------------------------------- IF (kpds(3) == 255) THEN PRINT *,'NON STANDARD GRID, CENTER = ',kpds(1) j = kgds(2) * kgds(3) kptr(10) = j DO i = 1, j kbms(i) = .true. END DO RETURN END IF ! ------------------------------------------------------- ! CHECK INTERNATIONAL SET ! ------------------------------------------------------- IF (kpds(3) == 21.OR.kpds(3) == 22) THEN ! ----- INT'L GRIDS 21, 22 - MAP SIZE 1369 j = 1369 kptr(10) = j CALL fi637(*820,j,kpds,kgds,kret) DO i = 1, 1369 kbms(i) = grd21(i) END DO RETURN ELSE IF (kpds(3) == 23.OR.kpds(3) == 24) THEN ! ----- INT'L GRIDS 23, 24 - MAP SIZE 1369 j = 1369 kptr(10) = j CALL fi637(*820,j,kpds,kgds,kret) DO i = 1, 1369 kbms(i) = grd23(i) END DO RETURN ELSE IF (kpds(3) == 25) THEN ! ----- INT'L GRID 25 - MAP SIZE 1368 j = 1368 kptr(10) = j CALL fi637(*820,j,kpds,kgds,kret) DO i = 1, 1368 kbms(i) = grd25(i) END DO RETURN ELSE IF (kpds(3) == 26) THEN ! ----- INT'L GRID 26 - MAP SIZE 1368 j = 1368 kptr(10) = j CALL fi637(*820,j,kpds,kgds,kret) DO i = 1, 1368 kbms(i) = grd26(i) END DO RETURN ELSE IF (kpds(3) >= 37.AND.kpds(3) <= 44) THEN ! ----- INT'L GRID 37-44 - MAP SIZE 3447 j = 3447 GO TO 800 ELSE IF (kpds(1) == 7.AND.kpds(3) == 50) THEN ! ----- INT'L GRIDS 50 - MAP SIZE 964 j = 1188 kptr(10) = j CALL fi637(*890,j,kpds,kgds,kret) DO i = 1, j kbms(i) = grd50(i) END DO RETURN ELSE IF (kpds(3) == 61.OR.kpds(3) == 62) THEN ! ----- INT'L GRIDS 61, 62 - MAP SIZE 4186 j = 4186 kptr(10) = j CALL fi637(*820,j,kpds,kgds,kret) DO i = 1, 4186 kbms(i) = grd61(i) END DO RETURN ELSE IF (kpds(3) == 63.OR.kpds(3) == 64) THEN ! ----- INT'L GRIDS 63, 64 - MAP SIZE 4186 j = 4186 kptr(10) = j CALL fi637(*820,j,kpds,kgds,kret) DO i = 1, 4186 kbms(i) = grd63(i) END DO RETURN END IF ! ------------------------------------------------------- ! CHECK UNITED STATES SET ! ------------------------------------------------------- IF (kpds(1) == 7) THEN IF (kpds(3) < 100) THEN IF (kpds(3) == 1) THEN ! ----- U.S. GRID 1 - MAP SIZE 1679 j = 1679 GO TO 800 END IF IF (kpds(3) == 2) THEN ! ----- U.S. GRID 2 - MAP SIZE 10512 j = 10512 GO TO 800 ELSE IF (kpds(3) == 3) THEN ! ----- U.S. GRID 3 - MAP SIZE 65160 j = 65160 GO TO 800 ELSE IF (kpds(3) == 5) THEN ! ----- U.S. GRID 5 - MAP SIZE 3021 j = 3021 GO TO 800 ELSE IF (kpds(3) == 6) THEN ! ----- U.S. GRID 6 - MAP SIZE 2385 j = 2385 GO TO 800 ELSE IF (kpds(3) == 27.OR.kpds(3) == 28) THEN ! ----- U.S. GRIDS 27, 28 - MAP SIZE 4225 j = 4225 GO TO 800 ELSE IF (kpds(3) == 29.OR.kpds(3) == 30) THEN ! ----- U.S. GRIDS 29,30 - MAP SIZE 5365 j = 5365 GO TO 800 ELSE IF (kpds(3) == 33.OR.kpds(3) == 34) THEN ! ----- U.S GRID 33, 34 - MAP SIZE 8326 j = 8326 GO TO 800 ELSE IF (kpds(3) >= 37.AND.kpds(3) <= 44) THEN ! ----- U.S. GRID 37-44 - MAP SIZE 3447 j = 3447 GO TO 800 ELSE IF (kpds(3) == 45) THEN ! ----- U.S. GRID 45 - MAP SIZE 41760 j = 41760 GO TO 800 ELSE IF (kpds(3) == 55.OR.kpds(3) == 56) THEN ! ----- U.S GRID 55, 56 - MAP SIZE 6177 j = 6177 GO TO 800 ELSE IF (kpds(3) >= 67.AND.kpds(3) <= 71) THEN ! ----- U.S GRID 67-71 - MAP SIZE 13689 j = 13689 GO TO 800 ELSE IF (kpds(3) == 72) THEN ! ----- U.S GRID 72 - MAP SIZE 406 j = 406 GO TO 800 ELSE IF (kpds(3) == 73) THEN ! ----- U.S GRID 73 - MAP SIZE 13056 j = 13056 GO TO 800 ELSE IF (kpds(3) == 74) THEN ! ----- U.S GRID 74 - MAP SIZE 10800 j = 10800 GO TO 800 ELSE IF (kpds(3) >= 75.AND.kpds(3) <= 77) THEN ! ----- U.S GRID 75-77 - MAP SIZE 12321 j = 12321 GO TO 800 ELSE IF (kpds(3) == 85.OR.kpds(3) == 86) THEN ! ----- U.S GRID 85,86 - MAP SIZE 32400 j = 32400 GO TO 800 ELSE IF (kpds(3) == 87) THEN ! ----- U.S GRID 87 - MAP SIZE 5022 j = 5022 GO TO 800 ELSE IF (kpds(3) == 90) THEN ! ----- U.S GRID 90 - MAP SIZE 12902 j = 12902 GO TO 800 ELSE IF (kpds(3) == 91) THEN ! ----- U.S GRID 91 - MAP SIZE 25803 j = 25803 GO TO 800 ELSE IF (kpds(3) == 92) THEN ! ----- U.S GRID 92 - MAP SIZE 24162 j = 24162 GO TO 800 ELSE IF (kpds(3) == 93) THEN ! ----- U.S GRID 93 - MAP SIZE 48323 j = 48323 GO TO 800 ELSE IF (kpds(3) == 94) THEN ! ----- U.S GRID 94 - MAP SIZE XXXXX j = 48916 GO TO 800 ELSE IF (kpds(3) == 95) THEN ! ----- U.S GRID 95 - MAP SIZE XXXXX j = 97831 GO TO 800 END IF ELSE IF (kpds(3) >= 100.AND.kpds(3) < 200) THEN IF (kpds(3) == 100) THEN ! ----- U.S. GRID 100 - MAP SIZE 6889 j = 6889 GO TO 800 ELSE IF (kpds(3) == 101) THEN ! ----- U.S. GRID 101 - MAP SIZE 10283 j = 10283 GO TO 800 ELSE IF (kpds(3) == 103) THEN ! ----- U.S. GRID 103 - MAP SIZE 3640 j = 3640 GO TO 800 ELSE IF (kpds(3) == 104) THEN ! ----- U.S. GRID 104 - MAP SIZE 16170 j = 16170 GO TO 800 ELSE IF (kpds(3) == 105) THEN ! ----- U.S. GRID 105 - MAP SIZE 6889 j = 6889 GO TO 800 ELSE IF (kpds(3) == 106) THEN ! ----- U.S. GRID 106 - MAP SIZE 19305 j = 19305 GO TO 800 ELSE IF (kpds(3) == 107) THEN ! ----- U.S. GRID 107 - MAP SIZE 11040 j = 11040 GO TO 800 ELSE IF (IAND(kpds(4),128) == 128) THEN ! ----- U.S. NON-STANDARD GRID GO TO 895 END IF ELSE IF (kpds(3) >= 200) THEN IF (kpds(3) == 201) THEN j = 4225 GO TO 800 ELSE IF (kpds(3) == 202) THEN j = 2795 GO TO 800 ELSE IF (kpds(3) == 203.OR.kpds(3) == 205) THEN j = 1755 GO TO 800 ELSE IF (kpds(3) == 204) THEN j = 6324 GO TO 800 ELSE IF (kpds(3) == 206) THEN j = 2091 GO TO 800 ELSE IF (kpds(3) == 207) THEN j = 1715 GO TO 800 ELSE IF (kpds(3) == 208) THEN j = 783 GO TO 800 ELSE IF (kpds(3) == 209) THEN j = 8181 GO TO 800 ELSE IF (kpds(3) == 210) THEN j = 625 GO TO 800 ELSE IF (kpds(3) == 211) THEN j = 6045 GO TO 800 ELSE IF (kpds(3) == 212) THEN j = 23865 GO TO 800 ELSE IF (kpds(3) == 213) THEN j = 10965 GO TO 800 ELSE IF (kpds(3) == 214) THEN j = 6693 GO TO 800 ELSE IF (IAND(kpds(4),128) == 128) THEN GO TO 895 END IF kret = 5 RETURN END IF END IF ! ------------------------------------------------------- ! CHECK JAPAN METEOROLOGICAL AGENCY SET ! ------------------------------------------------------- IF (kpds(1) == 34) THEN IF (IAND(kpds(4),128) == 128) THEN PRINT *,'JMA MAP IS NOT PREDEFINED, THE GDS WILL' PRINT *,'BE USED TO UNPACK THE DATA, MAP = ',kpds(3) GO TO 900 END IF END IF ! ------------------------------------------------------- ! CHECK CANADIAN SET ! ------------------------------------------------------- IF (kpds(1) == 54) THEN IF (IAND(kpds(4),128) == 128) THEN PRINT *,'CANADIAN MAP IS NOT PREDEFINED, THE GDS WILL' PRINT *,'BE USED TO UNPACK THE DATA, MAP = ',kpds(3) GO TO 900 END IF END IF ! ------------------------------------------------------- ! CHECK FNOC SET ! ------------------------------------------------------- IF (kpds(1) == 58) THEN IF (kpds(3) == 220.OR.kpds(3) == 221) THEN ! FNOC GRID 220, 221 - MAPSIZE 3969 (63 * 63) j = 3969 kptr(10) = j DO i = 1, j kbms(i) = .true. END DO RETURN END IF IF (kpds(3) == 223) THEN ! FNOC GRID 223 - MAPSIZE 10512 (73 * 144) j = 10512 kptr(10) = j DO i = 1, j kbms(i) = .true. END DO RETURN END IF IF (IAND(kpds(4),128) == 128) THEN PRINT *,'FNOC MAP IS NOT PREDEFINED, THE GDS WILL' PRINT *,'BE USED TO UNPACK THE DATA, MAP = ',kpds(3) GO TO 900 END IF END IF ! ------------------------------------------------------- ! CHECK UKMET SET ! ------------------------------------------------------- IF (kpds(1) == 74) THEN IF (IAND(kpds(4),128) == 128) THEN GO TO 820 END IF END IF ! ------------------------------------------------------- ! CHECK ECMWF SET ! ------------------------------------------------------- IF (kpds(1) == 98) THEN IF (kpds(3) >= 1.AND.kpds(3) <= 12) THEN IF (kpds(3) >= 5.AND.kpds(3) <= 8) THEN j = 1073 ELSE j = 1369 END IF kptr(10) = j CALL fi637(*810,j,kpds,kgds,kret) DO i = 1, j kbms(i) = .true. END DO RETURN ELSE IF (kpds(3) >= 13.AND.kpds(3) <= 16) THEN j = 361 kptr(10) = j CALL fi637(*810,j,kpds,kgds,kret) DO i = 1, j kbms(i) = .true. END DO RETURN ELSE IF (IAND(kpds(4),128) == 128) THEN GO TO 810 ELSE kret = 5 RETURN END IF ELSE PRINT *,'CENTER',kpds(1),' IS NOT DEFINED' IF (IAND(kpds(4),128) == 128) THEN PRINT *,'GDS WILL BE USED TO UNPACK THE DATA', & ' MAP = ',kpds(3) GO TO 900 ELSE kret = 10 RETURN END IF END IF ! ======================================= ! 800 CONTINUE kptr(10) = j CALL fi637 (*801,j,kpds,kgds,kret) DO i = 1, j kbms(i) = .true. END DO RETURN 801 CONTINUE ! ! ----- THE MAP HAS A GDS, BYTE 7 OF THE (PDS) THE GRID IDENTIFICATION ! ----- IS NOT 255, THE SIZE OF THE GRID IS NOT THE SAME AS THE ! ----- PREDEFINED SIZES OF THE U.S. GRIDS, OR KNOWN GRIDS OF THE ! ----- OF THE OTHER CENTERS. THE GRID CAN BE UNKNOWN, OR FROM AN ! ----- UNKNOWN CENTER, WE WILL USE THE INFORMATION IN THE GDS TO MAKE ! ----- A BIT MAP. ! 810 CONTINUE PRINT *,'ECMWF PREDEFINED MAP SIZE DOES NOT MATCH, I WILL USE' GO TO 895 ! 820 CONTINUE PRINT *,'U.K. MET PREDEFINED MAP SIZE DOES NOT MATCH, I WILL USE' GO TO 895 ! 890 CONTINUE PRINT *,'PREDEFINED MAP SIZE DOES NOT MATCH, I WILL USE' 895 CONTINUE PRINT *,'THE GDS TO UNPACK THE DATA, MAP TYPE = ',kpds(3) ! 900 CONTINUE j = kgds(2) * kgds(3) ! AFOS AFOS AFOS SPECIAL CASE ! INVOLVES NEXT SINGLE STATEMENT ONLY IF (kpds(3) == 211) kret = 0 kptr(10) = j DO i = 1, j kbms(i) = .true. END DO ! PRINT *,'EXIT FI634' RETURN END SUBROUTINE fi634 SUBROUTINE fi635(msga,kptr,kpds,kgds,kbms,DATA,kret) 1,37 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: FI635 EXTRACT GRIB DATA ELEMENTS FROM BDS ! PRGMMR: BILL CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 ! ! ABSTRACT: EXTRACT GRIB DATA FROM BINARY DATA SECTION AND PLACE ! INTO OUTPUT ARRAY IN PROPER POSITION. ! ! PROGRAM HISTORY LOG: ! 91-09-13 CAVANAUGH ! 94-04-01 CAVANAUGH MODIFIED CODE TO INCLUDE DECIMAL SCALING WHEN ! CALCULATING THE VALUE OF DATA POINTS SPECIFIED ! AS BEING EQUAL TO THE REFERENCE VALUE ! 94-11-10 FARLEY INCREASED MXSIZE FROM 72960 TO 260000 ! FOR .5 DEGREE SST ANALYSIS FIELDS ! ! USAGE: CALL FI635(MSGA,KPTR,KPDS,KGDS,KBMS,DATA,KRET) ! INPUT ARGUMENT LIST: ! MSGA - ARRAY CONTAINING GRIB MESSAGE ! KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS ! (1) - TOTAL LENGTH OF GRIB MESSAGE ! (2) - LENGTH OF INDICATOR (SECTION 0) ! (3) - LENGTH OF PDS (SECTION 1) ! (4) - LENGTH OF GDS (SECTION 2) ! (5) - LENGTH OF BMS (SECTION 3) ! (6) - LENGTH OF BDS (SECTION 4) ! (7) - VALUE OF CURRENT BYTE ! (8) - BIT POINTER ! (9) - GRIB START BIT NR ! (10) - GRIB/GRID ELEMENT COUNT ! (11) - NR UNUSED BITS AT END OF SECTION 3 ! (12) - BIT MAP FLAG ! (13) - NR UNUSED BITS AT END OF SECTION 2 ! (14) - BDS FLAGS ! (15) - NR UNUSED BITS AT END OF SECTION 4 ! KPDS - ARRAY CONTAINING PDS ELEMENTS. ! SEE INITIAL ROUTINE ! KBMS - BITMAP DESCRIBING LOCATION OF OUTPUT ELEMENTS. ! ! OUTPUT ARGUMENT LIST: ! KBDS - INFORMATION EXTRACTED FROM BINARY DATA SECTION ! KBDS(1) - N1 ! KBDS(2) - N2 ! KBDS(3) - P1 ! KBDS(4) - P2 ! KBDS(5) - BIT POINTER TO 2ND ORDER WIDTHS ! KBDS(6) - " " " " " BIT MAPS ! KBDS(7) - " " " FIRST ORDER VALUES ! KBDS(8) - " " " SECOND ORDER VALUES ! KBDS(9) - " " START OF BDS ! KBDS(10) - " " MAIN BIT MAP ! KBDS(11) - BINARY SCALING ! KBDS(12) - DECIMAL SCALING ! KBDS(13) - BIT WIDTH OF FIRST ORDER VALUES ! KBDS(14) - BIT MAP FLAG ! 0 = NO SECOND ORDER BIT MAP ! 1 = SECOND ORDER BIT MAP PRESENT ! KBDS(15) - SECOND ORDER BIT WIDTH ! KBDS(16) - CONSTANT / DIFFERENT WIDTHS ! 0 = CONSTANT WIDTHS ! 1 = DIFFERENT WIDTHS ! KBDS(17) - SINGLE DATUM / MATRIX ! 0 = SINGLE DATUM AT EACH GRID POINT ! 1 = MATRIX OF VALUES AT EACH GRID POINT ! (18-20)- UNUSED ! ! DATA - REAL*4 ARRAY OF GRIDDED ELEMENTS IN GRIB MESSAGE. ! KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS ! SEE INPUT LIST ! KRET - ERROR RETURN ! ! REMARKS: ! ERROR RETURN ! 3 = UNPACKED FIELD IS LARGER THAN 65160 ! 6 = DOES NOT MATCH NR OF ENTRIES FOR THIS GRIB/GRID ! 7 = NUMBER OF BITS IN FILL TOO LARGE ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 77 ! MACHINE: NAS9000 ! !$$$ ! ************************************************************* ! ON A PC THIS CAN BE CHANGED TO A SMALLER SIZE TO BETTER FIT ! THE DOS MEMORY LIMIT OF 640K BYTES. YOU COULD DO THIS ! FOR MICROSOFT 5.0. A PC 32 BIT FORTRAN COMPILER ! WOULD NOT NEED THIS CHANGE. IF NONE OF YOUR GRIB RECORDS ! IS LARGER THAN 20000, SET MXSIZE TO 20000. ! ************************************************************* PARAMETER (mxsize=260000) CHARACTER :: msga(*) CHARACTER (LEN=1) :: kk(8) CHARACTER (LEN=1) :: ckref(8) ! LOGICAL*1 kbms(*) ! INTEGER :: kpds(*) INTEGER :: kgds(*) INTEGER :: kbds(20) INTEGER :: kptr(*) INTEGER :: nrbits INTEGER :: kref INTEGER :: kkk INTEGER :: ksave(mxsize) INTEGER :: kscale ! REAL :: DATA(*) REAL :: refnce REAL :: scale REAL :: realkk ! EQUIVALENCE (ckref(1),kref,refnce) EQUIVALENCE (kk(1),kkk,realkk) ! ! ! CHANGED HEX VALUES TO DECIMAL TO MAKE CODE MORE PORTABLE ! ! ************************************************************* SAVE ! PRINT *,'ENTER FI635' ! SET UP BIT POINTER kptr(8) = kptr(9) + (kptr(2)*8) + (kptr(3)*8) + (kptr(4)*8) & + (kptr(5)*8) + 24 ! ------------- EXTRACT FLAGS ! BYTE 4 CALL grbgbyte(msga,kptr(14),kptr(8),4) kptr(8) = kptr(8) + 4 ! --------- NR OF UNUSED BITS IN SECTION 4 CALL grbgbyte(msga,kptr(15),kptr(8),4) kptr(8) = kptr(8) + 4 kend = kptr(9) + (kptr(2)*8) + (kptr(3)*8) + (kptr(4)*8) & + (kptr(5)*8) + kptr(6) * 8 - kptr(15) ! ------------- GET SCALE FACTOR ! BYTES 5,6 ! CHECK SIGN CALL grbgbyte (msga,ksign,kptr(8),1) kptr(8) = kptr(8) + 1 ! GET ABSOLUTE SCALE VALUE CALL grbgbyte (msga,kscale,kptr(8),15) kptr(8) = kptr(8) + 15 IF (ksign > 0) THEN kscale = - kscale END IF scale = 2.0**kscale ! ------------ GET REFERENCE VALUE ! BYTES 7,10 CALL grbgbyte (msga,kref,kptr(8),32) kptr(8) = kptr(8) + 32 ! ! THE NEXT CODE WILL CONVERT THE IBM370 FLOATING POINT ! TO THE FLOATING POINT USED ON YOUR COMPUTER. ! ! 1ST TEST TO SEE IN ON 32 OR 64 BIT WORD MACHINE ! LW = 4 OR 8; IF 8 MAY BE A CRAY ! CALL w3fi01(lw) IF (lw == 4) THEN CALL grbgbyte (ckref,jsgn,0,1) CALL grbgbyte (ckref,jexp,1,7) CALL grbgbyte (ckref,ifr,8,24) ELSE CALL grbgbyte (ckref,jsgn,32,1) CALL grbgbyte (ckref,jexp,33,7) CALL grbgbyte (ckref,ifr,40,24) END IF ! PRINT *,109,JSGN,JEXP,IFR ! 109 FORMAT (' JSGN,JEXP,IFR = ',3(1X,Z8)) IF (ifr == 0) THEN refnce = 0.0 ELSE IF (jexp == 0.AND.ifr == 0) THEN refnce = 0.0 ELSE refnce = FLOAT(ifr) * 16.0 ** (jexp - 64 - 6) IF (jsgn /= 0) refnce = - refnce END IF ! PRINT *,'SCALE ',SCALE,' REF VAL ',KREF,REFNCE ! ------------- NUMBER OF BITS SPECIFIED FOR EACH ENTRY ! BYTE 11 CALL grbgbyte (msga,kbits,kptr(8),8) kptr(8) = kptr(8) + 8 kbds(4) = kbits ! KBDS(13) = KBITS ibyt12 = kptr(8) ! ------------------ IF THERE ARE NO EXTENDED FLAGS PRESENT ! THIS IS WHERE DATA BEGINS AND AND THE PROCESSING ! INCLUDED IN THE FOLLOWING IF...END IF ! WILL BE SKIPPED ! PRINT *,'BASIC FLAGS =',KPTR(14) ,IAND(KPTR(14),1) IF (IAND(kptr(14),1) == 0) THEN ! PRINT *,'NO EXTENDED FLAGS' ELSE ! BYTES 12,13 CALL grbgbyte (msga,koctet,kptr(8),16) kptr(8) = kptr(8) + 16 ! --------------------------- EXTENDED FLAGS ! BYTE 14 CALL grbgbyte (msga,kxflag,kptr(8),8) ! PRINT *,'HAVE EXTENDED FLAGS',KXFLAG kptr(8) = kptr(8) + 8 IF (IAND(kxflag,16) == 0) THEN ! SECOND ORDER VALUES CONSTANT WIDTHS kbds(16) = 0 ELSE ! SECOND ORDER VALUES DIFFERENT WIDTHS kbds(16) = 1 END IF IF (IAND (kxflag,32) == 0) THEN ! NO SECONDARY BIT MAP kbds(14) = 0 ELSE ! HAVE SECONDARY BIT MAP kbds(14) = 1 END IF IF (IAND (kxflag,64) == 0) THEN ! SINGLE DATUM AT GRID POINT kbds(17) = 0 ELSE ! MATRIX OF VALUES AT GRID POINT kbds(17) = 1 END IF ! ---------------------- NR - FIRST DIMENSION (ROWS) OF EACH MATRIX ! BYTES 15,16 CALL grbgbyte (msga,nr,kptr(8),16) kptr(8) = kptr(8) + 16 ! ---------------------- NC - SECOND DIMENSION (COLS) OF EACH MATRIX ! BYTES 17,18 CALL grbgbyte (msga,nc,kptr(8),16) kptr(8) = kptr(8) + 16 ! ---------------------- NRV - FIRST DIM COORD VALS ! BYTE 19 CALL grbgbyte (msga,nrv,kptr(8),8) kptr(8) = kptr(8) + 8 ! ---------------------- NC1 - NR COEFF'S OR VALUES ! BYTE 20 CALL grbgbyte (msga,nc1,kptr(8),8) kptr(8) = kptr(8) + 8 ! ---------------------- NCV - SECOND DIM COORD OR VALUE ! BYTE 21 CALL grbgbyte (msga,ncv,kptr(8),8) kptr(8) = kptr(8) + 8 ! ---------------------- NC2 - NR COEFF'S OR VALS ! BYTE 22 CALL grbgbyte (msga,nc2,kptr(8),8) kptr(8) = kptr(8) + 8 ! ---------------------- KPHYS1 - FIRST DIM PHYSICAL SIGNIF ! BYTE 23 CALL grbgbyte (msga,kphys1,kptr(8),8) kptr(8) = kptr(8) + 8 ! ---------------------- KPHYS2 - SECOND DIM PHYSICAL SIGNIF ! BYTE 24 CALL grbgbyte (msga,kphys2,kptr(8),8) kptr(8) = kptr(8) + 8 ! BYTES 25-N END IF IF (kbits == 0) THEN ! HAVE NO BDS ENTRIES, ALL ENTRIES = REFNCE scal10 = 10.0 ** kpds(22) scal10 = 1.0 / scal10 refn10 = refnce * scal10 kentry = kptr(10) DO i = 1, kentry DATA(i) = 0.0 IF (kbms(i)) THEN DATA(i) = refn10 END IF END DO GO TO 900 END IF ! PRINT *,'KEND ',KEND,' KPTR(8) ',KPTR(8),'KBITS ',KBITS knr = (kend - kptr(8)) / kbits ! PRINT *,'NUMBER OF ENTRIES IN DATA ARRAY',KNR ! -------------------- ! CYCLE THRU BDS UNTIL HAVE USED ALL (SPECIFIED NUMBER) ! ENTRIES. ! ------------- UNUSED BITS IN DATA AREA ! NUMBER OF BYTES IN DATA AREA nrbyte = kptr(6) - 11 ! ------------- TOTAL NR OF USABLE BITS nrbits = nrbyte * 8 - kptr(15) ! ------------- TOTAL NR OF ENTRIES kentry = nrbits / kbits ! MAX SIZE CHECK IF (kentry > mxsize) THEN kret = 3 RETURN END IF ! ! IF (IAND(KPTR(14),2).EQ.0) THEN ! PRINT *,'SOURCE VALUES IN FLOATING POINT' ! ELSE ! PRINT *,'SOURCE VALUES IN INTEGER' ! END IF ! IF (IAND(kptr(14),8) == 0) THEN ! PRINT *,'PROCESSING GRID POINT DATA' IF (IAND(kptr(14),4) == 0) THEN ! PRINT *,' WITH SIMPLE PACKING' IF (IAND(kptr(14),1) == 0) THEN ! PRINT *,' WITH NO ADDITIONAL FLAGS' GO TO 4000 ELSE IF (IAND(kptr(14),1) /= 0) THEN PRINT *,' WITH ADDITIONAL FLAGS',kxflag IF (kbds(17) == 0) THEN PRINT *,' SINGLE DATUM EACH GRID PT' IF (kbds(14) == 0) THEN PRINT *,' NO SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (kbds(14) /= 0) THEN PRINT *,' SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF END IF ELSE IF (kbds(17) /= 0) THEN PRINT *,' MATRIX OF VALS EACH PT' IF (kbds(14) == 0) THEN PRINT *,' NO SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (kbds(14) /= 0) THEN PRINT *,' SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF END IF END IF END IF ELSE IF (IAND(kptr(14),4) /= 0) THEN PRINT *,' WITH COMPLEX/SECOND ORDER PACKING' IF (IAND(kptr(14),1) == 0) THEN PRINT *,' WITH NO ADDITIONAL FLAGS' ELSE IF (IAND(kptr(14),1) /= 0) THEN PRINT *,' WITH ADDITIONAL FLAGS' IF (kbds(17) == 0) THEN PRINT *,' SINGLE DATUM AT EACH PT' IF (kbds(14) == 0) THEN PRINT *,' NO SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF ! ROW BY ROW - COL BY COL CALL fi636 (DATA,msga,kbms, & refnce,kptr,kpds,kgds) GO TO 900 ELSE IF (kbds(14) /= 0) THEN PRINT *,' SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF CALL fi636 (DATA,msga,kbms, & refnce,kptr,kpds,kgds) GO TO 900 END IF ELSE IF (kbds(17) /= 0) THEN PRINT *,' MATRIX OF VALS EACH PT' IF (kbds(14) == 0) THEN PRINT *,' NO SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (kbds(14) /= 0) THEN PRINT *,' SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF END IF END IF END IF END IF ELSE IF (IAND(kptr(14),8) /= 0) THEN PRINT *,'PROCESSING SPHERICAL HARMONIC COEFFICIENTS' IF (IAND(kptr(14),4) == 0) THEN PRINT *,' WITH SIMPLE PACKING' IF (IAND(kptr(14),1) == 0) THEN PRINT *,' WITH NO ADDITIONAL FLAGS' GO TO 5000 ELSE IF (IAND(kptr(14),1) /= 0) THEN PRINT *,' WITH ADDITIONAL FLAGS' IF (kbds(17) == 0) THEN PRINT *,' SINGLE DATUM EACH GRID PT' IF (kbds(14) == 0) THEN PRINT *,' NO SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (kbds(14) /= 0) THEN PRINT *,' SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF END IF ELSE IF (kbds(17) /= 0) THEN PRINT *,' MATRIX OF VALS EACH PT' IF (kbds(14) == 0) THEN PRINT *,' NO SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (kbds(14) /= 0) THEN PRINT *,' SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF END IF END IF END IF ELSE IF (IAND(kptr(14),4) /= 0) THEN ! COMPLEX/SECOND ORDER PACKING PRINT *,' WITH COMPLEX/SECOND ORDER PACKING' IF (IAND(kptr(14),1) == 0) THEN PRINT *,' WITH NO ADDITIONAL FLAGS' ELSE IF (IAND(kptr(14),1) /= 0) THEN PRINT *,' WITH ADDITIONAL FLAGS' IF (kbds(17) == 0) THEN PRINT *,' SINGLE DATUM EACH GRID PT' IF (kbds(14) == 0) THEN PRINT *,' NO SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (kbds(14) /= 0) THEN PRINT *,' SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF END IF ELSE IF (kbds(17) /= 0) THEN PRINT *,' MATRIX OF VALS EACH PT' IF (kbds(14) == 0) THEN PRINT *,' NO SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (kbds(14) /= 0) THEN PRINT *,' SEC BIT MAP' IF (kbds(16) == 0) THEN PRINT *,' SECOND ORDER', & ' VALUES CONSTANT WIDTH' ELSE IF (kbds(16) /= 0) THEN PRINT *,' SECOND ORDER', & ' VALUES DIFFERENT WIDTHS' END IF END IF END IF END IF END IF END IF PRINT *,' NOT PROCESSED - NOT PROCESSED - NOT PROCESSED' kret = 11 RETURN 4000 CONTINUE ! **************************************************************** ! ! GRID POINT DATA, SIMPLE PACKING, FLOATING POINT, NO ADDN'L FLAGS ! scal10 = 10.0 ** kpds(22) scal10 = 1.0 / scal10 IF (kpds(3) == 23.OR.kpds(3) == 24.OR.kpds(3) == 26 & .OR.kpds(3) == 63.OR.kpds(3) == 64) THEN IF (kpds(3) == 26) THEN kadd = 72 ELSE IF (kpds(3) == 63.OR.kpds(3) == 64) THEN kadd = 91 ELSE kadd = 37 END IF CALL grbgbytes (msga,ksave,kptr(8),kbits,0,knr) kptr(8) = kptr(8) + kbits * knr ii = 1 kentry = kptr(10) DO i = 1, kentry IF (kbms(i)) THEN DATA(i) = (refnce+FLOAT(ksave(ii))*scale)*scal10 ii = ii + 1 ELSE DATA(i) = 0.0 END IF END DO DO i = 2, kadd DATA(i) = DATA(1) END DO ELSE IF (kpds(3) == 21.OR.kpds(3) == 22.OR.kpds(3) == 25 & .OR.kpds(3) == 61.OR.kpds(3) == 62) THEN CALL grbgbytes (msga,ksave,kptr(8),kbits,0,knr) ii = 1 kentry = kptr(10) DO i = 1, kentry IF (kbms(i)) THEN DATA(i) = (refnce + FLOAT(ksave(ii)) * scale) * scal10 ii = ii + 1 ELSE DATA(i) = 0.0 END IF END DO IF (kpds(3) == 25) THEN kadd = 71 ELSE IF (kpds(3) == 61.OR.kpds(3) == 62) THEN kadd = 90 ELSE kadd = 36 END IF lastp = kentry - kadd DO i = lastp+1, kentry DATA(i) = DATA(lastp) END DO ELSE CALL grbgbytes (msga,ksave,kptr(8),kbits,0,knr) ii = 1 kentry = kptr(10) DO i = 1, kentry IF (kbms(i)) THEN DATA(i) = (refnce + FLOAT(ksave(ii)) * scale) * scal10 ii = ii + 1 ELSE DATA(i) = 0.0 END IF END DO END IF GO TO 900 ! ------------- PROCESS SPHERICAL HARMONIC COEFFICIENTS, ! SIMPLE PACKING, FLOATING POINT, NO ADDN'L FLAGS 5000 CONTINUE ! PRINT *,'CHECK POINT SPECTRAL COEFF' kptr(8) = ibyt12 CALL grbgbyte (msga,kkk,kptr(8),32) kptr(8) = kptr(8) + 32 ! ! THE NEXT CODE WILL CONVERT THE IBM370 FOATING POINT ! TO THE FLOATING POINT USED ON YOUR MACHINE. ! ! 1ST TEST TO SEE IN ON 32 OR 64 BIT WORD MACHINE ! LW = 4 OR 8; IF 8 MAY BE A CRAY ! CALL w3fi01(lw) IF (lw == 4) THEN CALL grbgbyte (kk,jsgn,0,1) CALL grbgbyte (kk,jexp,1,7) CALL grbgbyte (kk,ifr,8,24) ELSE CALL grbgbyte (kk,jsgn,32,1) CALL grbgbyte (kk,jexp,33,7) CALL grbgbyte (kk,ifr,40,24) END IF ! IF (ifr == 0) THEN realkk = 0.0 ELSE IF (jexp == 0.AND.ifr == 0) THEN realkk = 0.0 ELSE realkk = FLOAT(ifr) * 16.0 ** (jexp - 64 - 6) IF (jsgn /= 0) realkk = -realkk END IF DATA(1) = realkk CALL grbgbytes (msga,ksave,kptr(8),kbits,0,knr) ! -------------- DO i = 1, kentry DATA(i+1) = refnce + FLOAT(ksave(i)) * scale END DO 900 CONTINUE ! PRINT *,'EXIT FI635' RETURN END SUBROUTINE fi635 SUBROUTINE fi636 (DATA,msga,kbms,refnce,kptr,kpds,kgds) 2,17 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: FI636 PROCESS SECOND ORDER PACKING ! PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 92-09-22 ! ! ABSTRACT: PROCESS SECOND ORDER PACKING FROM THE BINARY DATA SECTION ! (BDS) FOR SINGLE DATA ITEMS GRID POINT DATA ! ! PROGRAM HISTORY LOG: ! 93-06-08 CAVANAUGH ! 93-12-15 CAVANAUGH MODIFIED SECOND ORDER POINTERS TO FIRST ORDER ! VALUES AND SECOND ORDER VALUES CORRECTLY. ! ! USAGE: CALL FI636 (DATA,MSGA,KBMS,REFNCE,KPTR,KPDS,KGDS) ! INPUT ARGUMENT LIST: ! ! MSGA - ARRAY CONTAINING GRIB MESSAGE ! REFNCE - REFERENCE VALUE ! KPTR - WORK ARRAY ! ! OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) ! DATA - LOCATION OF OUTPUT ARRAY ! WORKING ARRAY ! KBDS(1) - N1 ! KBDS(2) - N2 ! KBDS(3) - P1 ! KBDS(4) - P2 ! KBDS(5) - BIT POINTER TO 2ND ORDER WIDTHS ! KBDS(6) - " " " " " BIT MAPS ! KBDS(7) - " " " FIRST ORDER VALUES ! KBDS(8) - " " " SECOND ORDER VALUES ! KBDS(9) - " " START OF BDS ! KBDS(10) - " " MAIN BIT MAP ! KBDS(11) - BINARY SCALING ! KBDS(12) - DECIMAL SCALING ! KBDS(13) - BIT WIDTH OF FIRST ORDER VALUES ! KBDS(14) - BIT MAP FLAG ! 0 = NO SECOND ORDER BIT MAP ! 1 = SECOND ORDER BIT MAP PRESENT ! KBDS(15) - SECOND ORDER BIT WIDTH ! KBDS(16) - CONSTANT / DIFFERENT WIDTHS ! 0 = CONSTANT WIDTHS ! 1 = DIFFERENT WIDTHS ! KBDS(17) - SINGLE DATUM / MATRIX ! 0 = SINGLE DATUM AT EACH GRID POINT ! 1 = MATRIX OF VALUES AT EACH GRID POINT ! (18-20)- UNUSED ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 77 ! MACHINE: NAS, CYBER ! !$$$ REAL :: DATA(*) REAL :: refn REAL :: refnce ! INTEGER :: kbds(20) INTEGER :: kptr(*) INTEGER :: jref,bmap2(12500) INTEGER :: i,ibds INTEGER :: kbit,ifoval,isoval INTEGER :: kpds(*),kgds(*) ! LOGICAL*1 kbms(*) ! CHARACTER :: msga(*) ! EQUIVALENCE (jref,refn) ! ******************* SETUP ****************************** SAVE ! PRINT *,'ENTER FI636' ij = 0 ! ======================================================== IF (kbds(14) == 0) THEN ! NO BIT MAP, MUST CONSTRUCT ONE IF (kgds(2) == 65535) THEN IF (kgds(20) == 255) THEN PRINT *,'CANNOT BE USED HERE' ELSE ! POINT TO PL lp = kptr(9) + kptr(2)*8 + kptr(3)*8 + kgds(20)*8 - 8 jt = 0 DO jz = 1, kgds(3) ! GET NUMBER IN CURRENT ROW CALL grbgbyte (msga,NUMBER,lp,16) ! INCREMENT TO NEXT ROW NUMBER lp = lp + 16 ! PRINT *,'NUMBER IN ROW',JZ,' = ',NUMBER DO jq = 1, NUMBER jt = jt + 1 IF (jq == 1) THEN bmap2(jt) = 1 ELSE bmap2(jt) = 0 END IF END DO END DO END IF ELSE IF (IAND(kgds(11),32) == 0) THEN ! ROW BY ROW ! PRINT *,' ROW BY ROW' kout = kgds(3) kin = kgds(2) ELSE ! COL BY COL ! PRINT *,' COL BY COL' kin = kgds(3) kout = kgds(2) END IF ! PRINT *,'KIN=',KIN,' KOUT= ',KOUT DO i = 1, kout DO j = 1, kin IF (j == 1) THEN CALL grbsbyte (bmap2,1,ij,1) ELSE CALL grbsbyte (bmap2,0,ij,1) END IF ij = ij + 1 END DO END DO END IF END IF ! ======================================================== ! CALL BINARY (BMAP2,2) ! START OF BMS (BIT POINTER) ! KBDS(10) = 0 ! BYTE START OF BDS ibds = kptr(2) + kptr(3) + kptr(4) + kptr(5) ! PRINT *,'KPTR(2-5) ',KPTR(2),KPTR(3),KPTR(4),KPTR(5) ! BIT START OF BDS jptr = ibds * 8 ! PRINT *,'JPTR ',JPTR kbds(9) = jptr ! PRINT *,'START OF BDS ',KBDS(9) ! BINARY SCALE VALUE BDS BYTES 5-6 CALL grbgbyte (msga,ISIGN,jptr+32,1) CALL grbgbyte (msga,kbds(11),jptr+33,15) IF (ISIGN > 0) THEN kbds(11) = - kbds(11) END IF ! PRINT *,'BINARY SCALE VALUE =',KBDS(11) ! EXTRACT REFERENCE VALUE CALL grbgbyte(msga,jref,jptr+48,32) ! PRINT *,'DECODED REFERENCE VALUE =',REFN,REFNCE ! F O BIT WIDTH CALL grbgbyte(msga,kbds(13),jptr+80,8) jptr = jptr + 88 ! AT START OF BDS BYTE 12 ! EXTRACT N1 CALL grbgbyte (msga,kbds(1),jptr,16) ! PRINT *,'N1 = ',KBDS(1) jptr = jptr + 16 ! EXTENDED FLAGS CALL grbgbyte (msga,kflag,jptr,8) ! ISOLATE BIT MAP FLAG kbds(14) = IAND(kflag,32) jptr = jptr + 8 ! EXTRACT N2 CALL grbgbyte (msga,kbds(2),jptr,16) ! PRINT *,'N2 = ',KBDS(2) jptr = jptr + 16 ! EXTRACT P1 CALL grbgbyte (msga,kbds(3),jptr,16) ! PRINT *,'P1 = ',KBDS(3) jptr = jptr + 16 ! EXTRACT P2 CALL grbgbyte (msga,kbds(4),jptr,16) ! PRINT *,'P2 = ',KBDS(4) jptr = jptr + 16 ! SKIP RESERVED BYTE jptr = jptr + 8 ! START OF SECOND ORDER BIT WIDTHS kbds(5) = jptr ! COMPUTE START OF SECONDARY BIT MAP IF (kbds(14) /= 0) THEN ! FOR INCLUDED SECONDARY BIT MAP jptr = jptr + (kbds(3) * 8) kbds(6) = jptr ELSE ! FOR CONSTRUCTED SECONDARY BIT MAP kbds(6) = 0 END IF ! CREATE POINTER TO START OF FIRST ORDER VALUES kbds(7) = kbds(9) + kbds(1) * 8 - 8 PRINT *,'BIT POINTER TO START OF FOVALS',kbds(7) ! CREATE POINTER TO START OF SECOND ORDER VALUES ! KBDS(8) = KBDS(9) + KBDS(2) * 8 - 8 ! PRINT *,'BIT POINTER TO START OF SOVALS',KBDS(8) ! PRINT *,'KBDS( 1) - N1 ',KBDS( 1) ! PRINT *,'KBDS( 2) - N2 ',KBDS( 2) ! PRINT *,'KBDS( 3) - P1 ',KBDS( 3) ! PRINT *,'KBDS( 4) - P2 ',KBDS( 4) ! PRINT *,'KBDS( 5) - BIT PTR - 2ND ORDER WIDTHS ',KBDS( 5) ! PRINT *,'KBDS( 6) - " " " " BIT MAPS ',KBDS( 6) ! PRINT *,'KBDS( 7) - " " F O VALS ',KBDS( 7) ! PRINT *,'KBDS( 8) - " " S O VALS ',KBDS( 8) ! PRINT *,'KBDS( 9) - " " START OF BDS ',KBDS( 9) ! PRINT *,'KBDS(10) - " " MAIN BIT MAP ',KBDS(10) ! PRINT *,'KBDS(11) - BINARY SCALING ',KBDS(11) ! PRINT *,'KPDS(22) - DECIMAL SCALING ',KPDS(22) ! PRINT *,'KBDS(13) - FO BIT WIDTH ',KBDS(13) ! PRINT *,'KBDS(14) - 2ND ORDER BIT MAP FLAG ',KBDS(14) ! PRINT *,'KBDS(15) - 2ND ORDER BIT WIDTH ',KBDS(15) ! PRINT *,'KBDS(16) - CONSTANT/DIFFERENT WIDTHS ',KBDS(16) ! PRINT *,'KBDS(17) - SINGLE DATUM/MATRIX ',KBDS(17) ! PRINT *,'REFNCE VAL ',REFNCE ! ************************* PROCESS DATA ********************** ! FOR EACH GRID POINT ENTRY ! scale2 = 2.0**kbds(11) scal10 = 10.0**kpds(22) ! PRINT *,'SCALE VALUES - ',SCALE2,SCAL10 DO i = 1, kptr(10) ! GET NEXT MASTER BIT MAP BIT POSITION ! IF NEXT MASTER BIT MAP BIT POSITION IS 'ON' (1) IF (kbms(i)) THEN ! WRITE(6,900)I,KBMS(I) ! 900 FORMAT (1X,I4,3X,14HMAIN BIT IS ON,3X,L4) IF (kbds(14) /= 0) THEN CALL grbgbyte (msga,kbit,kbds(6),1) ELSE CALL grbgbyte (bmap2,kbit,kbds(6),1) END IF ! PRINT *,'KBDS(6) =',KBDS(6),' KBIT =',KBIT kbds(6) = kbds(6) + 1 IF (kbit /= 0) THEN ! PRINT *,' SOB ON' ! GET NEXT FIRST ORDER PACKED VALUE CALL grbgbyte (msga,ifoval,kbds(7),kbds(13)) kbds(7) = kbds(7) + kbds(13) ! PRINT *,'FOVAL =',IFOVAL ! GET SECOND ORDER BIT WIDTH CALL grbgbyte (msga,kbds(15),kbds(5),8) kbds(5) = kbds(5) + 8 ! PRINT *,KBDS(7)-KBDS(13),' FOVAL =',IFOVAL,' KBDS(5)=', ! * ,KBDS(5), 'ISOWID =',KBDS(15) ELSE ! PRINT *,' SOB NOT ON' END IF isoval = 0 IF (kbds(15) == 0) THEN ! IF SECOND ORDER BIT WIDTH = 0 ! THEN SECOND ORDER VALUE IS 0 ! SO CALCULATE DATA VALUE FOR THIS POINT ! DATA(I) = (REFNCE + (FLOAT(IFOVAL) * SCALE2)) / SCAL10 ELSE CALL grbgbyte (msga,isoval,kbds(8),kbds(15)) kbds(8) = kbds(8) + kbds(15) END IF DATA(i) = (refnce + (FLOAT(ifoval + isoval) * & scale2)) / scal10 ! PRINT *,I,DATA(I),REFNCE,IFOVAL,ISOVAL,SCALE2,SCAL10 ELSE ! WRITE(6,901) I,KBMS(I) ! 901 FORMAT (1X,I4,3X,15HMAIN BIT NOT ON,3X,L4) DATA(i) = 0.0 END IF ! PRINT *,I,DATA(I),IFOVAL,ISOVAL,KBDS(5),KBDS(15) END DO ! ************************************************************** ! PRINT *,'EXIT FI636' RETURN END SUBROUTINE fi636 SUBROUTINE fi637(*,j,kpds,kgds,kret) 10 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: FI637 GRIB GRID/SIZE TEST ! PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 ! ! ABSTRACT: TO TEST WHEN GDS IS AVAILABLE TO SEE IF SIZE MISMATCH ! ON EXISTING GRIDS (BY CENTER) IS INDICATED ! ! PROGRAM HISTORY LOG: ! 91-09-13 CAVANAUGH ! ! USAGE: CALL FI637(*,J,KPDS,KGDS,KRET) ! INPUT ARGUMENT LIST: ! J - SIZE FOR INDICATED GRID ! KPDS - ! KGDS - ! ! OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) ! KRET - ERROR RETURN ! ! REMARKS: ! KRET - ! = 9 - GDS INDICATES SIZE MISMATCH WITH STD GRID ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 77 ! MACHINE: NAS ! !$$$ INTEGER :: kpds(*) INTEGER :: kgds(*) INTEGER :: j INTEGER :: i ! --------------------------------------- SAVE ! --------------------------------------- ! IF GDS NOT INDICATED, RETURN ! ---------------------------------------- IF (IAND(kpds(4),128) == 0) RETURN ! --------------------------------------- ! GDS IS INDICATED, PROCEED WITH TESTING ! --------------------------------------- IF (kgds(2) == 65535) THEN RETURN END IF i = kgds(2) * kgds(3) ! --------------------------------------- ! INTERNATIONAL SET ! --------------------------------------- IF (kpds(3) >= 21.AND.kpds(3) <= 26) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) >= 37.AND.kpds(3) <= 44) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) == 50) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) >= 61.AND.kpds(3) <= 64) THEN IF (i /= j) THEN RETURN 1 END IF ! --------------------------------------- ! TEST ECMWF CONTENT ! --------------------------------------- ELSE IF (kpds(1) == 98) THEN kret = 9 IF (kpds(3) >= 1.AND.kpds(3) <= 16) THEN IF (i /= j) THEN RETURN 1 END IF ELSE kret = 5 RETURN 1 END IF ! --------------------------------------- ! U.K. MET OFFICE, BRACKNELL ! --------------------------------------- ELSE IF (kpds(1) == 74) THEN kret = 9 IF (kpds(3) >= 25.AND.kpds(3) <= 26) THEN IF (i /= j) THEN RETURN 1 END IF ELSE kret = 5 RETURN 1 END IF ! --------------------------------------- ! CANADA ! --------------------------------------- ELSE IF (kpds(1) == 54) THEN PRINT *,' NO CURRENT LISTING OF CANADIAN GRIDS' RETURN 1 ! --------------------------------------- ! JAPAN METEOROLOGICAL AGENCY ! --------------------------------------- ELSE IF (kpds(1) == 34) THEN PRINT *,' NO CURRENT LISTING OF JMA GRIDS' RETURN 1 ! --------------------------------------- ! NAVY - FNOC ! --------------------------------------- ELSE IF (kpds(1) == 58) THEN IF (kpds(3) >= 220.AND.kpds(3) <= 221) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) == 223) THEN IF (i /= j) THEN RETURN 1 END IF END IF ! --------------------------------------- ! U.S. GRIDS ! --------------------------------------- ELSE IF (kpds(1) == 7) THEN kret = 9 IF (kpds(3) >= 1.AND.kpds(3) <= 3) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) == 5.OR.kpds(3) == 6) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) >= 27.AND.kpds(3) <= 30) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) >= 33.AND.kpds(3) <= 34) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) >= 37.AND.kpds(3) <= 45) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) >= 55.AND.kpds(3) <= 56) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) >= 67.AND.kpds(3) <= 77) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) >= 85.AND.kpds(3) <= 86) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) == 87) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) >= 90.AND.kpds(3) <= 95) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) == 100.OR.kpds(3) == 101) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) >= 103.AND.kpds(3) <= 107) THEN IF (i /= j) THEN RETURN 1 END IF ELSE IF (kpds(3) >= 201.AND.kpds(3) <= 214) THEN IF (i /= j) THEN RETURN 1 END IF ELSE kret = 5 RETURN 1 END IF ELSE kret = 10 RETURN 1 END IF ! ------------------------------------ ! NORMAL EXIT ! ------------------------------------ kret = 0 RETURN END SUBROUTINE fi637 SUBROUTINE w3fi83 (DATA,npts,fval1,fdiff1,iscal2, & 1 isc10,kpds,kgds) !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: W3FI83 RESTORE DELTA PACKED DATA TO ORIGINAL ! PRGMMR: CAVANAUGH ORG: NMC421 DATE:94-03-09 ! ! ABSTRACT: RESTORE DELTA PACKED DATA TO ORIGINAL VALUES ! RESTORE FROM BOUSTREPHEDONIC ALIGNMENT ! ! PROGRAM HISTORY LOG: ! 93-07-14 CAVANAUGH ! 93-07-22 STACKPOLE ADDITIONS TO FIX SCALING ! 94-01-27 CAVANAUGH ADDED REVERSAL OF EVEN NUMBERED ROWS ! (BOUSTROPHEDONIC PROCESSING) TO RESTORE ! DATA TO ORIGINAL SEQUENCE. ! 94-03-02 CAVANAUGH CORRECTED REVERSAL OF EVEN NUMBERED ROWS ! ! USAGE: CALL W3FI83(DATA,NPTS,FVAL1,FDIFF1,ISCAL2, ! * ISC10,KPDS,KGDS) ! INPUT ARGUMENT LIST: ! DATA - SECOND ORDER DIFFERENCES ! NPTS - NUMBER OF POINTS IN ARRAY ! FVAL1 - ORIGINAL FIRST ENTRY IN ARRAY ! FDIFF1 - ORIGINAL FIRST FIRST-DIFFERENCE ! ISCAL2 - POWER-OF-TWO EXPONENT FOR UNSCALING ! ISC10 - POWER-OF-TEN EXPONENT FOR UNSCALING ! KPDS - ARRAY OF INFORMATION FOR PDS ! KGDS - ARRAY OF INFORMATION FOR GDS ! ! OUTPUT ARGUMENT LIST: ! DATA - EXPANDED ORIGINAL DATA VALUES ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 77 ! MACHINE: NAS ! !$$$ REAL :: fval1,fdiff1 REAL :: DATA(*),boust(200) INTEGER :: npts,nrow,ncol,kpds(*),kgds(*),isc10 ! --------------------------------------- SAVE ! ! REMOVE DECIMAL UN-SCALING INTRODUCED DURING UNPACKING ! dscal = 10.0 ** isc10 IF (dscal == 0.0) THEN DO i=1,npts DATA(i) = 1.0 END DO ELSE IF (dscal == 1.0) THEN ELSE DO i=1,npts DATA(i) = DATA(i) * dscal END DO END IF ! DATA(1) = fval1 DATA(2) = fdiff1 DO j = 3,2,-1 DO k = j, npts DATA(k) = DATA(k) + DATA(k-1) END DO END DO ! ! NOW REMOVE THE BINARY SCALING FROM THE RECONSTRUCTED FIELD ! AND THE DECIMAL SCALING TOO ! IF (dscal == 0) THEN scale = 0.0 ELSE scale =(2.0**iscal2)/dscal END IF DO i=1,npts DATA(i) = DATA(i) * scale END DO ! ========================================================== IF (IAND(kpds(4),128) /= 0) THEN nrow = kgds(3) ncol = kgds(2) ! ! DATA LAID OUT BOUSTROPHEDONIC STYLE ! ! ! PRINT*, ' REVERSE BOUSTROPHEDON' DO i = 2, nrow, 2 ! ! REVERSE THE EVEN NUMBERED ROWS ! DO j = 1, ncol npos = i * ncol - j + 1 boust(j) = DATA(npos) END DO DO j = 1, ncol npos = ncol * (i-1) + j DATA(npos) = boust(j) END DO END DO ! ! END IF ! ================================================================= RETURN END SUBROUTINE w3fi83 SUBROUTINE w3fi01(lw) 5 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: W3FI01 DETERMINES MACHINE WORD LENGTH IN BYTES ! PRGMMR: KEYSER ORG: NMC22 DATE:92-07-30 ! ! ABSTRACT: DETERMINES THE NUMBER OF BYTES IN A FULL WORD FOR THE ! PARTICULAR MACHINE (IBM, CRAY, WORKSTATION, PC). ! ! PROGRAM HISTORY LOG: ! 92-01-10 R. KISTLER (W/NMC23) ! 92-05-22 D. A. KEYSER -- DOCBLOCKED/COMMENTED ! 93-03-29 R.E.JONES -- ADD SAVE STATEMENT ! ! USAGE: CALL W3FI01(LW) ! OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) ! LW - MACHINE WORD LENGTH IN BYTES ! ! OUTPUT FILES: ! FT06F001 - PRINTOUT ! ! REMARKS: NONE. ! ! ATTRIBUTES: ! LANGUAGE: IBM AIX XL FORTRAN Compiler/6000 ! MACHINE: IBM RS6000 model 530 ! !$$$ ! CHARACTER (LEN=8) :: ctest1,ctest2 ! INTEGER :: itest1,itest2 ! EQUIVALENCE (ctest1,itest1),(ctest2,itest2) ! SAVE ! DATA ctest1/'12345678'/ ! itest2 = itest1 !CCCC PRINT *,' CTEST1 = ',CTEST1,' CTEST2 = ',CTEST2 IF (ctest1 == ctest2) THEN !CCCC PRINT*,' MACHINE WORD LENGTH IS 8 BYTES' lw = 8 ELSE !CCCC PRINT*,' MACHINE WORD LENGTH IS 4 BYTES' lw = 4 END IF RETURN END SUBROUTINE w3fi01 SUBROUTINE w3fi04(iendn,itypec,lw) !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: W3FI04 FIND WORD SIZE, ENDIAN, CHARACTER SET ! PRGMNR: JONES,R.E. ORG: W/NMC42 DATE: 94-10-07 ! ! ABSTRACT: SUBROUTINE COMPUTES WORD SIZE, THE TYPE OF CHARACTER ! SET, ASCII OR EBCDIC, AND IF THE COMPUTER IS BIG-ENDIAN, OR ! LITTLE-ENDIAN. ! ! PROGRAM HISTORY LOG: ! 94-10-07 R.E.JONES ! ! USAGE: CALL W3FI04 (IENDN, ITYPEC, LW) ! ! OUTPUT ARGUMENT LIST: ! IENDN - INTEGER FOR BIG-ENDIAN OR LITTLE-ENDIAN ! = 0 BIG-ENDIAN ! = 1 LITTLE-ENDIAN ! = 2 CANNOT COMPUTE ! ITYPEC - INTEGER FOR TYPE OF CHARACTER SET ! = 0 ASCII CHARACTER SET ! = 1 EBCDIC CHARACTER SET ! = 2 NOT ASCII OR EBCDIC ! LW - INTEGER FOR WORDS SIZE OF COMPUTER IN BYTES ! = 4 FOR 32 BIT COMPUTERS ! = 8 FOR 64 BIT COMPUTERS ! ! ATTRIBUTES: ! LANGUAGE: SiliconGraphics 3.5 FORTRAN 77 ! MACHINE: SiliconGraphics IRIS-4D/25, 35, INDIGO, INDY ! !$$$ ! INTEGER :: itest1 INTEGER :: itest2 INTEGER :: itest3 INTEGER :: iendn INTEGER :: itypec INTEGER :: lw ! CHARACTER (LEN=8) :: ctest1 CHARACTER (LEN=8) :: ctest2 CHARACTER (LEN=1) :: ctest3(8) CHARACTER (LEN=1) :: BLANK ! EQUIVALENCE (ctest1,itest1),(ctest2,itest2) ! EQUIVALENCE (itest3,ctest3(1)) ! DATA ctest1/'12345678'/ DATA itest3/z'01020304'/ DATA BLANK /' '/ ! SAVE ! ! TEST FOR TYPE OF CHARACTER SET ! BLANK IS 32 (20 HEX) IN ASCII, 64 (40 HEX) IN EBCDEC ! IF (ICHAR(BLANK) == 32) THEN itypec = 0 ELSE IF (ICHAR(BLANK) == 64) THEN ! ! COMPUTER IS PROBABLY AN IBM360, 370, OR 390 WITH ! A 32 BIT WORD SIZE, AND BIG-ENDIAN. ! itypec = 1 ELSE itypec = 2 END IF ! ! TEST FOR WORD SIZE, SET LW TO 4 FOR 32 BIT COMPUTER, ! 8 FOR FOR 64 BIT COMPUTERS ! itest2 = itest1 IF (ctest1 == ctest2) THEN ! ! COMPUTER MAY BE A CRAY, OR COULD BE DEC VAX ALPHA ! OR SGI WITH R4000, R4400, R8800 AFTER THEY CHANGE ! FORTRAN COMPILERS FOR 64 BIT INTEGER. ! lw = 8 ELSE lw = 4 END IF ! ! IF CHARACTER SET IS NOT ASCII OR EBCDIC SET IENDN = 2 ! CAN NOT TEST FOR ENDIAN TYPE ! IF (itypec == 2) THEN iendn = 2 RETURN END IF ! ! USING ITEST3 WITH Z'01020304' EQUIVALNCED TO CTEST3 ! ON A 32 BIT BIG-ENDIAN COMPUTER 03 IS IN THE 3RD ! BYTE OF A 4 BYTE WORD. ON A 32 BIT LITTLE-ENDIAN ! COMPUTER IT IS IN 2ND BYTE. ! ON A 64 BIT COMPUTER Z'01020304' IS RIGHT ADJUSTED IN ! A 64 BIT WORD, 03 IS IN THE 7TH BYTE. ON A LITTLE- ! ENDIAN 64 BIT COMPUTER IT IS IN THE 2ND BYTE. ! IF (lw == 4) THEN IF (ICHAR(ctest3(3)) == 3) THEN iendn = 0 ELSE IF (ICHAR(ctest3(3)) == 2) THEN iendn = 1 ELSE iendn = 2 END IF ELSE IF (lw == 8) THEN IF (ICHAR(ctest3(7)) == 3) THEN iendn = 0 ELSE IF (ICHAR(ctest3(2)) == 3) THEN iendn = 1 ELSE iendn = 2 END IF ELSE iendn = 2 END IF ! RETURN END SUBROUTINE w3fi04