PROGRAM PLOTPP C$$$ MAIN PROGRAM DOCUMENTATION BLOCK C C C THIS CODE IS NOT "READY TO GO" RIGHT OUT OF THE BLOCKS C THE USER WILL PROBABLY HAVE TO MODIFY IT TO WORK ON THEIR MACHINE C PARTICULARLY IN THE AREA OF READING IN THE GRIB MESSAGE C FOR CRAYS, MAKE SURE UNIT 11 IS ASSIGNED AS -s unblocked C C THIS CODE USES NCAR GRAPHICS FOR PLOTTING C THIS CODE USES NCEP GRIB DECODERS, ACCESS THOSE ON THE nic C IN /pub/info IN THE SUBDIRECTORY gribYYYY WHERE YYYY IS C THE MACHINE YOU ARE USING. COMPILE THE W3LIB FOR YOUR C MACHINE AND MAKE SURE THAT MXSIZE HAS BEEN INCREASED IN W3FI63 C TO AT LEAST 11000000. C C C MAIN PROGRAM: PLOTPP PLOT GRIB PRECIP ANALYSIS FILE C AUTHOR: M. BALDWIN ORG: GSC DATE: 96-10-08 C C ABSTRACT: ASSUME ONE GRIB RECORD IN FILE. READ GRIB MESSAGE, C UNPACK, PLOT ENTIRE GRID. C C PROGRAM HISTORY LOG: C 96-10-08 M. BALDWIN ORIG C 97-02-19 M. BALDWIN FIX BUGS WITH OUT-OF-BOUNDS REFERENCES C (THANKS NANCY WESTCOTT, ILL WATER SURVEY) C C USAGE: C INPUT FILES: C 5 - CONTROL FILE C INDICATES NAME OF INPUT FILE C 11 - GRIB FILE USED FOR INPUT C C OUTPUT FILES: C 6 - PRINT OUTPUT (STANDARD FORTRAN OUTPUT FILE) C 8 - NCAR GRAPHICS GMETA FILE C C SUBPROGRAMS CALLED: C UNIQUE: NEWPEN, W3FB07 C C LIBRARY: C W3LIB - W3FI63, GBYTES, SBYTES, OTHER STUFF FOR DEGRIB C NCARGRAPHICS 4.0 - C C EXIT STATES: C C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C MACHINE: Cray J90 C C$$$ C C MXSIZE IS SET FOR LARGEST GRIB RECORD, 4KM PRECIP GRID C (1160 * 880) GRID. THIS MAY BE CHANGED TO THE C LARGEST RECORD SIZE YOU USE. C PARAMETER (MXSIZE=1100000,IMX=1160,JMX=880) PARAMETER (MSIZE=MXSIZE*3) C CHARACTER * 1 MSGA(MSIZE) CHARACTER * 3 STR3 CHARACTER * 80 FILNAM CHARACTER * 44 TITLE C INTEGER KPTR(20) INTEGER KPDS(100) INTEGER KGDS(200) INTEGER JSTAT(13) C REAL DATA(MXSIZE),THRESH(15) C LOGICAL*1 LBMS(MXSIZE) LOGICAL*1 VAX,IMISS(IMX,JMX) C SAVE C EQUIVALENCE (MSGA(1),IMESG) C DATA KGDS / 200 * 0 / DATA KPDS / 100 * 0 / DATA KPTR / 20 * 0 / DATA OUTPUT/6/ DATA THRESH/0.1,2,5,10,15,20,25,50,75,100,125,150,200,250,999/ COMMON /MAPCM2/ UMIN,UMAX,VMIN,VMAX,UEPS,VEPS,UCEN,VCEN,URNG,VRNG, + BLAM,SLAM,BLOM,SLOM,ISSL C C READ INPUT FILE NAME C FIND OUT HOW BIG IT IS C READ GRIB MESSAGE C MOST OF THIS CODE WAS STOLEN FROM UNPKGRB.F C READ(5,333,END=999) FILNAM 333 FORMAT(A80) PRINT *,' READING ',FILNAM C C CALL FUNCTION STAT TO FIND NUMBER OF BYTES IN FILE C WORD 8 OF ARRAY JSTAT HAS THE NUMBER OF BYTES IN FILE C IF YOUR COMPUTER DOES NOT HAVE THE C FUNCTION STAT C COMMENT OUT THE NEXT 7 LINES C IF (STAT(FILNAM,JSTAT).NE.0) THEN PRINT *,'ERROR IN FUNCTION STAT GETTING FILE STATS' STOP 99 ELSE KBYTES = JSTAT(8) PRINT *,'NUMBER OF BYTES IN GRIB FILE = ',KBYTES END IF C C TEST TO SEE IF INPUT GRIB FILE IS TO BIG C IF (MSIZE.LT.KBYTES) THEN PRINT *,'GRIB INPUT FILE IS TO BIG ' PRINT *,'CHANGE PROGRAM SO PARAMETER MXSIZE IS LARGER' PRINT *,'THAN THE NUMBER OF BYTES IN THE FILE' STOP 40 END IF C C OPEN FILE WITH GRIB COM SEQUENTIAL FILE OR AWIPS FILE C On some computers RECL can be set to KBYTES, the number C of bytes in the file. This is both compiler and computer C dependent. On the 1st read, the entire file is read into C array MSGA. SET RECL=KBYTES FOR IBM RS/6000, HP 700. C C OPEN (UNIT=11,FILE=FILNAM,STATUS='OLD',ACCESS='DIRECT', C & FORM='UNFORMATTED',IOSTAT=MERR,RECL=KBYTES) C OPEN (UNIT=11,FILE=FILNAM,STATUS='OLD',ACCESS='DIRECT', & FORM='UNFORMATTED',IOSTAT=MERR,RECL=1) IF (MERR.NE.0) THEN PRINT *,'OPEN INPUT FILE ERROR ON FILE = ', FILNAM PRINT *,'ERROR = ',MERR STOP 20 END IF C C Read the entire GRIB file into array MSGA C On SUN and SGI the stat function works, when RECL=1 C REC is a pointer at what byte in the file to start C reading at. The SUN and SGI compilers compile this as C one big read, other compilers, such as the INTERGRAPH, C give a compile error, or run time error. C If so, comment out this line, uncomment the code that C that applies to your computer. This works, but is very C slow because you are doing one read for each byte in C file. Replace this if your computer has a better way of C doing it. This should also work if RECL=KBYTES. C READ (11) (MSGA(I),I=1,KBYTES) C C C REVERSE ORDER OF BYTES IF VAX, PC, OR INTERGRAPH C IF (VAX) THEN IWORDS = (MBYTES+3) / 4 CALL SWAP32(MSGA,IWORDS) END IF C C DECODE GRIB MESSAGE C KPTR(10) = 0 CALL W3FI63(MSGA,KPDS,KGDS,LBMS,DATA,KPTR,IRET) C C SET UP GRID INFO C ALAT1=KGDS(4)*0.001 ALON1=KGDS(5)*0.001 DX=KGDS(8) ALONV=KGDS(7)*0.001 IM=KGDS(2) JM=KGDS(3) IYR=KPDS(8) IMN=KPDS(9) IDA=KPDS(10) IHR=KPDS(11) C C DO A PLOT NCAR GRAPHICS BEGINS HERE C CALL GOPKS(6,0) CALL GOPWK(8,9,1) CALL GACWK(8) call agseti('FRAME.',2) C C SET UP MAP CORNERS C CALL MAPSTI ('MV',8) CALL MAPSTC ('OU','PS') CALL MAPROJ('ST',90.,alonv,0.) call W3FB07(1.,1.,ALAT1,ALON1,DX,ALONV,ALATL,ALONL) xi=im*1.0 xj=jm*1.0 call W3FB07(xi,xj,ALAT1,ALON1,DX,ALONV,ALATR,ALONR) call newpen(2) call gstxci(1) call gsplci(1) CALL MAPSET('CO',alatl,alonl,alatr,alonr) CALL MAPSTI ('GR',0) CALL MAPINT CALL MAPPOS(.01,.99,.01,.90) CALL MAPDRW CALL MAPTRN(ALATL,ALONL,UVAL1,VVAL1) CALL MAPTRN(ALATR,ALONR,UVAL2,VVAL2) CALL SET(uval1,uval2,vval1,vval2,0.,1.,0.,1.,1) C C PLOT PRECIP, POINT BY POINT C write(6,3201) (kpds(kjk),kjk=1,22) write(6,3202) (kgds(kjk),kjk=1,20) 3201 format(' PDS ',9i8) 3202 format(' GDS ',9i8) xinc=1./(IM-1) yinc=1./(JM-1) do ijcol=1,14 icol=24-ijcol call gstxci(icol) call gsplci(icol) k=0 DO J=1,JM DO I=1,IM xi=(i-1)*xinc xj=(j-1)*yinc k=k+1 imiss(i,j)=lbms(k) vpre=data(k) if (imiss(i,j).and.vpre.lt.thresh(ijcol+1).and. & vpre.ge.thresh(ijcol)) CALL point(xi,xj) ENDDO ENDDO enddo c c plot edge of radar extent c call gstxci(7) call gsplci(7) DO J=2,JM-1 DO I=1,IM-1 xi=(i-1)*xinc xj=(j-1)*yinc if (imiss(i,j).and.(.not.imiss(i+1,j))) call point(xi,xj) if ((.not.imiss(i,j)).and.imiss(i+1,j)) call point(xi,xj) if (imiss(i,j).and.imiss(i+1,j).and.(.not.imiss(i,j+1))) & call point(xi,xj) if (imiss(i,j).and.imiss(i+1,j).and.(.not.imiss(i,j-1))) & call point(xi,xj) ENDDO ENDDO c c do title c CALL SET(.01,.99,.01,.99,0.,1.,0.,1.,1) call gstxci(1) call gsplci(1) 7887 CALL SETUSV('LW',1000) write(title,733) imn,ida,iyr,ihr+1 733 format(' pcp (mm) valid ', & i2.2,'/',i2.2,'/',i2.2,1x,i2.2,'00 GMT') CALL WTSTR(0.05,0.98,title,2,0,-1) c c do legend c udiff=0.001 vdiff=0.015 do ijcol=1,14 icol=24-ijcol call gstxci(icol) call gsplci(icol) y=0.051+((ijcol-1)/18)*0.03 ii=mod(ijcol,18) if (ii.eq.0) ii=18 x=ii*0.06 do kx=1,15 u=x-kx*udiff vb=y-vdiff*0.5 vt=y+vdiff*0.5 call line(u,vb,u,vt) enddo ighgt=nint(thresh(ijcol)) write(str3,3099) ighgt 3099 format(i3) call wtstr(x+0.018,y,str3,1,0,0) enddo CALL FRAME C 999 CALL GDAWK(8) CALL GCLWK(8) CALL GCLKS C END SUBROUTINE newpen(inp) COMMON/SPLT/ LNOP,NPLTDV common /compen/ ipx,xold,yold character*8 ctype common /comgks/ ctype DIMENSION red(32),green(32),blue(32) data red /0.000,1.000,0.000,0.000,1.000,0.000,1.000,0.545, & 1.000,1.000,1.000,1.000,0.933,0.545,0.804,0.933, & 1.000,0.804,1.000,0.933,0.498,0.000,0.000,0.063, & 0.118,0.000,0.000,0.537,0.569,0.545,1.000,0.000/ data green /0.106,0.000,1.000,0.000,1.000,1.000,0.000,0.278, & 0.510,0.647,0.682,0.416,0.173,0.000,0.000,0.251, & 0.498,0.522,0.843,0.933,1.000,0.804,0.545,0.306, & 0.565,0.698,0.933,0.408,0.173,0.000,1.000,0.000/ data blue /0.137,0.000,0.000,1.000,0.000,1.000,1.000,0.149, & 0.278,0.310,0.725,0.416,0.173,0.000,0.000,0.000, & 0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.545, & 1.000,0.933,0.933,0.804,0.933,0.545,1.000,0.000/ DATA ICALL/0/,ctype/'GKSCOLOR'/,pi/3.14159265/ rfac=2.0*pi/18. ICALL=1 call gsfais(1) call gscr(8,0,1.,1.,1.) do 100 ic = 1, 64 if (ic.le.32) then rr = red(ic) gg = green(ic) bb = blue(ic) else rx = (ic-8)*1.0 if (rx.le.18.) then rr = 0.25* (cos (rx*rfac))+0.75 gg = rx/18. bb = 0. else if (rx.gt.18.and.rx.le.24) then rr = 0.67 - (rx-18.)/9.0 gg = 1.0 - (rx-18.)/18. bb = 0. else if (rx.gt.24.and.rx.le.36) then rr = 0. gg = 1.0 - (rx-18.)/18. bb = (rx-18.0)*(rx-18.0)/324.0 else rr = (rx-36.)/20. gg = 0. bb = 0.25* (cos ((rx-36.)*rfac))+0.75 endif endif call gscr(8,ic,rr,gg,bb) 100 continue RETURN END SUBROUTINE W3FB07(XI,XJ,ALAT1,ALON1,DX,ALONV,ALAT,ALON) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: W3FB07 GRID COORDS TO LAT/LON FOR GRIB C PRGMMR: STACKPOLE ORG: NMC42 DATE:88-04-05 C C ABSTRACT: CONVERTS THE COORDINATES OF A LOCATION ON EARTH GIVEN IN A C GRID COORDINATE SYSTEM OVERLAID ON A POLAR STEREOGRAPHIC MAP PRO- C JECTION TRUE AT 60 DEGREES N OR S LATITUDE TO THE C NATURAL COORDINATE SYSTEM OF LATITUDE/LONGITUDE C W3FB07 IS THE REVERSE OF W3FB06. C USES GRIB SPECIFICATION OF THE LOCATION OF THE GRID C C PROGRAM HISTORY LOG: C 88-01-01 ORIGINAL AUTHOR: STACKPOLE, W/NMC42 C 90-04-12 R.E.JONES CONVERT TO CRAY CFT77 FORTRAN C C USAGE: CALL W3FB07(XI,XJ,ALAT1,ALON1,DX,ALONV,ALAT,ALON) C INPUT ARGUMENT LIST: C XI - I COORDINATE OF THE POINT REAL*4 C XJ - J COORDINATE OF THE POINT REAL*4 C ALAT1 - LATITUDE OF LOWER LEFT POINT OF GRID (POINT 1,1) C LATITUDE <0 FOR SOUTHERN HEMISPHERE; REAL*4 C ALON1 - LONGITUDE OF LOWER LEFT POINT OF GRID (POINT 1,1) C EAST LONGITUDE USED THROUGHOUT; REAL*4 C DX - MESH LENGTH OF GRID IN METERS AT 60 DEG LAT C MUST BE SET NEGATIVE IF USING C SOUTHERN HEMISPHERE PROJECTION; REAL*4 C 190500.0 LFM GRID, C 381000.0 NH PE GRID, -381000.0 SH PE GRID, ETC. C ALONV - THE ORIENTATION OF THE GRID. I.E., C THE EAST LONGITUDE VALUE OF THE VERTICAL MERIDIAN C WHICH IS PARALLEL TO THE Y-AXIS (OR COLUMNS OF C THE GRID) ALONG WHICH LATITUDE INCREASES AS C THE Y-COORDINATE INCREASES. REAL*4 C FOR EXAMPLE: C 255.0 FOR LFM GRID, C 280.0 NH PE GRID, 100.0 SH PE GRID, ETC. C C OUTPUT ARGUMENT LIST: C ALAT - LATITUDE IN DEGREES (NEGATIVE IN SOUTHERN HEMI.) C ALON - EAST LONGITUDE IN DEGREES, REAL*4 C C REMARKS: FORMULAE AND NOTATION LOOSELY BASED ON HOKE, HAYES, C AND RENNINGER'S "MAP PROJECTIONS AND GRID SYSTEMS...", MARCH 1981 C AFGWC/TN-79/003 C C ATTRIBUTES: C LANGUAGE: CRAY CFT77 FORTRAN C MACHINE: CRAY Y-MP8/832 C C$$$ C DATA RERTH /6.3712E+6/,PI/3.1416/ DATA SS60 /1.86603/ C C PRELIMINARY VARIABLES AND REDIFINITIONS C C H = 1 FOR NORTHERN HEMISPHERE; = -1 FOR SOUTHERN C C REFLON IS LONGITUDE UPON WHICH THE POSITIVE X-COORDINATE C DRAWN THROUGH THE POLE AND TO THE RIGHT LIES C ROTATED AROUND FROM ORIENTATION (Y-COORDINATE) LONGITUDE C DIFFERENTLY IN EACH HEMISPHERE C IF (DX.LT.0) THEN H = -1.0 DXL = -DX REFLON = ALONV - 90.0 ELSE H = 1.0 DXL = DX REFLON = ALONV - 270.0 ENDIF C RADPD = PI / 180.0 DEGPRD = 180.0 / PI REBYDX = RERTH / DXL C C RADIUS TO LOWER LEFT HAND (LL) CORNER C ALA1 = ALAT1 * RADPD RMLL = REBYDX * COS(ALA1) * SS60/(1. + H * SIN(ALA1)) C C USE LL POINT INFO TO LOCATE POLE POINT C ALO1 = (ALON1 - REFLON) * RADPD POLEI = 1. - RMLL * COS(ALO1) POLEJ = 1. - H * RMLL * SIN(ALO1) C C RADIUS TO THE I,J POINT (IN GRID UNITS) C XX = XI - POLEI YY = (XJ - POLEJ) * H R2 = XX**2 + YY**2 C C NOW THE MAGIC FORMULAE C IF (R2.EQ.0) THEN ALAT = H * 90. ALON = REFLON ELSE GI2 = (REBYDX * SS60)**2 ALAT = DEGPRD * H * ASIN((GI2 - R2)/(GI2 + R2)) ARCCOS = ACOS(XX/SQRT(R2)) IF (YY.GT.0) THEN ALON = REFLON + DEGPRD * ARCCOS ELSE ALON = REFLON - DEGPRD * ARCCOS ENDIF ENDIF IF (ALON.LT.0) ALON = ALON + 360. C RETURN END SUBROUTINE W3FI63(MSGA,KPDS,KGDS,KBMS,DATA,KPTR,KRET) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: W3FI63 UNPK GRIB FIELD TO GRIB GRID C PRGMMR: FARLEY ORG: NMC421 DATE:94-11-22 C C ABSTRACT: UNPACK A GRIB (EDITION 1) FIELD TO THE EXACT GRID C SPECIFIED IN THE GRIB MESSAGE, ISOLATE THE BIT MAP, AND MAKE C THE VALUES OF THE PRODUCT DESCRIPTON SECTION (PDS) AND THE C GRID DESCRIPTION SECTION (GDS) AVAILABLE IN RETURN ARRAYS. C C WHEN DECODING IS COMPLETED, DATA AT EACH GRID POINT HAS BEEN C RETURNED IN THE UNITS SPECIFIED IN THE GRIB MANUAL. C C PROGRAM HISTORY LOG: C 91-09-13 CAVANAUGH C 91-11-12 CAVANAUGH MODIFIED SIZE OF ECMWF GRIDS 5-8 C 91-12-22 CAVANAUGH CORRECTED PROCESSING OF MERCATOR PROJECTIONS C IN GRID DEFINITION SECTION (GDS) IN C ROUTINE FI633 C 92-08-05 CAVANAUGH CORRECTED MAXIMUM GRID SIZE TO ALLOW FOR C ONE DEGREE BY ONE DEGREE GLOBAL GRIDS C 92-08-27 CAVANAUGH CORRECTED TYPO ERROR, ADDED CODE TO COMPARE C TOTAL BYTE SIZE FROM SECTION 0 WITH SUM OF C SECTION SIZES. C 92-10-21 CAVANAUGH CORRECTIONS WERE MADE (IN FI634) TO REDUCE C PROCESSING TIME FOR INTERNATIONAL GRIDS. C REMOVED A TYPOGRAPHICAL ERROR IN FI635. C 93-01-07 CAVANAUGH CORRECTIONS WERE MADE (IN FI635) TO C FACILITATE USE OF THESE ROUTINES ON A PC. C A TYPOGRAPHICAL ERROR WAS ALSO CORRECTED C 93-01-13 CAVANAUGH CORRECTIONS WERE MADE (IN FI632) TO C PROPERLY HANDLE CONDITION WHEN C TIME RANGE INDICATOR = 10. C ADDED U.S.GRID 87. C 93-02-04 CAVANAUGH ADDED U.S.GRIDS 85 AND 86 C 93-02-26 CAVANAUGH ADDED GRIDS 2, 3, 37 THRU 44,AND C GRIDS 55, 56, 90, 91, 92, AND 93 TO C LIST OF U.S. GRIDS. C 93-04-07 CAVANAUGH ADDED GRIDS 67 THRU 77 TO C LIST OF U.S. GRIDS. C 93-04-20 CAVANAUGH INCREASED MAX SIZE TO ACCOMODATE C GAUSSIAN GRIDS. C 93-05-26 CAVANAUGH CORRECTED GRID RANGE SELECTION IN FI634 C FOR RANGES 67-71 & 75-77 C 93-06-08 CAVANAUGH CORRECTED FI635 TO ACCEPT GRIB MESSAGES C WITH SECOND ORDER PACKING. ADDED ROUTINE FI636 C TO PROCESS MESSAGES WITH SECOND ORDER PACKING. C 93-09-22 CAVANAUGH MODIFIED TO EXTRACT SUB-CENTER NUMBER FROM C PDS BYTE 26 C 93-10-13 CAVANAUGH MODIFIED FI634 TO CORRECT GRID SIZES FOR C GRIDS 204 AND 208 C 93-10-14 CAVANAUGH INCREASED SIZE OF KGDS TO INCLUDE ENTRIES FOR C NUMBER OF POINTS IN GRID AND NUMBER OF WORDS C IN EACH ROW C 93-12-08 CAVANAUGH CORRECTED TEST FOR EDITION NUMBER INSTEAD C OF VERSION NUMBER C 93-12-15 CAVANAUGH MODIFIED SECOND ORDER POINTERS TO FIRST ORDER C VALUES AND SECOND ORDER VALUES CORRECTLY C IN ROUTINE FI636 C 94-03-02 CAVANAUGH ADDED CALL TO W3FI83 WITHIN DECODER. USER C NO LONGER NEEDS TO MAKE CALL TO THIS ROUTINE C 94-04-22 CAVANAUGH MODIFIED FI635, FI636 TO PROCESS ROW BY ROW C SECOND ORDER PACKING, ADDED SCALING CORRECTION C TO FI635, AND CORRECTED TYPOGRAPHICAL ERRORS C IN COMMENT FIELDS IN FI634 C 94-05-17 CAVANAUGH CORRECTED ERROR IN FI633 TO EXTRACT RESOLUTION C FOR LAMBERT-CONFORMAL GRIDS. ADDED CLARIFYING C INFORMATION TO DOCBLOCK ENTRIES C 94-05-25 CAVANAUGH ADDED CODE TO PROCESS COLUMN BY COLUMN AS WELL C AS ROW BY ROW ORDERING OF SECOND ORDER DATA C 94-06-27 CAVANAUGH ADDED PROCESSING FOR GRIDS 45, 94 AND 95. C INCLUDES CONSTRUCTION OF SECOND ORDER BIT MAPS C FOR THINNED GRIDS IN FI636. C 94-07-08 CAVANAUGH COMMENTED OUT PRINT OUTS USED FOR DEBUGGING C 94-09-08 CAVANAUGH ADDED GRIDS 220, 221, 223 FOR FNOC C 94-11-10 FARLEY INCREASED MXSIZE FROM 72960 TO 260000 C FOR .5 DEGREE SST ANALYSIS FIELDS C 94-12-06 R.E.JONES CHANGES IN FI632 FOR PDS GREATER THAN 28 C 95-02-14 R.E.JONES CORRECT IN FI633 FOR NAVY WAFS GRIB C 95-03-20 M.BALDWIN FI633 QUICK AN DIRTY FIX MODIFICATION TO GET C DATA REP TYPE [KGDS(1)] 201 AND 202 TO WORK. C 95-04-10 E.ROGERS ADDED GRIDS 96 AND 97 FOR ETA MODEL IN FI634. C 95-04-26 R.E.JONES FI636 CORECTION FOR 2ND ORDER COMPLEX C UNPACKING. R C 95-05-19 R.E.JONES ADDED GRID 215, 20 KM AWIPS GRID C 95-07-06 R.E.JONES ADDED GAUSSIAN T62, T126 GRID 98, 126 C 95-10-19 R.E.JONES ADDED GRID 216, 45 KM ETA AWIPS ALASKA GRID C 95-10-31 IREDELL REMOVED SAVES AND PRINTS C 96-03-07 R.E.JONES CONTINUE UNPACK WITH KRET ERROR 9 IN FI631. C 96-08-19 R.E.JONES ADDED MERCATOR GRIDS 8 AND 53, AND GRID 196 C C USAGE: CALL W3FI63(MSGA,KPDS,KGDS,KBMS,DATA,KPTR,KRET) C INPUT ARGUMENT LIST: C MSGA - GRIB FIELD - "GRIB" THRU "7777" CHAR*1 C (MESSAGE CAN BE PRECEDED BY JUNK CHARS) C C OUTPUT ARGUMENT LIST: C DATA - ARRAY CONTAINING DATA ELEMENTS C KPDS - ARRAY CONTAINING PDS ELEMENTS. (EDITION 1) C (1) - ID OF CENTER C (2) - GENERATING PROCESS ID NUMBER C (3) - GRID DEFINITION C (4) - GDS/BMS FLAG (RIGHT ADJ COPY OF OCTET 8) C (5) - INDICATOR OF PARAMETER C (6) - TYPE OF LEVEL C (7) - HEIGHT/PRESSURE , ETC OF LEVEL C (8) - YEAR INCLUDING (CENTURY-1) C (9) - MONTH OF YEAR C (10) - DAY OF MONTH C (11) - HOUR OF DAY C (12) - MINUTE OF HOUR C (13) - INDICATOR OF FORECAST TIME UNIT C (14) - TIME RANGE 1 C (15) - TIME RANGE 2 C (16) - TIME RANGE FLAG C (17) - NUMBER INCLUDED IN AVERAGE C (18) - VERSION NR OF GRIB SPECIFICATION C (19) - VERSION NR OF PARAMETER TABLE C (20) - NR MISSING FROM AVERAGE/ACCUMULATION C (21) - CENTURY OF REFERENCE TIME OF DATA C (22) - UNITS DECIMAL SCALE FACTOR C (23) - SUBCENTER NUMBER C (24) - PDS BYTE 29, FOR NMC ENSEMBLE PRODUCTS C 128 IF FORECAST FIELD ERROR C 64 IF BIAS CORRECTED FCST FIELD C 32 IF SMOOTHED FIELD C WARNING: CAN BE COMBINATION OF MORE THAN 1 C (25) - PDS BYTE 30, NOT USED C (26-35) - RESERVED C (36-N) - CONSECUTIVE BYTES EXTRACTED FROM PROGRAM C DEFINITION SECTION (PDS) OF GRIB MESSAGE C KGDS - ARRAY CONTAINING GDS ELEMENTS. C (1) - DATA REPRESENTATION TYPE C (19) - NUMBER OF VERTICAL COORDINATE PARAMETERS C (20) - OCTET NUMBER OF THE LIST OF VERTICAL COORDINATE C PARAMETERS C OR C OCTET NUMBER OF THE LIST OF NUMBERS OF POINTS C IN EACH ROW C OR C 255 IF NEITHER ARE PRESENT C (21) - FOR GRIDS WITH PL, NUMBER OF POINTS IN GRID C (22) - NUMBER OF WORDS IN EACH ROW C LATITUDE/LONGITUDE GRIDS C (2) - N(I) NR POINTS ON LATITUDE CIRCLE C (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN C (4) - LA(1) LATITUDE OF ORIGIN C (5) - LO(1) LONGITUDE OF ORIGIN C (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17) C (7) - LA(2) LATITUDE OF EXTREME POINT C (8) - LO(2) LONGITUDE OF EXTREME POINT C (9) - DI LATITUDINAL DIRECTION OF INCREMENT C (10) - DJ LONGITUDINAL DIRECTION INCREMENT C (11) - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28) C GAUSSIAN GRIDS C (2) - N(I) NR POINTS ON LATITUDE CIRCLE C (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN C (4) - LA(1) LATITUDE OF ORIGIN C (5) - LO(1) LONGITUDE OF ORIGIN C (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17) C (7) - LA(2) LATITUDE OF EXTREME POINT C (8) - LO(2) LONGITUDE OF EXTREME POINT C (9) - DI LATITUDINAL DIRECTION OF INCREMENT C (10) - N - NR OF CIRCLES POLE TO EQUATOR C (11) - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28) C (12) - NV - NR OF VERT COORD PARAMETERS C (13) - PV - OCTET NR OF LIST OF VERT COORD PARAMETERS C OR C PL - LOCATION OF THE LIST OF NUMBERS OF POINTS IN C EACH ROW (IF NO VERT COORD PARAMETERS C ARE PRESENT C OR C 255 IF NEITHER ARE PRESENT C POLAR STEREOGRAPHIC GRIDS C (2) - N(I) NR POINTS ALONG LAT CIRCLE C (3) - N(J) NR POINTS ALONG LON CIRCLE C (4) - LA(1) LATITUDE OF ORIGIN C (5) - LO(1) LONGITUDE OF ORIGIN C (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17) C (7) - LOV GRID ORIENTATION C (8) - DX - X DIRECTION INCREMENT C (9) - DY - Y DIRECTION INCREMENT C (10) - PROJECTION CENTER FLAG C (11) - SCANNING MODE (RIGHT ADJ COPY OF OCTET 28) C SPHERICAL HARMONIC COEFFICIENTS C (2) - J PENTAGONAL RESOLUTION PARAMETER C (3) - K " " " C (4) - M " " " C (5) - REPRESENTATION TYPE C (6) - COEFFICIENT STORAGE MODE C MERCATOR GRIDS C (2) - N(I) NR POINTS ON LATITUDE CIRCLE C (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN C (4) - LA(1) LATITUDE OF ORIGIN C (5) - LO(1) LONGITUDE OF ORIGIN C (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17) C (7) - LA(2) LATITUDE OF LAST GRID POINT C (8) - LO(2) LONGITUDE OF LAST GRID POINT C (9) - LATIT - LATITUDE OF PROJECTION INTERSECTION C (10) - RESERVED C (11) - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28) C (12) - LONGITUDINAL DIR GRID LENGTH C (13) - LATITUDINAL DIR GRID LENGTH C LAMBERT CONFORMAL GRIDS C (2) - NX NR POINTS ALONG X-AXIS C (3) - NY NR POINTS ALONG Y-AXIS C (4) - LA1 LAT OF ORIGIN (LOWER LEFT) C (5) - LO1 LON OF ORIGIN (LOWER LEFT) C (6) - RESOLUTION (RIGHT ADJ COPY OF OCTET 17) C (7) - LOV - ORIENTATION OF GRID C (8) - DX - X-DIR INCREMENT C (9) - DY - Y-DIR INCREMENT C (10) - PROJECTION CENTER FLAG C (11) - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28) C (12) - LATIN 1 - FIRST LAT FROM POLE OF SECANT CONE INTER C (13) - LATIN 2 - SECOND LAT FROM POLE OF SECANT CONE INTER C KBMS - BITMAP DESCRIBING LOCATION OF OUTPUT ELEMENTS. C (ALWAYS CONSTRUCTED) C KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS C (1) - TOTAL LENGTH OF GRIB MESSAGE C (2) - LENGTH OF INDICATOR (SECTION 0) C (3) - LENGTH OF PDS (SECTION 1) C (4) - LENGTH OF GDS (SECTION 2) C (5) - LENGTH OF BMS (SECTION 3) C (6) - LENGTH OF BDS (SECTION 4) C (7) - VALUE OF CURRENT BYTE C (8) - BIT POINTER C (9) - GRIB START BIT NR C (10) - GRIB/GRID ELEMENT COUNT C (11) - NR UNUSED BITS AT END OF SECTION 3 C (12) - BIT MAP FLAG (COPY OF BMS OCTETS 5,6) C (13) - NR UNUSED BITS AT END OF SECTION 2 C (14) - BDS FLAGS (RIGHT ADJ COPY OF OCTET 4) C (15) - NR UNUSED BITS AT END OF SECTION 4 C KRET - FLAG INDICATING QUALITY OF COMPLETION C C REMARKS: WHEN DECODING IS COMPLETED, DATA AT EACH GRID POINT HAS BEEN C RETURNED IN THE UNITS SPECIFIED IN THE GRIB MANUAL. C C VALUES FOR RETURN FLAG (KRET) C KRET = 0 - NORMAL RETURN, NO ERRORS C = 1 - 'GRIB' NOT FOUND IN FIRST 100 CHARS C = 2 - '7777' NOT IN CORRECT LOCATION C = 3 - UNPACKED FIELD IS LARGER THAN 260000 C = 4 - GDS/ GRID NOT ONE OF CURRENTLY ACCEPTED VALUES C = 5 - GRID NOT CURRENTLY AVAIL FOR CENTER INDICATED C = 8 - TEMP GDS INDICATED, BUT GDS FLAG IS OFF C = 9 - GDS INDICATES SIZE MISMATCH WITH STD GRID C =10 - INCORRECT CENTER INDICATOR C =11 - BINARY DATA SECTION (BDS) NOT COMPLETELY PROCESSED. C PROGRAM IS NOT SET TO PROCESS FLAG COMBINATIONS C SHOWN IN OCTETS 4 AND 14. C =12 - BINARY DATA SECTION (BDS) NOT COMPLETELY PROCESSED. C PROGRAM IS NOT SET TO PROCESS FLAG COMBINATIONS C C SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT. C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C MACHINE: CRAY, WORKSTATIONS C C$$$ C 4 AUG 1988 C W3FI63 C C C GRIB UNPACKING ROUTINE C C C THIS ROUTINE WILL UNPACK A 'GRIB' FIELD TO THE EXACT GRID C TYPE SPECIFIED IN THE MESSAGE, RETURN A BIT MAP AND MAKE THE C VALUES OF THE PRODUCT DEFINITION SEC (PDS) AND THE GRID C DESCRIPTION SEC (GDS) AVAILABLE IN RETURN ARRAYS. C SEE "GRIB - THE WMO FORMAT FOR THE STORAGE OF WEATHER PRODUCT C INFORMATION AND THE EXCHANGE OF WEATHER PRODUCT MESSAGES IN C GRIDDED BINARY FORM" DATED JULY 1, 1988 BY JOHN D. STACKPOLE C DOC, NOAA, NWS, NATIONAL METEOROLOGICAL CENTER. C C THE CALL TO THE GRIB UNPACKING ROUTINE IS AS FOLLOWS: C C CALL W3FI63(MSGA,KPDS,KGDS,LBMS,DATA,KPTR,KRET) C C INPUT: C C MSGA = CONTAINS THE GRIB MESSAGE TO BE UNPACKED. CHARACTERS C "GRIB" MAY BEGIN ANYWHERE WITHIN FIRST 100 BYTES. C C OUTPUT: C C KPDS(100) INTEGER*4 C ARRAY TO CONTAIN THE ELEMENTS OF THE PRODUCT C DEFINITION SEC . C (VERSION 1) C KPDS(1) - ID OF CENTER C KPDS(2) - MODEL IDENTIFICATION (SEE "GRIB" TABLE 1) C KPDS(3) - GRID IDENTIFICATION (SEE "GRIB" TABLE 2) C KPDS(4) - GDS/BMS FLAG C BIT DEFINITION C 25 0 - GDS OMITTED C 1 - GDS INCLUDED C 26 0 - BMS OMITTED C 1 - BMS INCLUDED C NOTE:- LEFTMOST BIT = 1, C RIGHTMOST BIT = 32 C KPDS(5) - INDICATOR OF PARAMETER (SEE "GRIB" TABLE 5) C KPDS(6) - TYPE OF LEVEL (SEE "GRIB" TABLES 6 & 7) C KPDS(7) - HEIGHT,PRESSURE,ETC OF LEVEL C KPDS(8) - YEAR INCLUDING CENTURY C KPDS(9) - MONTH OF YEAR C KPDS(10) - DAY OF MONTH C KPDS(11) - HOUR OF DAY C KPDS(12) - MINUTE OF HOUR C KPDS(13) - INDICATOR OF FORECAST TIME UNIT (SEE "GRIB" C TABLE 8) C KPDS(14) - TIME 1 (SEE "GRIB" TABLE 8A) C KPDS(15) - TIME 2 (SEE "GRIB" TABLE 8A) C KPDS(16) - TIME RANGE INDICATOR (SEE "GRIB" TABLE 8A) C KPDS(17) - NUMBER INCLUDED IN AVERAGE C KPDS(18) - EDITION NR OF GRIB SPECIFICATION C KPDS(19) - VERSION NR OF PARAMETER TABLE C C KGDS(13) INTEGER*4 C ARRAY CONTAINING GDS ELEMENTS. C C KGDS(1) - DATA REPRESENTATION TYPE C C LATITUDE/LONGITUDE GRIDS (SEE "GRIB" TABLE 10) C KGDS(2) - N(I) NUMBER OF POINTS ON LATITUDE C CIRCLE C KGDS(3) - N(J) NUMBER OF POINTS ON LONGITUDE C CIRCLE C KGDS(4) - LA(1) LATITUDE OF ORIGIN C KGDS(5) - LO(1) LONGITUDE OF ORIGIN C KGDS(6) - RESOLUTION FLAG C BIT MEANING C 25 0 - DIRECTION INCREMENTS NOT C GIVEN C 1 - DIRECTION INCREMENTS GIVEN C KGDS(7) - LA(2) LATITUDE OF EXTREME POINT C KGDS(8) - LO(2) LONGITUDE OF EXTREME POINT C KGDS(9) - DI LATITUDINAL DIRECTION INCREMENT C KGDS(10) - REGULAR LAT/LON GRID C DJ - LONGITUDINAL DIRECTION C INCREMENT C GAUSSIAN GRID C N - NUMBER OF LATITUDE CIRCLES C BETWEEN A POLE AND THE EQUATOR C KGDS(11) - SCANNING MODE FLAG C BIT MEANING C 25 0 - POINTS ALONG A LATITUDE C SCAN FROM WEST TO EAST C 1 - POINTS ALONG A LATITUDE C SCAN FROM EAST TO WEST C 26 0 - POINTS ALONG A MERIDIAN C SCAN FROM NORTH TO SOUTH C 1 - POINTS ALONG A MERIDIAN C SCAN FROM SOUTH TO NORTH C 27 0 - POINTS SCAN FIRST ALONG C CIRCLES OF LATITUDE, THEN C ALONG MERIDIANS C (FORTRAN: (I,J)) C 1 - POINTS SCAN FIRST ALONG C MERIDIANS THEN ALONG C CIRCLES OF LATITUDE C (FORTRAN: (J,I)) C C POLAR STEREOGRAPHIC GRIDS (SEE GRIB TABLE 12) C KGDS(2) - N(I) NR POINTS ALONG LAT CIRCLE C KGDS(3) - N(J) NR POINTS ALONG LON CIRCLE C KGDS(4) - LA(1) LATITUDE OF ORIGIN C KGDS(5) - LO(1) LONGITUDE OF ORIGIN C KGDS(6) - RESERVED C KGDS(7) - LOV GRID ORIENTATION C KGDS(8) - DX - X DIRECTION INCREMENT C KGDS(9) - DY - Y DIRECTION INCREMENT C KGDS(10) - PROJECTION CENTER FLAG C KGDS(11) - SCANNING MODE C C SPHERICAL HARMONIC COEFFICIENTS (SEE "GRIB" TABLE 14) C KGDS(2) - J PENTAGONAL RESOLUTION PARAMETER C KGDS(3) - K PENTAGONAL RESOLUTION PARAMETER C KGDS(4) - M PENTAGONAL RESOLUTION PARAMETER C KGDS(5) - REPRESENTATION TYPE C KGDS(6) - COEFFICIENT STORAGE MODE C C MERCATOR GRIDS C KGDS(2) - N(I) NR POINTS ON LATITUDE CIRCLE C KGDS(3) - N(J) NR POINTS ON LONGITUDE MERIDIAN C KGDS(4) - LA(1) LATITUDE OF ORIGIN C KGDS(5) - LO(1) LONGITUDE OF ORIGIN C KGDS(6) - RESOLUTION FLAG C KGDS(7) - LA(2) LATITUDE OF LAST GRID POINT C KGDS(8) - LO(2) LONGITUDE OF LAST GRID POINT C KGDS(9) - LATIN - LATITUDE OF PROJECTION INTERSECTION C KGDS(10) - RESERVED C KGDS(11) - SCANNING MODE FLAG C KGDS(12) - LONGITUDINAL DIR GRID LENGTH C KGDS(13) - LATITUDINAL DIR GRID LENGTH C LAMBERT CONFORMAL GRIDS C KGDS(2) - NX NR POINTS ALONG X-AXIS C KGDS(3) - NY NR POINTS ALONG Y-AXIS C KGDS(4) - LA1 LAT OF ORIGIN (LOWER LEFT) C KGDS(5) - LO1 LON OF ORIGIN (LOWER LEFT) C KGDS(6) - RESOLUTION (RIGHT ADJ COPY OF OCTET 17) C KGDS(7) - LOV - ORIENTATION OF GRID C KGDS(8) - DX - X-DIR INCREMENT C KGDS(9) - DY - Y-DIR INCREMENT C KGDS(10) - PROJECTION CENTER FLAG C KGDS(11) - SCANNING MODE FLAG C KGDS(12) - LATIN 1 - FIRST LAT FROM POLE OF C SECANT CONE INTERSECTION C KGDS(13) - LATIN 2 - SECOND LAT FROM POLE OF C SECANT CONE INTERSECTION C C LBMS(260000) LOGICAL C ARRAY TO CONTAIN THE BIT MAP DESCRIBING THE C PLACEMENT OF DATA IN THE OUTPUT ARRAY. IF A C BIT MAP IS NOT INCLUDED IN THE SOURCE MESSAGE, C ONE WILL BE GENERATED AUTOMATICALLY BY THE C UNPACKING ROUTINE. C C C DATA(260000) REAL*4 C THIS ARRAY WILL CONTAIN THE UNPACKED DATA POINTS. C C NOTE:- 65160 IS MAXIMUN FIELD SIZE ALLOWABLE C C KPTR(10) INTEGER*4 C ARRAY CONTAINING STORAGE FOR THE FOLLOWING C PARAMETERS. C C (1) - UNUSED C (2) - UNUSED C (3) - LENGTH OF PDS (IN BYTES) C (4) - LENGTH OF GDS (IN BYTES) C (5) - LENGTH OF BMS (IN BYTES) C (6) - LENGTH OF BDS (IN BYTES) C (7) - USED BY UNPACKING ROUTINE C (8) - NUMBER OF DATA POINTS FOR GRID C (9) - "GRIB" CHARACTERS START IN BYTE NUMBER C (10) - USED BY UNPACKING ROUTINE C C C KRET INTEGER*4 C THIS VARIABLE WILL CONTAIN THE RETURN INDICATOR. C C 0 - NO ERRORS DETECTED. C C 1 - 'GRIB' NOT FOUND IN FIRST 100 C CHARACTERS. C C 2 - '7777' NOT FOUND, EITHER MISSING OR C TOTAL OF SEC COUNTS OF INDIVIDUAL C SECTIONS IS INCORRECT. C C 3 - UNPACKED FIELD IS LARGER THAN 65160. C C 4 - IN GDS, DATA REPRESENTATION TYPE C NOT ONE OF THE CURRENTLY ACCEPTABLE C VALUES. SEE "GRIB" TABLE 9. VALUE C OF INCORRECT TYPE RETURNED IN KGDS(1). C C 5 - GRID INDICATED IN KPDS(3) IS NOT C AVAILABLE FOR THE CENTER INDICATED IN C KPDS(1) AND NO GDS SENT. C C 7 - EDITION INDICATED IN KPDS(18) HAS NOT C YET BEEN INCLUDED IN THE DECODER. C C 8 - GRID IDENTIFICATION = 255 (NOT STANDARD C GRID) BUT FLAG INDICATING PRESENCE OF C GDS IS TURNED OFF. NO METHOD OF C GENERATING PROPER GRID. C C 9 - PRODUCT OF KGDS(2) AND KGDS(3) DOES NOT C MATCH STANDARD NUMBER OF POINTS FOR THIS C GRID (FOR OTHER THAN SPECTRALS). THIS C WILL OCCUR ONLY IF THE GRID. C IDENTIFICATION, KPDS(3), AND A C TRANSMITTED GDS ARE INCONSISTENT. C C 10 - CENTER INDICATOR WAS NOT ONE INDICATED C IN "GRIB" TABLE 1. PLEASE CONTACT AD C PRODUCTION MANAGEMENT BRANCH (W/NMC42) C IF THIS ERROR IS ENCOUNTERED. C C 11 - BINARY DATA SECTION (BDS) NOT COMPLETELY C PROCESSED. PROGRAM IS NOT SET TO PROCESS C FLAG COMBINATIONS AS SHOWN IN C OCTETS 4 AND 14. C C C LIST OF TEXT MESSAGES FROM CODE C C C W3FI63/FI632 C C 'HAVE ENCOUNTERED A NEW GRID FOR NMC, PLEASE NOTIFY C AUTOMATION DIVISION, PRODUCTION MANAGEMENT BRANCH C (W/NMC42)' C C 'HAVE ENCOUNTERED A NEW GRID FOR ECMWF, PLEASE NOTIFY C AUTOMATION DIVISION, PRODUCTION MANAGEMENT BRANCH C (W/NMC42)' C C 'HAVE ENCOUNTERED A NEW GRID FOR U.K. METEOROLOGICAL C OFFICE, BRACKNELL. PLEASE NOTIFY AUTOMATION DIVISION, C PRODUCTION MANAGEMENT BRANCH (W/NMC42)' C C 'HAVE ENCOUNTERED A NEW GRID FOR FNOC, PLEASE NOTIFY C AUTOMATION DIVISION, PRODUCTION MANAGEMENT BRANCH C (W/NMC42)' C C C W3FI63/FI633 C C 'POLAR STEREO PROCESSING NOT AVAILABLE' * C C W3FI63/FI634 C C 'WARNING - BIT MAP MAY NOT BE ASSOCIATED WITH SPHERICAL C COEFFICIENTS' C C C W3FI63/FI637 C C 'NO CURRENT LISTING OF FNOC GRIDS' * C C C * WILL BE AVAILABLE IN NEXT UPDATE C *************************************************************** C C INCOMING MESSAGE HOLDER CHARACTER*1 MSGA(*) C BIT MAP LOGICAL*1 KBMS(*) C C ELEMENTS OF PRODUCT DESCRIPTION SEC (PDS) INTEGER KPDS(*) C ELEMENTS OF GRID DESCRIPTION SEC (PDS) INTEGER KGDS(*) C C CONTAINER FOR GRIB GRID REAL DATA(*) C C ARRAY OF POINTERS AND COUNTERS INTEGER KPTR(*) C C ***************************************************************** INTEGER KKK,JSGN,JEXP,IFR,NPTS CHARACTER KK(8) REAL REALKK,FVAL1,FDIFF1 EQUIVALENCE (KK(1),KKK) C ***************************************************************** C 1.0 LOCATE BEGINNING OF 'GRIB' MESSAGE C FIND 'GRIB' CHARACTERS C 2.0 USE COUNTS IN EACH DESCRIPTION SEC TO DETERMINE C IF '7777' IS IN PROPER PLACE. C 3.0 PARSE PRODUCT DEFINITION SECTION. C 4.0 PARSE GRID DESCRIPTION SEC (IF INCLUDED) C 5.0 PARSE BIT MAP SEC (IF INCLUDED) C 6.0 USING INFORMATION FROM PRODUCT DEFINITION, GRID C DESCRIPTION, AND BIT MAP SECTIONS.. EXTRACT C DATA AND PLACE INTO PROPER ARRAY. C ******************************************************************* C C MAIN DRIVER C C ******************************************************************* KPTR(10) = 0 C SEE IF PROPER 'GRIB' KEY EXISTS, THEN C USING SEC COUNTS, DETERMINE IF '7777' C IS IN THE PROPER LOCATION C CALL FI631(MSGA,KPTR,KPDS,KRET) IF(KRET.NE.0) THEN GO TO 900 END IF C PRINT *,'FI631 KPTR',(KPTR(I),I=1,16) C C PARSE PARAMETERS FROM PRODUCT DESCRIPTION SECTION C CALL FI632(MSGA,KPTR,KPDS,KRET) IF(KRET.NE.0) THEN GO TO 900 END IF C PRINT *,'FI632 KPTR',(KPTR(I),I=1,16) C C IF AVAILABLE, EXTRACT NEW GRID DESCRIPTION C IF (IAND(KPDS(4),128).NE.0) THEN CALL FI633(MSGA,KPTR,KGDS,KRET) IF(KRET.NE.0) THEN GO TO 900 END IF C PRINT *,'FI633 KPTR',(KPTR(I),I=1,16) END IF C C EXTRACT OR GENERATE BIT MAP C CALL FI634(MSGA,KPTR,KPDS,KGDS,KBMS,KRET) IF (KRET.NE.0) THEN IF (KRET.NE.9) THEN GO TO 900 END IF END IF C PRINT *,'FI634 KPTR',(KPTR(I),I=1,16) C C USING INFORMATION FROM PDS, BMS AND BIT DATA SEC , C EXTRACT AND SAVE IN GRIB GRID, ALL DATA ENTRIES. C IF (KPDS(18).EQ.1) THEN CALL FI635(MSGA,KPTR,KPDS,KGDS,KBMS,DATA,KRET) IF (KPTR(3).EQ.50) THEN C C PDS EQUAL 50 BYTES C THEREFORE SOMETHING SPECIAL IS GOING ON C C IN THIS CASE 2ND DIFFERENCE PACKING C NEEDS TO BE UNDONE. C C EXTRACT FIRST VALUE FROM BYTE 41-44 PDS C KPTR(9) CONTAINS OFFSET TO START OF C GRIB MESSAGE. C EXTRACT FIRST FIRST-DIFFERENCE FROM BYTES 45-48 PDS C C AND EXTRACT SCALE FACTOR (E) TO UNDO 2**E C THAT WAS APPLIED PRIOR TO 2ND ORDER PACKING C AND PLACED IN PDS BYTES 49-51 C FACTOR IS A SIGNED TWO BYTE INTEGER C C ALSO NEED THE DECIMAL SCALING FROM PDS(27-28) C (AVAILABLE IN KPDS(22) FROM UNPACKER) C TO UNDO THE DECIMAL SCALING APPLIED TO THE C SECOND DIFFERENCES DURING UNPACKING. C SECOND DIFFS ALWAYS PACKED WITH 0 DECIMAL SCALE C BUT UNPACKER DOESNT KNOW THAT. C C CALL GBYTE (MSGA,FVAL1,KPTR(9)+384,32) C C NOTE INTEGERS, CHARACTERS AND EQUIVALENCES C DEFINED ABOVE TO MAKE THIS KKK EXTRACTION C WORK AND LINE UP ON WORD BOUNDARIES C CALL GBYTE (MSGA,KKK,KPTR(9)+384,32) C C THE NEXT CODE WILL CONVERT THE IBM370 FOATING POINT C TO THE FLOATING POINT USED ON YOUR MACHINE. C C 1ST TEST TO SEE IN ON 32 OR 64 BIT WORD MACHINE C LW = 4 OR 8; IF 8 MAY BE A CRAY C CALL W3FI01(LW) IF (LW.EQ.4) THEN CALL GBYTE (KK,JSGN,0,1) CALL GBYTE (KK,JEXP,1,7) CALL GBYTE (KK,IFR,8,24) ELSE CALL GBYTE (KK,JSGN,32,1) CALL GBYTE (KK,JEXP,33,7) CALL GBYTE (KK,IFR,40,24) ENDIF C IF (IFR.EQ.0) THEN REALKK = 0.0 ELSE IF (JEXP.EQ.0.AND.IFR.EQ.0) THEN REALKK = 0.0 ELSE REALKK = FLOAT(IFR) * 16.0 ** (JEXP - 64 - 6) IF (JSGN.NE.0) REALKK = -REALKK END IF FVAL1 = REALKK C C CALL GBYTE (MSGA,FDIFF1,KPTR(9)+416,32) C (REPLACED BY FOLLOWING EXTRACTION) C CALL GBYTE (MSGA,KKK,KPTR(9)+416,32) C C THE NEXT CODE WILL CONVERT THE IBM370 FOATING POINT C TO THE FLOATING POINT USED ON YOUR MACHINE. C C 1ST TEST TO SEE IN ON 32 OR 64 BIT WORD MACHINE C LW = 4 OR 8; IF 8 MAY BE A CRAY C CALL W3FI01(LW) IF (LW.EQ.4) THEN CALL GBYTE (KK,JSGN,0,1) CALL GBYTE (KK,JEXP,1,7) CALL GBYTE (KK,IFR,8,24) ELSE CALL GBYTE (KK,JSGN,32,1) CALL GBYTE (KK,JEXP,33,7) CALL GBYTE (KK,IFR,40,24) ENDIF C IF (IFR.EQ.0) THEN REALKK = 0.0 ELSE IF (JEXP.EQ.0.AND.IFR.EQ.0) THEN REALKK = 0.0 ELSE REALKK = FLOAT(IFR) * 16.0 ** (JEXP - 64 - 6) IF (JSGN.NE.0) REALKK = -REALKK END IF FDIFF1 = REALKK C CALL GBYTE (MSGA,ISIGN,KPTR(9)+448,1) CALL GBYTE (MSGA,ISCAL2,KPTR(9)+449,15) IF(ISIGN.GT.0) THEN ISCAL2 = - ISCAL2 ENDIF C PRINT *,'DELTA POINT 1-',FVAL1 C PRINT *,'DELTA POINT 2-',FDIFF1 C PRINT *,'DELTA POINT 3-',ISCAL2 NPTS = KPTR(10) C WRITE (6,FMT='('' 2ND DIFF POINTS IN FIELD = '',/, C & 10(3X,10F12.2,/))') (DATA(I),I=1,NPTS) C PRINT *,'DELTA POINT 4-',KPDS(22) CALL W3FI83 (DATA,NPTS,FVAL1,FDIFF1, & ISCAL2,KPDS(22),KPDS,KGDS) C WRITE (6,FMT='('' 2ND DIFF EXPANDED POINTS IN FIELD = '', C & /,10(3X,10F12.2,/))') (DATA(I),I=1,NPTS) C WRITE (6,FMT='('' END OF ARRAY IN FIELD = '',/, C & 10(3X,10F12.2,/))') (DATA(I),I=NPTS-5,NPTS) END IF ELSE C PRINT *,'FI635 NOT PROGRAMMED FOR EDITION NR',KPDS(18) KRET = 7 END IF C 900 RETURN END SUBROUTINE FI631(MSGA,KPTR,KPDS,KRET) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: FI631 FIND 'GRIB' CHARS & RESET POINTERS C PRGMMR: BILL CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 C C ABSTRACT: FIND 'GRIB; CHARACTERS AND SET POINTERS TO THE NEXT C BYTE FOLLOWING 'GRIB'. IF THEY EXIST EXTRACT COUNTS FROM GDS AND C BMS. EXTRACT COUNT FROM BDS. DETERMINE IF SUM OF COUNTS ACTUALLY C PLACES TERMINATOR '7777' AT THE CORRECT LOCATION. C C PROGRAM HISTORY LOG: C 91-09-13 CAVANAUGH C 95-10-31 IREDELL REMOVED SAVES AND PRINTS C C USAGE: CALL FI631(MSGA,KPTR,KPDS,KRET) C INPUT ARGUMENT LIST: C MSGA - GRIB FIELD - "GRIB" THRU "7777" C KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS C (1) - TOTAL LENGTH OF GRIB MESSAGE C (2) - LENGTH OF INDICATOR (SECTION 0) C (3) - LENGTH OF PDS (SECTION 1) C (4) - LENGTH OF GDS (SECTION 2) C (5) - LENGTH OF BMS (SECTION 3) C (6) - LENGTH OF BDS (SECTION 4) C (7) - VALUE OF CURRENT BYTE C (8) - BIT POINTER C (9) - GRIB START BIT NR C (10) - GRIB/GRID ELEMENT COUNT C (11) - NR UNUSED BITS AT END OF SECTION 3 C (12) - BIT MAP FLAG C (13) - NR UNUSED BITS AT END OF SECTION 2 C (14) - BDS FLAGS C (15) - NR UNUSED BITS AT END OF SECTION 4 C C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) C KPDS - ARRAY CONTAINING PDS ELEMENTS. C (1) - ID OF CENTER C (2) - MODEL IDENTIFICATION C (3) - GRID IDENTIFICATION C (4) - GDS/BMS FLAG C (5) - INDICATOR OF PARAMETER C (6) - TYPE OF LEVEL C (7) - HEIGHT/PRESSURE , ETC OF LEVEL C (8) - YEAR OF CENTURY C (9) - MONTH OF YEAR C (10) - DAY OF MONTH C (11) - HOUR OF DAY C (12) - MINUTE OF HOUR C (13) - INDICATOR OF FORECAST TIME UNIT C (14) - TIME RANGE 1 C (15) - TIME RANGE 2 C (16) - TIME RANGE FLAG C (17) - NUMBER INCLUDED IN AVERAGE C KPTR - SEE INPUT LIST C KRET - ERROR RETURN C C REMARKS: C ERROR RETURNS C KRET = 1 - NO 'GRIB' C 2 - NO '7777' OR MISLOCATED (BY COUNTS) C C SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT. C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C MACHINE: HDS9000 C C$$$ C C INCOMING MESSAGE HOLDER CHARACTER*1 MSGA(*) C ARRAY OF POINTERS AND COUNTERS INTEGER KPTR(*) C PRODUCT DESCRIPTION SECTION DATA. INTEGER KPDS(*) C INTEGER KRET C C ****************************************************************** KRET = 0 C ------------------- FIND 'GRIB' KEY DO 50 I = 0, 839, 8 CALL GBYTE (MSGA,MGRIB,I,32) IF (MGRIB.EQ.1196575042) THEN KPTR(9) = I GO TO 60 END IF 50 CONTINUE KRET = 1 RETURN 60 CONTINUE C -------------FOUND 'GRIB' C SKIP GRIB CHARACTERS C PRINT *,'FI631 GRIB AT',I KPTR(8) = KPTR(9) + 32 CALL GBYTE (MSGA,ITOTAL,KPTR(8),24) C HAVE LIFTED WHAT MAY BE A MSG TOTAL BYTE COUNT IPOINT = KPTR(9) + ITOTAL * 8 - 32 CALL GBYTE (MSGA,I7777,IPOINT,32) IF (I7777.EQ.926365495) THEN C HAVE FOUND END OF MESSAGE '7777' IN PROPER LOCATION C MARK AND PROCESS AS GRIB VERSION 1 OR HIGHER C PRINT *,'FI631 7777 AT',IPOINT KPTR(8) = KPTR(8) + 24 KPTR(1) = ITOTAL KPTR(2) = 8 CALL GBYTE (MSGA,KPDS(18),KPTR(8),8) KPTR(8) = KPTR(8) + 8 ELSE C CANNOT FIND END OF GRIB EDITION 1 MESSAGE KRET = 2 RETURN END IF C ------------------- PROCESS SECTION 1 C EXTRACT COUNT FROM PDS C PRINT *,'START OF PDS',KPTR(8) CALL GBYTE (MSGA,KPTR(3),KPTR(8),24) LOOK = KPTR(8) + 56 C EXTRACT GDS/BMS FLAG CALL GBYTE (MSGA,KPDS(4),LOOK,8) KPTR(8) = KPTR(8) + KPTR(3) * 8 C PRINT *,'START OF GDS',KPTR(8) IF (IAND(KPDS(4),128).NE.0) THEN C EXTRACT COUNT FROM GDS CALL GBYTE (MSGA,KPTR(4),KPTR(8),24) KPTR(8) = KPTR(8) + KPTR(4) * 8 ELSE KPTR(4) = 0 END IF C PRINT *,'START OF BMS',KPTR(8) IF (IAND(KPDS(4),64).NE.0) THEN C EXTRACT COUNT FROM BMS CALL GBYTE (MSGA,KPTR(5),KPTR(8),24) ELSE KPTR(5) = 0 END IF KPTR(8) = KPTR(8) + KPTR(5) * 8 C PRINT *,'START OF BDS',KPTR(8) C EXTRACT COUNT FROM BDS CALL GBYTE (MSGA,KPTR(6),KPTR(8),24) C --------------- TEST FOR '7777' C PRINT *,(KPTR(KJ),KJ=1,10) KPTR(8) = KPTR(8) + KPTR(6) * 8 C EXTRACT FOUR BYTES FROM THIS LOCATION C PRINT *,'FI631 LOOKING FOR 7777 AT',KPTR(8) CALL GBYTE (MSGA,K7777,KPTR(8),32) MATCH = KPTR(2) + KPTR(3) + KPTR(4) + KPTR(5) + KPTR(6) + 4 IF (K7777.NE.926365495.OR.MATCH.NE.KPTR(1)) THEN KRET = 2 ELSE C PRINT *,'FI631 7777 AT',KPTR(8) IF (KPDS(18).EQ.0) THEN KPTR(1) = KPTR(2) + KPTR(3) + KPTR(4) + KPTR(5) + * KPTR(6) + 4 END IF END IF C PRINT *,'KPTR',(KPTR(I),I=1,16) RETURN END SUBROUTINE FI632(MSGA,KPTR,KPDS,KRET) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: FI632 GATHER INFO FROM PRODUCT DEFINITION SEC C PRGMMR: BILL CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 C C ABSTRACT: EXTRACT INFORMATION FROM THE PRODUCT DESCRIPTION C SEC , AND GENERATE LABEL INFORMATION TO PERMIT STORAGE C IN OFFICE NOTE 84 FORMAT. C C PROGRAM HISTORY LOG: C 91-09-13 CAVANAUGH C 93-12-08 CAVANAUGH CORRECTED TEST FOR EDITION NUMBER INSTEAD C OF VERSION NUMBER C 95-10-31 IREDELL REMOVED SAVES AND PRINTS C C USAGE: CALL FI632(MSGA,KPTR,KPDS,KRET) C INPUT ARGUMENT LIST: C MSGA - ARRAY CONTAINING GRIB MESSAGE C KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS C (1) - TOTAL LENGTH OF GRIB MESSAGE C (2) - LENGTH OF INDICATOR (SECTION 0) C (3) - LENGTH OF PDS (SECTION 1) C (4) - LENGTH OF GDS (SECTION 2) C (5) - LENGTH OF BMS (SECTION 3) C (6) - LENGTH OF BDS (SECTION 4) C (7) - VALUE OF CURRENT BYTE C (8) - BIT POINTER C (9) - GRIB START BIT NR C (10) - GRIB/GRID ELEMENT COUNT C (11) - NR UNUSED BITS AT END OF SECTION 3 C (12) - BIT MAP FLAG C (13) - NR UNUSED BITS AT END OF SECTION 2 C (14) - BDS FLAGS C (15) - NR UNUSED BITS AT END OF SECTION 4 C C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) C KPDS - ARRAY CONTAINING PDS ELEMENTS. C (1) - ID OF CENTER C (2) - MODEL IDENTIFICATION C (3) - GRID IDENTIFICATION C (4) - GDS/BMS FLAG C (5) - INDICATOR OF PARAMETER C (6) - TYPE OF LEVEL C (7) - HEIGHT/PRESSURE , ETC OF LEVEL C (8) - YEAR OF CENTURY C (9) - MONTH OF YEAR C (10) - DAY OF MONTH C (11) - HOUR OF DAY C (12) - MINUTE OF HOUR C (13) - INDICATOR OF FORECAST TIME UNIT C (14) - TIME RANGE 1 C (15) - TIME RANGE 2 C (16) - TIME RANGE FLAG C (17) - NUMBER INCLUDED IN AVERAGE C (18) - C (19) - C (20) - NUMBER MISSING FROM AVGS/ACCUMULATIONS C (21) - CENTURY C (22) - UNITS DECIMAL SCALE FACTOR C (23) - SUBCENTER C KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS C SEE INPUT LIST C KRET - ERROR RETURN C C REMARKS: C ERROR RETURN = 0 - NO ERRORS C = 8 - TEMP GDS INDICATED, BUT NO GDS C C SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT. C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C MACHINE: HDS9000 C C$$$ C C INCOMING MESSAGE HOLDER CHARACTER*1 MSGA(*) C C ARRAY OF POINTERS AND COUNTERS INTEGER KPTR(*) C PRODUCT DESCRIPTION SECTION ENTRIES INTEGER KPDS(*) C INTEGER KRET C ------------------- PROCESS SECTION 1 KPTR(8) = KPTR(9) + KPTR(2) * 8 + 24 C BYTE 4 C PARAMETER TABLE VERSION NR CALL GBYTE (MSGA,KPDS(19),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 5 IDENTIFICATION OF CENTER CALL GBYTE (MSGA,KPDS(1),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 6 C GET GENERATING PROCESS ID NR CALL GBYTE (MSGA,KPDS(2),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 7 C GRID DEFINITION CALL GBYTE (MSGA,KPDS(3),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 8 C GDS/BMS FLAGS C CALL GBYTE (MSGA,KPDS(4),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 9 C INDICATOR OF PARAMETER CALL GBYTE (MSGA,KPDS(5),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 10 C TYPE OF LEVEL CALL GBYTE (MSGA,KPDS(6),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 11,12 C HEIGHT/PRESSURE CALL GBYTE (MSGA,KPDS(7),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C BYTE 13 C YEAR OF CENTURY CALL GBYTE (MSGA,KPDS(8),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 14 C MONTH OF YEAR CALL GBYTE (MSGA,KPDS(9),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 15 C DAY OF MONTH CALL GBYTE (MSGA,KPDS(10),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 16 C HOUR OF DAY CALL GBYTE (MSGA,KPDS(11),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 17 C MINUTE CALL GBYTE (MSGA,KPDS(12),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 18 C INDICATOR TIME UNIT RANGE CALL GBYTE (MSGA,KPDS(13),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 19 C P1 - PERIOD OF TIME CALL GBYTE (MSGA,KPDS(14),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 20 C P2 - PERIOD OF TIME CALL GBYTE (MSGA,KPDS(15),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 21 C TIME RANGE INDICATOR CALL GBYTE (MSGA,KPDS(16),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C C IF TIME RANGE INDICATOR IS 10, P1 IS PACKED IN C PDS BYTES 19-20 C IF (KPDS(16).EQ.10) THEN KPDS(14) = KPDS(14) * 256 + KPDS(15) KPDS(15) = 0 END IF C BYTE 22,23 C NUMBER INCLUDED IN AVERAGE CALL GBYTE (MSGA,KPDS(17),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C BYTE 24 C NUMBER MISSING FROM AVERAGES/ACCUMULATIONS CALL GBYTE (MSGA,KPDS(20),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 25 C IDENTIFICATION OF CENTURY CALL GBYTE (MSGA,KPDS(21),KPTR(8),8) KPTR(8) = KPTR(8) + 8 IF (KPTR(3).GT.25) THEN C BYTE 26 SUB CENTER NUMBER CALL GBYTE (MSGA,KPDS(23),KPTR(8),8) KPTR(8) = KPTR(8) + 8 IF (KPTR(3).GE.28) THEN C BYTE 27-28 C UNITS DECIMAL SCALE FACTOR CALL GBYTE (MSGA,ISIGN,KPTR(8),1) KPTR(8) = KPTR(8) + 1 CALL GBYTE (MSGA,IDEC,KPTR(8),15) KPTR(8) = KPTR(8) + 15 IF (ISIGN.GT.0) THEN KPDS(22) = - IDEC ELSE KPDS(22) = IDEC END IF ISIZ = KPTR(3) - 28 IF (ISIZ.LE.12) THEN C BYTE 29 CALL GBYTE (MSGA,KPDS(24),KPTR(8)+8,8) C BYTE 30 CALL GBYTE (MSGA,KPDS(25),KPTR(8)+16,8) C BYTES 31-40 CURRENTLY RESERVED FOR FUTURE USE KPTR(8) = KPTR(8) + ISIZ * 8 ELSE C BYTE 29 CALL GBYTE (MSGA,KPDS(24),KPTR(8)+8,8) C BYTE 30 CALL GBYTE (MSGA,KPDS(25),KPTR(8)+16,8) C BYTES 31-40 CURRENTLY RESERVED FOR FUTURE USE KPTR(8) = KPTR(8) + 12 * 8 C BYTES 41 - N LOCAL USE DATA CALL W3FI01(LW) MWDBIT = LW * 8 ISIZ = KPTR(3) - 40 ITER = ISIZ / LW IF (MOD(ISIZ,LW).NE.0) ITER = ITER + 1 CALL GBYTES (MSGA,KPDS(36),KPTR(8),MWDBIT,0,ITER) KPTR(8) = KPTR(8) + ISIZ * 8 END IF END IF END IF C ----------- TEST FOR NEW GRID IF (IAND(KPDS(4),128).NE.0) THEN IF (IAND(KPDS(4),64).NE.0) THEN IF (KPDS(3).NE.255) THEN IF (KPDS(3).GE.21.AND.KPDS(3).LE.26)THEN RETURN ELSE IF (KPDS(3).GE.37.AND.KPDS(3).LE.44)THEN RETURN ELSE IF (KPDS(3).GE.61.AND.KPDS(3).LE.64) THEN RETURN END IF IF (KPDS(1).EQ.7) THEN IF (KPDS(3).GE.2.AND.KPDS(3).LE.3) THEN ELSE IF (KPDS(3).GE.5.AND.KPDS(3).LE.6) THEN ELSE IF (KPDS(3).EQ.8) THEN ELSE IF (KPDS(3).GE.27.AND.KPDS(3).LE.34) THEN ELSE IF (KPDS(3).EQ.50) THEN ELSE IF (KPDS(3).EQ.53) THEN ELSE IF (KPDS(3).GE.70.AND.KPDS(3).LE.77) THEN ELSE IF (KPDS(3).EQ.98) THEN ELSE IF (KPDS(3).GE.100.AND.KPDS(3).LE.105) THEN ELSE IF (KPDS(3).EQ.126) THEN ELSE IF (KPDS(3).EQ.196) THEN ELSE IF (KPDS(3).GE.201.AND.KPDS(3).LE.216) THEN ELSE C PRINT *,' HAVE ENCOUNTERED A NEW GRID FOR', C * ' NMC WITHOUT A GRID DESCRIPTION SECTION' C PRINT *,' PLEASE NOTIFY AUTOMATION DIVISION' C PRINT *,' PRODUCTION MANAGEMENT BRANCH' C PRINT *,' W/NMC42)' END IF ELSE IF (KPDS(1).EQ.98) THEN IF (KPDS(3).GE.1.AND.KPDS(3).LE.16) THEN ELSE C PRINT *,' HAVE ENCOUNTERED A NEW GRID FOR', C * ' ECMWF WITHOUT A GRID DESCRIPTION SECTION' C PRINT *,' PLEASE NOTIFY AUTOMATION DIVISION' C PRINT *,' PRODUCTION MANAGEMENT BRANCH' C PRINT *,' W/NMC42)' END IF ELSE IF (KPDS(1).EQ.74) THEN IF (KPDS(3).GE.1.AND.KPDS(3).LE.12) THEN ELSE IF (KPDS(3).GE.21.AND.KPDS(3).LE.26)THEN ELSE IF (KPDS(3).GE.61.AND.KPDS(3).LE.64) THEN ELSE IF (KPDS(3).GE.70.AND.KPDS(3).LE.77) THEN ELSE C PRINT *,' HAVE ENCOUNTERED A NEW GRID FOR', C * ' U.K. MET OFFICE, BRACKNELL', C * ' WITHOUT A GRID DESCRIPTION SECTION' C PRINT *,' PLEASE NOTIFY AUTOMATION DIVISION' C PRINT *,' PRODUCTION MANAGEMENT BRANCH' C PRINT *,' W/NMC42)' END IF ELSE IF (KPDS(1).EQ.58) THEN IF (KPDS(3).GE.1.AND.KPDS(3).LE.12) THEN ELSE C PRINT *,' HAVE ENCOUNTERED A NEW GRID FOR', C * ' FNOC WITHOUT A GRID DESCRIPTION SECTION' C PRINT *,' PLEASE NOTIFY AUTOMATION DIVISION' C PRINT *,' PRODUCTION MANAGEMENT BRANCH' C PRINT *,' W/NMC42)' END IF END IF END IF END IF END IF RETURN END SUBROUTINE FI633(MSGA,KPTR,KGDS,KRET) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: FI633 EXTRACT INFO FROM GRIB-GDS C PRGMMR: BILL CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 C C ABSTRACT: EXTRACT INFORMATION ON UNLISTED GRID TO ALLOW C CONVERSION TO OFFICE NOTE 84 FORMAT. C C PROGRAM HISTORY LOG: C 91-09-13 CAVANAUGH C 95-03-20 M.BALDWIN FI633 QUICK AN DIRTY FIX MODIFICATION TO GET C DATA REP TYPE [KGDS(1)] 201 AND 202 TO WORK. C 95-10-31 IREDELL REMOVED SAVES AND PRINTS C C USAGE: CALL FI633(MSGA,KPTR,KGDS,KRET) C INPUT ARGUMENT LIST: C MSGA - ARRAY CONTAINING GRIB MESSAGE C KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS C (1) - TOTAL LENGTH OF GRIB MESSAGE C (2) - LENGTH OF INDICATOR (SECTION 0) C (3) - LENGTH OF PDS (SECTION 1) C (4) - LENGTH OF GDS (SECTION 2) C (5) - LENGTH OF BMS (SECTION 3) C (6) - LENGTH OF BDS (SECTION 4) C (7) - VALUE OF CURRENT BYTE C (8) - BIT POINTER C (9) - GRIB START BIT NR C (10) - GRIB/GRID ELEMENT COUNT C (11) - NR UNUSED BITS AT END OF SECTION 3 C (12) - BIT MAP FLAG C (13) - NR UNUSED BITS AT END OF SECTION 2 C (14) - BDS FLAGS C (15) - NR UNUSED BITS AT END OF SECTION 4 C C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) C KGDS - ARRAY CONTAINING GDS ELEMENTS. C (1) - DATA REPRESENTATION TYPE C (19) - NUMBER OF VERTICAL COORDINATE PARAMETERS C (20) - OCTET NUMBER OF THE LIST OF VERTICAL COORDINATE C PARAMETERS C OR C OCTET NUMBER OF THE LIST OF NUMBERS OF POINTS C IN EACH ROW C OR C 255 IF NEITHER ARE PRESENT C (21) - FOR GRIDS WITH PL, NUMBER OF POINTS IN GRID C (22) - NUMBER OF WORDS IN EACH ROW C LATITUDE/LONGITUDE GRIDS C (2) - N(I) NR POINTS ON LATITUDE CIRCLE C (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN C (4) - LA(1) LATITUDE OF ORIGIN C (5) - LO(1) LONGITUDE OF ORIGIN C (6) - RESOLUTION FLAG C (7) - LA(2) LATITUDE OF EXTREME POINT C (8) - LO(2) LONGITUDE OF EXTREME POINT C (9) - DI LATITUDINAL DIRECTION OF INCREMENT C (10) - DJ LONGITUDINAL DIRECTION INCREMENT C (11) - SCANNING MODE FLAG C POLAR STEREOGRAPHIC GRIDS C (2) - N(I) NR POINTS ALONG LAT CIRCLE C (3) - N(J) NR POINTS ALONG LON CIRCLE C (4) - LA(1) LATITUDE OF ORIGIN C (5) - LO(1) LONGITUDE OF ORIGIN C (6) - RESERVED C (7) - LOV GRID ORIENTATION C (8) - DX - X DIRECTION INCREMENT C (9) - DY - Y DIRECTION INCREMENT C (10) - PROJECTION CENTER FLAG C (11) - SCANNING MODE C SPHERICAL HARMONIC COEFFICIENTS C (2) - J PENTAGONAL RESOLUTION PARAMETER C (3) - K " " " C (4) - M " " " C (5) - REPRESENTATION TYPE C (6) - COEFFICIENT STORAGE MODE C MERCATOR GRIDS C (2) - N(I) NR POINTS ON LATITUDE CIRCLE C (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN C (4) - LA(1) LATITUDE OF ORIGIN C (5) - LO(1) LONGITUDE OF ORIGIN C (6) - RESOLUTION FLAG C (7) - LA(2) LATITUDE OF LAST GRID POINT C (8) - LO(2) LONGITUDE OF LAST GRID POINT C (9) - LATIN - LATITUDE OF PROJECTION INTERSECTION C (10) - RESERVED C (11) - SCANNING MODE FLAG C (12) - LONGITUDINAL DIR GRID LENGTH C (13) - LATITUDINAL DIR GRID LENGTH C LAMBERT CONFORMAL GRIDS C (2) - NX NR POINTS ALONG X-AXIS C (3) - NY NR POINTS ALONG Y-AXIS C (4) - LA1 LAT OF ORIGIN (LOWER LEFT) C (5) - LO1 LON OF ORIGIN (LOWER LEFT) C (6) - RESOLUTION (RIGHT ADJ COPY OF OCTET 17) C (7) - LOV - ORIENTATION OF GRID C (8) - DX - X-DIR INCREMENT C (9) - DY - Y-DIR INCREMENT C (10) - PROJECTION CENTER FLAG C (11) - SCANNING MODE FLAG C (12) - LATIN 1 - FIRST LAT FROM POLE OF SECANT CONE INTER C (13) - LATIN 2 - SECOND LAT FROM POLE OF SECANT CONE INTER C KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS C SEE INPUT LIST C KRET - ERROR RETURN C C REMARKS: C KRET = 0 C = 4 - DATA REPRESENTATION TYPE NOT CURRENTLY ACCEPTABLE C C SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT. C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C MACHINE: HDS9000 C C$$$ C ************************************************************ C INCOMING MESSAGE HOLDER CHARACTER*1 MSGA(*) C C ARRAY GDS ELEMENTS INTEGER KGDS(*) C ARRAY OF POINTERS AND COUNTERS INTEGER KPTR(*) C INTEGER KRET C --------------------------------------------------------------- KRET = 0 C PROCESS GRID DEFINITION SECTION (IF PRESENT) C MAKE SURE BIT POINTER IS PROPERLY SET KPTR(8) = KPTR(9) + (KPTR(2)*8) + (KPTR(3)*8) + 24 NSAVE = KPTR(8) - 24 C BYTE 4 C NV - NR OF VERT COORD PARAMETERS CALL GBYTE (MSGA,KGDS(19),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 5 C PV - LOCATION - SEE FM92 MANUAL CALL GBYTE (MSGA,KGDS(20),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTE 6 C DATA REPRESENTATION TYPE CALL GBYTE (MSGA,KGDS(1),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTES 7-32 ARE GRID DEFINITION DEPENDING ON C DATA REPRESENTATION TYPE IF (KGDS(1).EQ.0) THEN GO TO 1000 ELSE IF (KGDS(1).EQ.1) THEN GO TO 4000 ELSE IF (KGDS(1).EQ.2.OR.KGDS(1).EQ.5) THEN GO TO 2000 ELSE IF (KGDS(1).EQ.3) THEN GO TO 5000 ELSE IF (KGDS(1).EQ.4) THEN GO TO 1000 C ELSE IF (KGDS(1).EQ.10) THEN C ELSE IF (KGDS(1).EQ.14) THEN C ELSE IF (KGDS(1).EQ.20) THEN C ELSE IF (KGDS(1).EQ.24) THEN C ELSE IF (KGDS(1).EQ.30) THEN C ELSE IF (KGDS(1).EQ.34) THEN ELSE IF (KGDS(1).EQ.50) THEN GO TO 3000 C ELSE IF (KGDS(1).EQ.60) THEN C ELSE IF (KGDS(1).EQ.70) THEN C ELSE IF (KGDS(1).EQ.80) THEN ELSE IF (KGDS(1).EQ.201.OR.KGDS(1).EQ.202) THEN GO TO 1000 ELSE C MARK AS GDS/ UNKNOWN DATA REPRESENTATION TYPE KRET = 4 RETURN END IF C BYTE 33-N VERTICAL COORDINATE PARAMETERS C ----------- C BYTES 33-42 EXTENSIONS OF GRID DEFINITION FOR ROTATION C OR STRETCHING OF THE COORDINATE SYSTEM OR C LAMBERT CONFORMAL PROJECTION. C BYTE 43-N VERTICAL COORDINATE PARAMETERS C ----------- C BYTES 33-52 EXTENSIONS OF GRID DEFINITION FOR STRETCHED C AND ROTATED COORDINATE SYSTEM C BYTE 53-N VERTICAL COORDINATE PARAMETERS C ----------- C ************************************************************ C ------------------- LATITUDE/LONGITUDE GRIDS C C ------------------- BYTE 7-8 NR OF POINTS ALONG LATITUDE CIRCLE 1000 CONTINUE CALL GBYTE (MSGA,KGDS(2),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 9-10 NR OF POINTS ALONG LONG MERIDIAN CALL GBYTE (MSGA,KGDS(3),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 11-13 LATITUDE OF ORIGIN CALL GBYTE (MSGA,KGDS(4),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(4),8388608).NE.0) THEN KGDS(4) = IAND(KGDS(4),8388607) * (-1) END IF C ------------------- BYTE 14-16 LONGITUDE OF ORIGIN CALL GBYTE (MSGA,KGDS(5),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(5),8388608).NE.0) THEN KGDS(5) = - IAND(KGDS(5),8388607) END IF C ------------------- BYTE 17 RESOLUTION FLAG CALL GBYTE (MSGA,KGDS(6),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 18-20 LATITUDE OF LAST GRID POINT CALL GBYTE (MSGA,KGDS(7),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(7),8388608).NE.0) THEN KGDS(7) = - IAND(KGDS(7),8388607) END IF C ------------------- BYTE 21-23 LONGITUDE OF LAST GRID POINT CALL GBYTE (MSGA,KGDS(8),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(8),8388608).NE.0) THEN KGDS(8) = - IAND(KGDS(8),8388607) END IF C ------------------- BYTE 24-25 LATITUDINAL DIR INCREMENT CALL GBYTE (MSGA,KGDS(9),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 26-27 IF REGULAR LAT/LON GRID C HAVE LONGIT DIR INCREMENT C ELSE IF GAUSSIAN GRID C HAVE NR OF LAT CIRCLES C BETWEEN POLE AND EQUATOR CALL GBYTE (MSGA,KGDS(10),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 28 SCANNING MODE FLAGS CALL GBYTE (MSGA,KGDS(11),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 29-32 RESERVED C SKIP TO START OF BYTE 33 CALL GBYTE (MSGA,KGDS(12),KPTR(8),32) KPTR(8) = KPTR(8) + 32 C ------------------- GO TO 900 C ****************************************************************** C ' POLAR STEREO PROCESSING ' C C ------------------- BYTE 7-8 NR OF POINTS ALONG X=AXIS 2000 CONTINUE CALL GBYTE (MSGA,KGDS(2),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 9-10 NR OF POINTS ALONG Y-AXIS CALL GBYTE (MSGA,KGDS(3),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 11-13 LATITUDE OF ORIGIN CALL GBYTE (MSGA,KGDS(4),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(4),8388608).NE.0) THEN KGDS(4) = - IAND(KGDS(4),8388607) END IF C ------------------- BYTE 14-16 LONGITUDE OF ORIGIN CALL GBYTE (MSGA,KGDS(5),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(5),8388608).NE.0) THEN KGDS(5) = - IAND(KGDS(5),8388607) END IF C ------------------- BYTE 17 RESERVED CALL GBYTE (MSGA,KGDS(6),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 18-20 LOV ORIENTATION OF THE GRID CALL GBYTE (MSGA,KGDS(7),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(7),8388608).NE.0) THEN KGDS(7) = - IAND(KGDS(7),8388607) END IF C ------------------- BYTE 21-23 DX - THE X DIRECTION INCREMENT CALL GBYTE (MSGA,KGDS(8),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(8),8388608).NE.0) THEN KGDS(8) = - IAND(KGDS(8),8388607) END IF C ------------------- BYTE 24-26 DY - THE Y DIRECTION INCREMENT CALL GBYTE (MSGA,KGDS(9),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(9),8388608).NE.0) THEN KGDS(9) = - IAND(KGDS(9),8388607) END IF C ------------------- BYTE 27 PROJECTION CENTER FLAG CALL GBYTE (MSGA,KGDS(10),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 28 SCANNING MODE CALL GBYTE (MSGA,KGDS(11),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 29-32 RESERVED C SKIP TO START OF BYTE 33 CALL GBYTE (MSGA,KGDS(12),KPTR(8),32) KPTR(8) = KPTR(8) + 32 C C ------------------- GO TO 900 C C ****************************************************************** C ------------------- GRID DESCRIPTION FOR SPHERICAL HARMONIC COEFF. C C ------------------- BYTE 7-8 J PENTAGONAL RESOLUTION PARAMETER 3000 CONTINUE CALL GBYTE (MSGA,KGDS(2),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 9-10 K PENTAGONAL RESOLUTION PARAMETER CALL GBYTE (MSGA,KGDS(3),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 11-12 M PENTAGONAL RESOLUTION PARAMETER CALL GBYTE (MSGA,KGDS(4),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 13 REPRESENTATION TYPE CALL GBYTE (MSGA,KGDS(5),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 14 COEFFICIENT STORAGE MODE CALL GBYTE (MSGA,KGDS(6),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- EMPTY FIELDS - BYTES 15 - 32 C SET TO START OF BYTE 33 KPTR(8) = KPTR(8) + 18 * 8 GO TO 900 C ****************************************************************** C PROCESS MERCATOR GRIDS C C ------------------- BYTE 7-8 NR OF POINTS ALONG LATITUDE CIRCLE 4000 CONTINUE CALL GBYTE (MSGA,KGDS(2),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 9-10 NR OF POINTS ALONG LONG MERIDIAN CALL GBYTE (MSGA,KGDS(3),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 11-13 LATITUE OF ORIGIN CALL GBYTE (MSGA,KGDS(4),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(4),8388608).NE.0) THEN KGDS(4) = - IAND(KGDS(4),8388607) END IF C ------------------- BYTE 14-16 LONGITUDE OF ORIGIN CALL GBYTE (MSGA,KGDS(5),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(5),8388608).NE.0) THEN KGDS(5) = - IAND(KGDS(5),8388607) END IF C ------------------- BYTE 17 RESOLUTION FLAG CALL GBYTE (MSGA,KGDS(6),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 18-20 LATITUDE OF EXTREME POINT CALL GBYTE (MSGA,KGDS(7),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(7),8388608).NE.0) THEN KGDS(7) = - IAND(KGDS(7),8388607) END IF C ------------------- BYTE 21-23 LONGITUDE OF EXTREME POINT CALL GBYTE (MSGA,KGDS(8),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(8),8388608).NE.0) THEN KGDS(8) = - IAND(KGDS(8),8388607) END IF C ------------------- BYTE 24-26 LATITUDE OF PROJECTION INTERSECTION CALL GBYTE (MSGA,KGDS(9),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(9),8388608).NE.0) THEN KGDS(9) = - IAND(KGDS(9),8388607) END IF C ------------------- BYTE 27 RESERVED CALL GBYTE (MSGA,KGDS(10),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 28 SCANNING MODE CALL GBYTE (MSGA,KGDS(11),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 29-31 LONGITUDINAL DIR INCREMENT CALL GBYTE (MSGA,KGDS(12),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(12),8388608).NE.0) THEN KGDS(12) = - IAND(KGDS(12),8388607) END IF C ------------------- BYTE 32-34 LATITUDINAL DIR INCREMENT CALL GBYTE (MSGA,KGDS(13),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(13),8388608).NE.0) THEN KGDS(13) = - IAND(KGDS(13),8388607) END IF C ------------------- BYTE 35-42 RESERVED C SKIP TO START OF BYTE 43 KPTR(8) = KPTR(8) + 8 * 8 C ------------------- GO TO 900 C ****************************************************************** C PROCESS LAMBERT CONFORMAL C C ------------------- BYTE 7-8 NR OF POINTS ALONG X-AXIS 5000 CONTINUE CALL GBYTE (MSGA,KGDS(2),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 9-10 NR OF POINTS ALONG Y-AXIS CALL GBYTE (MSGA,KGDS(3),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- BYTE 11-13 LATITUDE OF ORIGIN CALL GBYTE (MSGA,KGDS(4),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(4),8388608).NE.0) THEN KGDS(4) = - IAND(KGDS(4),8388607) END IF C ------------------- BYTE 14-16 LONGITUDE OF ORIGIN (LOWER LEFT) CALL GBYTE (MSGA,KGDS(5),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(5),8388608).NE.0) THEN KGDS(5) = - IAND(KGDS(5),8388607) END IF C ------------------- BYTE 17 RESOLUTION CALL GBYTE (MSGA,KGDS(6),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 18-20 LOV -ORIENTATION OF GRID CALL GBYTE (MSGA,KGDS(7),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(7),8388608).NE.0) THEN KGDS(7) = - IAND(KGDS(7),8388607) END IF C ------------------- BYTE 21-23 DX - X-DIR INCREMENT CALL GBYTE (MSGA,KGDS(8),KPTR(8),24) KPTR(8) = KPTR(8) + 24 C ------------------- BYTE 24-26 DY - Y-DIR INCREMENT CALL GBYTE (MSGA,KGDS(9),KPTR(8),24) KPTR(8) = KPTR(8) + 24 C ------------------- BYTE 27 PROJECTION CENTER FLAG CALL GBYTE (MSGA,KGDS(10),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 28 SCANNING MODE CALL GBYTE (MSGA,KGDS(11),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ------------------- BYTE 29-31 LATIN1 - 1ST LAT FROM POLE CALL GBYTE (MSGA,KGDS(12),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(12),8388608).NE.0) THEN KGDS(12) = - IAND(KGDS(12),8388607) END IF C ------------------- BYTE 32-34 LATIN2 - 2ND LAT FROM POLE CALL GBYTE (MSGA,KGDS(13),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(13),8388608).NE.0) THEN KGDS(13) = - IAND(KGDS(13),8388607) END IF C ------------------- BYTE 35-37 LATITUDE OF SOUTHERN POLE CALL GBYTE (MSGA,KGDS(14),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(14),8388608).NE.0) THEN KGDS(14) = - IAND(KGDS(14),8388607) END IF C ------------------- BYTE 38-40 LONGITUDE OF SOUTHERN POLE CALL GBYTE (MSGA,KGDS(15),KPTR(8),24) KPTR(8) = KPTR(8) + 24 IF (IAND(KGDS(15),8388608).NE.0) THEN KGDS(15) = - IAND(KGDS(15),8388607) END IF C ------------------- BYTE 41-42 RESERVED CALL GBYTE (MSGA,KGDS(16),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ------------------- 900 CONTINUE C C MORE CODE FOR GRIDS WITH PL C IF (KGDS(19).EQ.0.OR.KGDS(19).EQ.255) THEN IF (KGDS(20).NE.255) THEN ISUM = 0 KPTR(8) = NSAVE + (KGDS(20) - 1) * 8 CALL GBYTES (MSGA,KGDS(22),KPTR(8),16,0,KGDS(3)) DO 910 J = 1, KGDS(3) ISUM = ISUM + KGDS(21+J) 910 CONTINUE KGDS(21) = ISUM END IF END IF RETURN END SUBROUTINE FI634(MSGA,KPTR,KPDS,KGDS,KBMS,KRET) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: FI634 EXTRACT OR GENERATE BIT MAP FOR OUTPUT C PRGMMR: BILL CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 C C ABSTRACT: IF BIT MAP SEC IS AVAILABLE IN GRIB MESSAGE, EXTRACT C FOR PROGRAM USE, OTHERWISE GENERATE AN APPROPRIATE BIT MAP. C C PROGRAM HISTORY LOG: C 91-09-13 CAVANAUGH C 91-11-12 CAVANAUGH MODIFIED SIZE OF ECMWF GRIDS 5 - 8. C 95-10-31 IREDELL REMOVED SAVES AND PRINTS C C USAGE: CALL FI634(MSGA,KPTR,KPDS,KGDS,KBMS,KRET) C INPUT ARGUMENT LIST: C MSGA - BUFR MESSAGE C KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS C (1) - TOTAL LENGTH OF GRIB MESSAGE C (2) - LENGTH OF INDICATOR (SECTION 0) C (3) - LENGTH OF PDS (SECTION 1) C (4) - LENGTH OF GDS (SECTION 2) C (5) - LENGTH OF BMS (SECTION 3) C (6) - LENGTH OF BDS (SECTION 4) C (7) - VALUE OF CURRENT BYTE C (8) - BIT POINTER C (9) - GRIB START BIT NR C (10) - GRIB/GRID ELEMENT COUNT C (11) - NR UNUSED BITS AT END OF SECTION 3 C (12) - BIT MAP FLAG C (13) - NR UNUSED BITS AT END OF SECTION 2 C (14) - BDS FLAGS C (15) - NR UNUSED BITS AT END OF SECTION 4 C KPDS - ARRAY CONTAINING PDS ELEMENTS. C (1) - ID OF CENTER C (2) - MODEL IDENTIFICATION C (3) - GRID IDENTIFICATION C (4) - GDS/BMS FLAG C (5) - INDICATOR OF PARAMETER C (6) - TYPE OF LEVEL C (7) - HEIGHT/PRESSURE , ETC OF LEVEL C (8) - YEAR OF CENTURY C (9) - MONTH OF YEAR C (10) - DAY OF MONTH C (11) - HOUR OF DAY C (12) - MINUTE OF HOUR C (13) - INDICATOR OF FORECAST TIME UNIT C (14) - TIME RANGE 1 C (15) - TIME RANGE 2 C (16) - TIME RANGE FLAG C (17) - NUMBER INCLUDED IN AVERAGE C C OUTPUT ARGUMENT LIST: C KBMS - BITMAP DESCRIBING LOCATION OF OUTPUT ELEMENTS. C KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS C SEE INPUT LIST C KRET - ERROR RETURN C C REMARKS: C KRET = 0 - NO ERROR C = 5 - GRID NOT AVAIL FOR CENTER INDICATED C =10 - INCORRECT CENTER INDICATOR C =12 - BYTES 5-6 ARE NOT ZERO IN BMS, PREDEFINED BIT MAP C NOT PROVIDED BY THIS CENTER C C SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT. C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C MACHINE: HDS9000 C C$$$ C C INCOMING MESSAGE HOLDER CHARACTER*1 MSGA(*) C C BIT MAP LOGICAL*1 KBMS(*) C C ARRAY OF POINTERS AND COUNTERS INTEGER KPTR(*) C ARRAY OF POINTERS AND COUNTERS INTEGER KPDS(*) INTEGER KGDS(*) C INTEGER KRET INTEGER MASK(8) C ----------------------GRID 21 AND GRID 22 ARE THE SAME LOGICAL*1 GRD21( 1369) C ----------------------GRID 23 AND GRID 24 ARE THE SAME LOGICAL*1 GRD23( 1369) LOGICAL*1 GRD25( 1368) LOGICAL*1 GRD26( 1368) C ----------------------GRID 27 AND GRID 28 ARE THE SAME C ----------------------GRID 29 AND GRID 30 ARE THE SAME C ----------------------GRID 33 AND GRID 34 ARE THE SAME LOGICAL*1 GRD50( 1188) C -----------------------GRID 61 AND GRID 62 ARE THE SAME LOGICAL*1 GRD61( 4186) C -----------------------GRID 63 AND GRID 64 ARE THE SAME LOGICAL*1 GRD63( 4186) C LOGICAL*1 GRD70(16380)/16380*.TRUE./ C ------------------------------------------------------------- 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/ C LINE 1-4 & 7*.FALSE.,22*.TRUE.,14*.FALSE.,22*.TRUE., & 14*.FALSE.,22*.TRUE.,14*.FALSE.,22*.TRUE.,7*.FALSE., C LINE 5-8 & 6*.FALSE.,24*.TRUE.,12*.FALSE.,24*.TRUE., & 12*.FALSE.,24*.TRUE.,12*.FALSE.,24*.TRUE.,6*.FALSE., C LINE 9-12 & 5*.FALSE.,26*.TRUE.,10*.FALSE.,26*.TRUE., & 10*.FALSE.,26*.TRUE.,10*.FALSE.,26*.TRUE.,5*.FALSE., C LINE 13-16 & 4*.FALSE.,28*.TRUE., 8*.FALSE.,28*.TRUE., & 8*.FALSE.,28*.TRUE., 8*.FALSE.,28*.TRUE.,4*.FALSE., C LINE 17-20 & 3*.FALSE.,30*.TRUE., 6*.FALSE.,30*.TRUE., & 6*.FALSE.,30*.TRUE., 6*.FALSE.,30*.TRUE.,3*.FALSE., C LINE 21-24 & 2*.FALSE.,32*.TRUE., 4*.FALSE.,32*.TRUE., & 4*.FALSE.,32*.TRUE., 4*.FALSE.,32*.TRUE.,2*.FALSE., C LINE 25-28 & .FALSE.,34*.TRUE., 2*.FALSE.,34*.TRUE., & 2*.FALSE.,34*.TRUE., 2*.FALSE.,34*.TRUE., .FALSE., C 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/ C C PRINT *,'FI634' IF (IAND(KPDS(4),64).EQ.64) THEN C C SET UP BIT POINTER C SECTION 0 SECTION 1 SECTION 2 KPTR(8) = KPTR(9) + (KPTR(2)*8) + (KPTR(3)*8) + (KPTR(4)*8) + 24 C C BYTE 4 NUMBER OF UNUSED BITS AT END OF SECTION 3 C CALL GBYTE (MSGA,KPTR(11),KPTR(8),8) KPTR(8) = KPTR(8) + 8 C C BYTE 5,6 TABLE REFERENCE IF 0, BIT MAP FOLLOWS C CALL GBYTE (MSGA,KPTR(12),KPTR(8),16) KPTR(8) = KPTR(8) + 16 C IF TABLE REFERENCE = 0, EXTRACT BIT MAP IF (KPTR(12).EQ.0) THEN C CALCULATE NR OF BITS IN BIT MAP IBITS = (KPTR(5) - 6) * 8 - KPTR(11) KPTR(10) = IBITS IF (KPDS(3).EQ.21.OR.KPDS(3).EQ.22.OR.KPDS(3).EQ.25. * OR.KPDS(3).EQ.61.OR.KPDS(3).EQ.62) THEN C NORTHERN HEMISPHERE 21, 22, 25, 61, 62 DO 2122 I = 1, IBITS CALL GBYTE (MSGA,ICHK,KPTR(8),1) KPTR(8) = KPTR(8) + 1 IF (ICHK.NE.0) THEN KBMS(I) = .TRUE. ELSE KBMS(I) = .FALSE. END IF 2122 CONTINUE IF (KPDS(3).EQ.25) THEN KADD = 71 ELSE IF (KPDS(3).EQ.61.OR.KPDS(3).EQ.62) THEN KADD = 90 ELSE KADD = 36 END IF DO 25 I = 1, KADD KBMS(I+IBITS) = .FALSE. 25 CONTINUE KPTR(10) = KPTR(10) + KADD RETURN ELSE IF (KPDS(3).EQ.23.OR.KPDS(3).EQ.24.OR.KPDS(3).EQ.26. * OR.KPDS(3).EQ.63.OR.KPDS(3).EQ.64) THEN C SOUTHERN HEMISPHERE 23, 24, 26, 63, 64 IF (KPDS(3).EQ.26) THEN KADD = 72 ELSE IF (KPDS(3).EQ.63.OR.KPDS(3).EQ.64) THEN KADD = 91 ELSE KADD = 37 END IF DO 26 I = 1, KADD KBMS(I+IBITS) = .FALSE. 26 CONTINUE DO 2324 I = 1, IBITS CALL GBYTE (MSGA,ICHK,KPTR(8),1) KPTR(8) = KPTR(8) + 1 IF (ICHK.NE.0) THEN KBMS(I) = .TRUE. ELSE KBMS(I) = .FALSE. END IF 2324 CONTINUE KPTR(10) = KPTR(10) + KADD - 1 RETURN ELSE IF (KPDS(3).EQ.50) THEN KPAD = 7 KIN = 22 KBITS = 0 DO 55 I = 1, 7 DO 54 J = 1, 4 DO 51 K = 1, KPAD KBITS = KBITS + 1 KBMS(KBITS) = .FALSE. 51 CONTINUE DO 52 K = 1, KIN CALL GBYTE (MSGA,ICHK,KPTR(8),1) KPTR(8) = KPTR(8) + 1 KBITS = KBITS + 1 IF (ICHK.NE.0) THEN KBMS(KBITS) = .TRUE. ELSE KBMS(KBITS) = .FALSE. END IF 52 CONTINUE DO 53 K = 1, KPAD KBITS = KBITS + 1 KBMS(KBITS) = .FALSE. 53 CONTINUE 54 CONTINUE KIN = KIN + 2 KPAD = KPAD - 1 55 CONTINUE DO 57 II = 1, 5 DO 56 J = 1, KIN CALL GBYTE (MSGA,ICHK,KPTR(8),1) KPTR(8) = KPTR(8) + 1 KBITS = KBITS + 1 IF (ICHK.NE.0) THEN KBMS(KBITS) = .TRUE. ELSE KBMS(KBITS) = .FALSE. END IF 56 CONTINUE 57 CONTINUE ELSE C EXTRACT BIT MAP FROM BMS FOR OTHER GRIDS DO 100 I = 1, IBITS CALL GBYTE (MSGA,ICHK,KPTR(8),1) KPTR(8) = KPTR(8) + 1 IF (ICHK.NE.0) THEN KBMS(I) = .TRUE. ELSE KBMS(I) = .FALSE. END IF 100 CONTINUE END IF RETURN ELSE C PRINT *,'FI634-NO PREDEFINED BIT MAP PROVIDED BY THIS CENTER' KRET = 12 RETURN END IF C END IF KRET = 0 C ------------------------------------------------------- C PROCESS NON-STANDARD GRID C ------------------------------------------------------- IF (KPDS(3).EQ.255) THEN C PRINT *,'NON STANDARD GRID, CENTER = ',KPDS(1) J = KGDS(2) * KGDS(3) KPTR(10) = J DO 600 I = 1, J KBMS(I) = .TRUE. 600 CONTINUE RETURN END IF C ------------------------------------------------------- C CHECK INTERNATIONAL SET C ------------------------------------------------------- IF (KPDS(3).EQ.21.OR.KPDS(3).EQ.22) THEN C ----- INT'L GRIDS 21, 22 - MAP SIZE 1369 J = 1369 KPTR(10) = J CALL FI637(*820,J,KPDS,KGDS,KRET) DO 3021 I = 1, 1369 KBMS(I) = GRD21(I) 3021 CONTINUE RETURN ELSE IF (KPDS(3).EQ.23.OR.KPDS(3).EQ.24) THEN C ----- INT'L GRIDS 23, 24 - MAP SIZE 1369 J = 1369 KPTR(10) = J CALL FI637(*820,J,KPDS,KGDS,KRET) DO 3023 I = 1, 1369 KBMS(I) = GRD23(I) 3023 CONTINUE RETURN ELSE IF (KPDS(3).EQ.25) THEN C ----- INT'L GRID 25 - MAP SIZE 1368 J = 1368 KPTR(10) = J CALL FI637(*820,J,KPDS,KGDS,KRET) DO 3025 I = 1, 1368 KBMS(I) = GRD25(I) 3025 CONTINUE RETURN ELSE IF (KPDS(3).EQ.26) THEN C ----- INT'L GRID 26 - MAP SIZE 1368 J = 1368 KPTR(10) = J CALL FI637(*820,J,KPDS,KGDS,KRET) DO 3026 I = 1, 1368 KBMS(I) = GRD26(I) 3026 CONTINUE RETURN ELSE IF (KPDS(3).GE.37.AND.KPDS(3).LE.44) THEN C ----- INT'L GRID 37-44 - MAP SIZE 3447 J = 3447 GO TO 800 ELSE IF (KPDS(1).EQ.7.AND.KPDS(3).EQ.50) THEN C ----- INT'L GRIDS 50 - MAP SIZE 964 J = 1188 KPTR(10) = J CALL FI637(*890,J,KPDS,KGDS,KRET) DO 3050 I = 1, J KBMS(I) = GRD50(I) 3050 CONTINUE RETURN ELSE IF (KPDS(3).EQ.61.OR.KPDS(3).EQ.62) THEN C ----- INT'L GRIDS 61, 62 - MAP SIZE 4186 J = 4186 KPTR(10) = J CALL FI637(*820,J,KPDS,KGDS,KRET) DO 3061 I = 1, 4186 KBMS(I) = GRD61(I) 3061 CONTINUE RETURN ELSE IF (KPDS(3).EQ.63.OR.KPDS(3).EQ.64) THEN C ----- INT'L GRIDS 63, 64 - MAP SIZE 4186 J = 4186 KPTR(10) = J CALL FI637(*820,J,KPDS,KGDS,KRET) DO 3063 I = 1, 4186 KBMS(I) = GRD63(I) 3063 CONTINUE RETURN END IF C ------------------------------------------------------- C CHECK UNITED STATES SET C ------------------------------------------------------- IF (KPDS(1).EQ.7) THEN IF (KPDS(3).LT.100) THEN IF (KPDS(3).EQ.1) THEN C ----- U.S. GRID 1 - MAP SIZE 1679 J = 1679 GO TO 800 END IF IF (KPDS(3).EQ.2) THEN C ----- U.S. GRID 2 - MAP SIZE 10512 J = 10512 GO TO 800 ELSE IF (KPDS(3).EQ.3) THEN C ----- U.S. GRID 3 - MAP SIZE 65160 J = 65160 GO TO 800 ELSE IF (KPDS(3).EQ.4) THEN C ----- U.S. GRID 4 - MAP SIZE 259920 J = 259920 GO TO 800 ELSE IF (KPDS(3).EQ.5) THEN C ----- U.S. GRID 5 - MAP SIZE 3021 J = 3021 GO TO 800 ELSE IF (KPDS(3).EQ.6) THEN C ----- U.S. GRID 6 - MAP SIZE 2385 J = 2385 GO TO 800 ELSE IF (KPDS(3).EQ.8) THEN C ----- U.S. GRID 8 - MAP SIZE 5104 J = 5104 GO TO 800 ELSE IF (KPDS(3).EQ.27.OR.KPDS(3).EQ.28) THEN C ----- U.S. GRIDS 27, 28 - MAP SIZE 4225 J = 4225 GO TO 800 ELSE IF (KPDS(3).EQ.29.OR.KPDS(3).EQ.30) THEN C ----- U.S. GRIDS 29,30 - MAP SIZE 5365 J = 5365 GO TO 800 ELSE IF (KPDS(3).EQ.33.OR.KPDS(3).EQ.34) THEN C ----- U.S GRID 33, 34 - MAP SIZE 8326 J = 8326 GO TO 800 ELSE IF (KPDS(3).GE.37.AND.KPDS(3).LE.44) THEN C ----- U.S. GRID 37-44 - MAP SIZE 3447 J = 3447 GO TO 800 ELSE IF (KPDS(3).EQ.45) THEN C ----- U.S. GRID 45 - MAP SIZE 41760 J = 41760 GO TO 800 ELSE IF (KPDS(3).EQ.53) THEN C ----- U.S. GRID 53 - MAP SIZE 5967 J = 5967 GO TO 800 ELSE IF (KPDS(3).EQ.55.OR.KPDS(3).EQ.56) THEN C ----- U.S GRID 55, 56 - MAP SIZE 6177 J = 6177 GO TO 800 ELSE IF (KPDS(3).GE.67.AND.KPDS(3).LE.71) THEN C ----- U.S GRID 67-71 - MAP SIZE 13689 J = 13689 GO TO 800 ELSE IF (KPDS(3).EQ.72) THEN C ----- U.S GRID 72 - MAP SIZE 406 J = 406 GO TO 800 ELSE IF (KPDS(3).EQ.73) THEN C ----- U.S GRID 73 - MAP SIZE 13056 J = 13056 GO TO 800 ELSE IF (KPDS(3).EQ.74) THEN C ----- U.S GRID 74 - MAP SIZE 10800 J = 10800 GO TO 800 ELSE IF (KPDS(3).GE.75.AND.KPDS(3).LE.77) THEN C ----- U.S GRID 75-77 - MAP SIZE 12321 J = 12321 GO TO 800 ELSE IF (KPDS(3).EQ.85.OR.KPDS(3).EQ.86) THEN C ----- U.S GRID 85,86 - MAP SIZE 32400 J = 32400 GO TO 800 ELSE IF (KPDS(3).EQ.87) THEN C ----- U.S GRID 87 - MAP SIZE 5022 J = 5022 GO TO 800 ELSE IF (KPDS(3).EQ.90) THEN C ----- U.S GRID 90 - MAP SIZE 12902 J = 12902 GO TO 800 ELSE IF (KPDS(3).EQ.91) THEN C ----- U.S GRID 91 - MAP SIZE 25803 J = 25803 GO TO 800 ELSE IF (KPDS(3).EQ.92) THEN C ----- U.S GRID 92 - MAP SIZE 24162 J = 24162 GO TO 800 ELSE IF (KPDS(3).EQ.93) THEN C ----- U.S GRID 93 - MAP SIZE 48323 J = 48323 GO TO 800 ELSE IF (KPDS(3).EQ.94) THEN C ----- U.S GRID 94 - MAP SIZE 48916 J = 48916 GO TO 800 ELSE IF (KPDS(3).EQ.95) THEN C ----- U.S GRID 95 - MAP SIZE 97831 J = 97831 GO TO 800 ELSE IF (KPDS(3).EQ.96) THEN C ----- U.S GRID 96 - MAP SIZE 41630 J = 41630 GO TO 800 ELSE IF (KPDS(3).EQ.97) THEN C ----- U.S GRID 97 - MAP SIZE 83259 J = 83259 GO TO 800 ELSE IF (KPDS(3).EQ.98) THEN C ----- U.S GRID 98 - MAP SIZE 18048 J = 18048 GO TO 800 END IF ELSE IF (KPDS(3).GE.100.AND.KPDS(3).LT.200) THEN IF (KPDS(3).EQ.100) THEN C ----- U.S. GRID 100 - MAP SIZE 6889 J = 6889 GO TO 800 ELSE IF (KPDS(3).EQ.101) THEN C ----- U.S. GRID 101 - MAP SIZE 10283 J = 10283 GO TO 800 ELSE IF (KPDS(3).EQ.103) THEN C ----- U.S. GRID 103 - MAP SIZE 3640 J = 3640 GO TO 800 ELSE IF (KPDS(3).EQ.104) THEN C ----- U.S. GRID 104 - MAP SIZE 16170 J = 16170 GO TO 800 ELSE IF (KPDS(3).EQ.105) THEN C ----- U.S. GRID 105 - MAP SIZE 6889 J = 6889 GO TO 800 ELSE IF (KPDS(3).EQ.106) THEN C ----- U.S. GRID 106 - MAP SIZE 19305 J = 19305 GO TO 800 ELSE IF (KPDS(3).EQ.107) THEN C ----- U.S. GRID 107 - MAP SIZE 11040 J = 11040 GO TO 800 ELSE IF (KPDS(3).EQ.126) THEN C ----- U.S. GRID 126 - MAP SIZE 72960 J = 72960 GO TO 800 ELSE IF (KPDS(3).EQ.196) THEN C ----- U.S. GRID 196 - MAP SIZE 45903 J = 45903 GO TO 800 ELSE IF (IAND(KPDS(4),128).EQ.128) THEN C ----- U.S. NON-STANDARD GRID GO TO 895 END IF ELSE IF (KPDS(3).GE.200) THEN IF (KPDS(3).EQ.201) THEN J = 4225 GO TO 800 ELSE IF (KPDS(3).EQ.202) THEN J = 2795 GO TO 800 ELSE IF (KPDS(3).EQ.203.OR.KPDS(3).EQ.205) THEN J = 1755 GO TO 800 ELSE IF (KPDS(3).EQ.204) THEN J = 6324 GO TO 800 ELSE IF (KPDS(3).EQ.206) THEN J = 2091 GO TO 800 ELSE IF (KPDS(3).EQ.207) THEN J = 1715 GO TO 800 ELSE IF (KPDS(3).EQ.208) THEN J = 783 GO TO 800 ELSE IF (KPDS(3).EQ.209) THEN J = 8181 GO TO 800 ELSE IF (KPDS(3).EQ.210) THEN J = 625 GO TO 800 ELSE IF (KPDS(3).EQ.211) THEN J = 6045 GO TO 800 ELSE IF (KPDS(3).EQ.212) THEN J = 23865 GO TO 800 ELSE IF (KPDS(3).EQ.213) THEN J = 10965 GO TO 800 ELSE IF (KPDS(3).EQ.214) THEN J = 6693 GO TO 800 ELSE IF (KPDS(3).EQ.215) THEN J = 94833 GO TO 800 ELSE IF (KPDS(3).EQ.216) THEN J = 14873 GO TO 800 ELSE IF (IAND(KPDS(4),128).EQ.128) THEN GO TO 895 END IF KRET = 5 RETURN END IF END IF C ------------------------------------------------------- C CHECK JAPAN METEOROLOGICAL AGENCY SET C ------------------------------------------------------- IF (KPDS(1).EQ.34) THEN IF (IAND(KPDS(4),128).EQ.128) THEN C PRINT *,'JMA MAP IS NOT PREDEFINED, THE GDS WILL' C PRINT *,'BE USED TO UNPACK THE DATA, MAP = ',KPDS(3) GO TO 900 END IF END IF C ------------------------------------------------------- C CHECK CANADIAN SET C ------------------------------------------------------- IF (KPDS(1).EQ.54) THEN IF (IAND(KPDS(4),128).EQ.128) THEN C PRINT *,'CANADIAN MAP IS NOT PREDEFINED, THE GDS WILL' C PRINT *,'BE USED TO UNPACK THE DATA, MAP = ',KPDS(3) GO TO 900 END IF END IF C ------------------------------------------------------- C CHECK FNOC SET C ------------------------------------------------------- IF (KPDS(1).EQ.58) THEN IF (KPDS(3).EQ.220.OR.KPDS(3).EQ.221) THEN C 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).EQ.223) THEN C 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).EQ.128) THEN C PRINT *,'FNOC MAP IS NOT PREDEFINED, THE GDS WILL' C PRINT *,'BE USED TO UNPACK THE DATA, MAP = ',KPDS(3) GO TO 900 END IF END IF C ------------------------------------------------------- C CHECK UKMET SET C ------------------------------------------------------- IF (KPDS(1).EQ.74) THEN IF (IAND(KPDS(4),128).EQ.128) THEN GO TO 820 END IF END IF C ------------------------------------------------------- C CHECK ECMWF SET C ------------------------------------------------------- IF (KPDS(1).EQ.98) THEN IF (KPDS(3).GE.1.AND.KPDS(3).LE.12) THEN IF (KPDS(3).GE.5.AND.KPDS(3).LE.8) THEN J = 1073 ELSE J = 1369 END IF KPTR(10) = J CALL FI637(*810,J,KPDS,KGDS,KRET) DO 1000 I = 1, J KBMS(I) = .TRUE. 1000 CONTINUE RETURN ELSE IF (KPDS(3).GE.13.AND.KPDS(3).LE.16) THEN J = 361 KPTR(10) = J CALL FI637(*810,J,KPDS,KGDS,KRET) DO 1013 I = 1, J KBMS(I) = .TRUE. 1013 CONTINUE RETURN ELSE IF (IAND(KPDS(4),128).EQ.128) THEN GO TO 810 ELSE KRET = 5 RETURN END IF ELSE C PRINT *,'CENTER ',KPDS(1),' IS NOT DEFINED' IF (IAND(KPDS(4),128).EQ.128) THEN C PRINT *,'GDS WILL BE USED TO UNPACK THE DATA', C * ' MAP = ',KPDS(3) GO TO 900 ELSE KRET = 10 RETURN END IF END IF C ======================================= C 800 CONTINUE KPTR(10) = J CALL FI637 (*801,J,KPDS,KGDS,KRET) DO 2201 I = 1, J KBMS(I) = .TRUE. 2201 CONTINUE RETURN 801 CONTINUE C C ----- THE MAP HAS A GDS, BYTE 7 OF THE (PDS) THE GRID IDENTIFICATION C ----- IS NOT 255, THE SIZE OF THE GRID IS NOT THE SAME AS THE C ----- PREDEFINED SIZES OF THE U.S. GRIDS, OR KNOWN GRIDS OF THE C ----- OF THE OTHER CENTERS. THE GRID CAN BE UNKNOWN, OR FROM AN C ----- UNKNOWN CENTER, WE WILL USE THE INFORMATION IN THE GDS TO MAKE C ----- A BIT MAP. C 810 CONTINUE C PRINT *,'ECMWF PREDEFINED MAP SIZE DOES NOT MATCH, I WILL USE' GO TO 895 C 820 CONTINUE C PRINT *,'U.K. MET PREDEFINED MAP SIZE DOES NOT MATCH, I WILL USE' GO TO 895 C 890 CONTINUE C PRINT *,'PREDEFINED MAP SIZE DOES NOT MATCH, I WILL USE' 895 CONTINUE C PRINT *,'THE GDS TO UNPACK THE DATA, MAP TYPE = ',KPDS(3) C 900 CONTINUE J = KGDS(2) * KGDS(3) C AFOS AFOS AFOS SPECIAL CASE C INVOLVES NEXT SINGLE STATEMENT ONLY IF (KPDS(3).EQ.211) KRET = 0 KPTR(10) = J DO 2203 I = 1, J KBMS(I) = .TRUE. 2203 CONTINUE C PRINT *,'EXIT FI634' RETURN END SUBROUTINE FI635(MSGA,KPTR,KPDS,KGDS,KBMS,DATA,KRET) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: FI635 EXTRACT GRIB DATA ELEMENTS FROM BDS C PRGMMR: BILL CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 C C ABSTRACT: EXTRACT GRIB DATA FROM BINARY DATA SECTION AND PLACE C INTO OUTPUT ARRAY IN PROPER POSITION. C C PROGRAM HISTORY LOG: C 91-09-13 CAVANAUGH C 94-04-01 CAVANAUGH MODIFIED CODE TO INCLUDE DECIMAL SCALING WHEN C CALCULATING THE VALUE OF DATA POINTS SPECIFIED C AS BEING EQUAL TO THE REFERENCE VALUE C 94-11-10 FARLEY INCREASED MXSIZE FROM 72960 TO 260000 C FOR .5 DEGREE SST ANALYSIS FIELDS C 95-10-31 IREDELL REMOVED SAVES AND PRINTS C C USAGE: CALL FI635(MSGA,KPTR,KPDS,KGDS,KBMS,DATA,KRET) C INPUT ARGUMENT LIST: C MSGA - ARRAY CONTAINING GRIB MESSAGE C KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS C (1) - TOTAL LENGTH OF GRIB MESSAGE C (2) - LENGTH OF INDICATOR (SECTION 0) C (3) - LENGTH OF PDS (SECTION 1) C (4) - LENGTH OF GDS (SECTION 2) C (5) - LENGTH OF BMS (SECTION 3) C (6) - LENGTH OF BDS (SECTION 4) C (7) - VALUE OF CURRENT BYTE C (8) - BIT POINTER C (9) - GRIB START BIT NR C (10) - GRIB/GRID ELEMENT COUNT C (11) - NR UNUSED BITS AT END OF SECTION 3 C (12) - BIT MAP FLAG C (13) - NR UNUSED BITS AT END OF SECTION 2 C (14) - BDS FLAGS C (15) - NR UNUSED BITS AT END OF SECTION 4 C KPDS - ARRAY CONTAINING PDS ELEMENTS. C SEE INITIAL ROUTINE C KBMS - BITMAP DESCRIBING LOCATION OF OUTPUT ELEMENTS. C C OUTPUT ARGUMENT LIST: C KBDS - INFORMATION EXTRACTED FROM BINARY DATA SECTION C KBDS(1) - N1 C KBDS(2) - N2 C KBDS(3) - P1 C KBDS(4) - P2 C KBDS(5) - BIT POINTER TO 2ND ORDER WIDTHS C KBDS(6) - " " " " " BIT MAPS C KBDS(7) - " " " FIRST ORDER VALUES C KBDS(8) - " " " SECOND ORDER VALUES C KBDS(9) - " " START OF BDS C KBDS(10) - " " MAIN BIT MAP C KBDS(11) - BINARY SCALING C KBDS(12) - DECIMAL SCALING C KBDS(13) - BIT WIDTH OF FIRST ORDER VALUES C KBDS(14) - BIT MAP FLAG C 0 = NO SECOND ORDER BIT MAP C 1 = SECOND ORDER BIT MAP PRESENT C KBDS(15) - SECOND ORDER BIT WIDTH C KBDS(16) - CONSTANT / DIFFERENT WIDTHS C 0 = CONSTANT WIDTHS C 1 = DIFFERENT WIDTHS C KBDS(17) - SINGLE DATUM / MATRIX C 0 = SINGLE DATUM AT EACH GRID POINT C 1 = MATRIX OF VALUES AT EACH GRID POINT C (18-20)- UNUSED C C DATA - REAL*4 ARRAY OF GRIDDED ELEMENTS IN GRIB MESSAGE. C KPTR - ARRAY CONTAINING STORAGE FOR FOLLOWING PARAMETERS C SEE INPUT LIST C KRET - ERROR RETURN C C REMARKS: C ERROR RETURN C 3 = UNPACKED FIELD IS LARGER THAN 65160 C 6 = DOES NOT MATCH NR OF ENTRIES FOR THIS GRIB/GRID C 7 = NUMBER OF BITS IN FILL TOO LARGE C C SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT. C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C MACHINE: HDS9000 C C$$$ C ************************************************************* C ON A PC THIS CAN BE CHANGED TO A SMALLER SIZE TO BETTER FIT C THE DOS MEMORY LIMIT OF 640K BYTES. YOU COULD DO THIS C FOR MICROSOFT 5.0. A PC 32 BIT FORTRAN COMPILER C WOULD NOT NEED THIS CHANGE. IF NONE OF YOUR GRIB RECORDS C IS LARGER THAN 20000, SET MXSIZE TO 20000. C ************************************************************* C PARAMETER (MXSIZE=1200000) C CHARACTER*1 MSGA(*) CHARACTER*1 KK(8) CHARACTER*1 CKREF(8) C LOGICAL*1 KBMS(*) C INTEGER KPDS(*) INTEGER KGDS(*) INTEGER KBDS(20) INTEGER KPTR(*) INTEGER NRBITS INTEGER KREF INTEGER KKK INTEGER KSAVE(MXSIZE) INTEGER KSCALE C REAL DATA(*) REAL REFNCE REAL SCALE REAL REALKK C EQUIVALENCE (CKREF(1),KREF,REFNCE) EQUIVALENCE (KK(1),KKK,REALKK) C C C CHANGED HEX VALUES TO DECIMAL TO MAKE CODE MORE PORTABLE C C ************************************************************* C PRINT *,'ENTER FI635' C SET UP BIT POINTER KPTR(8) = KPTR(9) + (KPTR(2)*8) + (KPTR(3)*8) + (KPTR(4)*8) * + (KPTR(5)*8) + 24 C ------------- EXTRACT FLAGS C BYTE 4 CALL GBYTE(MSGA,KPTR(14),KPTR(8),4) KPTR(8) = KPTR(8) + 4 C --------- NR OF UNUSED BITS IN SECTION 4 CALL GBYTE(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) C ------------- GET SCALE FACTOR C BYTES 5,6 C CHECK SIGN CALL GBYTE (MSGA,KSIGN,KPTR(8),1) KPTR(8) = KPTR(8) + 1 C GET ABSOLUTE SCALE VALUE CALL GBYTE (MSGA,KSCALE,KPTR(8),15) KPTR(8) = KPTR(8) + 15 IF (KSIGN.GT.0) THEN KSCALE = - KSCALE END IF SCALE = 2.0**KSCALE C ------------ GET REFERENCE VALUE C BYTES 7,10 CALL GBYTE (MSGA,KREF,KPTR(8),32) KPTR(8) = KPTR(8) + 32 C C THE NEXT CODE WILL CONVERT THE IBM370 FLOATING POINT C TO THE FLOATING POINT USED ON YOUR COMPUTER. C C 1ST TEST TO SEE IN ON 32 OR 64 BIT WORD MACHINE C LW = 4 OR 8; IF 8 MAY BE A CRAY C CALL W3FI01(LW) IF (LW.EQ.4) THEN CALL GBYTE (CKREF,JSGN,0,1) CALL GBYTE (CKREF,JEXP,1,7) CALL GBYTE (CKREF,IFR,8,24) ELSE CALL GBYTE (CKREF,JSGN,32,1) CALL GBYTE (CKREF,JEXP,33,7) CALL GBYTE (CKREF,IFR,40,24) ENDIF C PRINT *,109,JSGN,JEXP,IFR C 109 FORMAT (' JSGN,JEXP,IFR = ',3(1X,Z8)) IF (IFR.EQ.0) THEN REFNCE = 0.0 ELSE IF (JEXP.EQ.0.AND.IFR.EQ.0) THEN REFNCE = 0.0 ELSE REFNCE = FLOAT(IFR) * 16.0 ** (JEXP - 64 - 6) IF (JSGN.NE.0) REFNCE = - REFNCE END IF C PRINT *,'SCALE ',SCALE,' REF VAL ',KREF,REFNCE C ------------- NUMBER OF BITS SPECIFIED FOR EACH ENTRY C BYTE 11 CALL GBYTE (MSGA,KBITS,KPTR(8),8) KPTR(8) = KPTR(8) + 8 KBDS(4) = KBITS C KBDS(13) = KBITS IBYT12 = KPTR(8) C ------------------ IF THERE ARE NO EXTENDED FLAGS PRESENT C THIS IS WHERE DATA BEGINS AND AND THE PROCESSING C INCLUDED IN THE FOLLOWING IF...END IF C WILL BE SKIPPED C PRINT *,'BASIC FLAGS =',KPTR(14) ,IAND(KPTR(14),1) IF (IAND(KPTR(14),1).EQ.0) THEN C PRINT *,'NO EXTENDED FLAGS' ELSE C BYTES 12,13 CALL GBYTE (MSGA,KOCTET,KPTR(8),16) KPTR(8) = KPTR(8) + 16 C --------------------------- EXTENDED FLAGS C BYTE 14 CALL GBYTE (MSGA,KXFLAG,KPTR(8),8) C PRINT *,'HAVE EXTENDED FLAGS',KXFLAG KPTR(8) = KPTR(8) + 8 IF (IAND(KXFLAG,16).EQ.0) THEN C SECOND ORDER VALUES CONSTANT WIDTHS KBDS(16) = 0 ELSE C SECOND ORDER VALUES DIFFERENT WIDTHS KBDS(16) = 1 END IF IF (IAND (KXFLAG,32).EQ.0) THEN C NO SECONDARY BIT MAP KBDS(14) = 0 ELSE C HAVE SECONDARY BIT MAP KBDS(14) = 1 END IF IF (IAND (KXFLAG,64).EQ.0) THEN C SINGLE DATUM AT GRID POINT KBDS(17) = 0 ELSE C MATRIX OF VALUES AT GRID POINT KBDS(17) = 1 END IF C ---------------------- NR - FIRST DIMENSION (ROWS) OF EACH MATRIX C BYTES 15,16 CALL GBYTE (MSGA,NR,KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ---------------------- NC - SECOND DIMENSION (COLS) OF EACH MATRIX C BYTES 17,18 CALL GBYTE (MSGA,NC,KPTR(8),16) KPTR(8) = KPTR(8) + 16 C ---------------------- NRV - FIRST DIM COORD VALS C BYTE 19 CALL GBYTE (MSGA,NRV,KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ---------------------- NC1 - NR COEFF'S OR VALUES C BYTE 20 CALL GBYTE (MSGA,NC1,KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ---------------------- NCV - SECOND DIM COORD OR VALUE C BYTE 21 CALL GBYTE (MSGA,NCV,KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ---------------------- NC2 - NR COEFF'S OR VALS C BYTE 22 CALL GBYTE (MSGA,NC2,KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ---------------------- KPHYS1 - FIRST DIM PHYSICAL SIGNIF C BYTE 23 CALL GBYTE (MSGA,KPHYS1,KPTR(8),8) KPTR(8) = KPTR(8) + 8 C ---------------------- KPHYS2 - SECOND DIM PHYSICAL SIGNIF C BYTE 24 CALL GBYTE (MSGA,KPHYS2,KPTR(8),8) KPTR(8) = KPTR(8) + 8 C BYTES 25-N END IF IF (KBITS.EQ.0) THEN C HAVE NO BDS ENTRIES, ALL ENTRIES = REFNCE SCAL10 = 10.0 ** KPDS(22) SCAL10 = 1.0 / SCAL10 REFN10 = REFNCE * SCAL10 KENTRY = KPTR(10) DO 210 I = 1, KENTRY DATA(I) = 0.0 IF (KBMS(I)) THEN DATA(I) = REFN10 END IF 210 CONTINUE GO TO 900 END IF C PRINT *,'KEND ',KEND,' KPTR(8) ',KPTR(8),'KBITS ',KBITS KNR = (KEND - KPTR(8)) / KBITS C PRINT *,'NUMBER OF ENTRIES IN DATA ARRAY',KNR C -------------------- C CYCLE THRU BDS UNTIL HAVE USED ALL (SPECIFIED NUMBER) C ENTRIES. C ------------- UNUSED BITS IN DATA AREA C NUMBER OF BYTES IN DATA AREA NRBYTE = KPTR(6) - 11 C ------------- TOTAL NR OF USABLE BITS NRBITS = NRBYTE * 8 - KPTR(15) C ------------- TOTAL NR OF ENTRIES KENTRY = NRBITS / KBITS C MAX SIZE CHECK IF (KENTRY.GT.MXSIZE) THEN KRET = 3 RETURN END IF C C IF (IAND(KPTR(14),2).EQ.0) THEN C PRINT *,'SOURCE VALUES IN FLOATING POINT' C ELSE C PRINT *,'SOURCE VALUES IN INTEGER' C END IF C IF (IAND(KPTR(14),8).EQ.0) THEN C PRINT *,'PROCESSING GRID POINT DATA' IF (IAND(KPTR(14),4).EQ.0) THEN C PRINT *,' WITH SIMPLE PACKING' IF (IAND(KPTR(14),1).EQ.0) THEN C PRINT *,' WITH NO ADDITIONAL FLAGS' GO TO 4000 ELSE IF (IAND(KPTR(14),1).NE.0) THEN C PRINT *,' WITH ADDITIONAL FLAGS',KXFLAG IF (KBDS(17).EQ.0) THEN C PRINT *,' SINGLE DATUM EACH GRID PT' IF (KBDS(14).EQ.0) THEN C PRINT *,' NO SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (KBDS(14).NE.0) THEN C PRINT *,' SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF END IF ELSE IF (KBDS(17).NE.0) THEN C PRINT *,' MATRIX OF VALS EACH PT' IF (KBDS(14).EQ.0) THEN C PRINT *,' NO SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (KBDS(14).NE.0) THEN C PRINT *,' SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF END IF END IF END IF ELSE IF (IAND(KPTR(14),4).NE.0) THEN C PRINT *,' WITH COMPLEX/SECOND ORDER PACKING' IF (IAND(KPTR(14),1).EQ.0) THEN C PRINT *,' WITH NO ADDITIONAL FLAGS' ELSE IF (IAND(KPTR(14),1).NE.0) THEN C PRINT *,' WITH ADDITIONAL FLAGS' IF (KBDS(17).EQ.0) THEN C PRINT *,' SINGLE DATUM AT EACH PT' IF (KBDS(14).EQ.0) THEN C PRINT *,' NO SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF C ROW BY ROW - COL BY COL CALL FI636 (DATA,MSGA,KBMS, * REFNCE,KPTR,KPDS,KGDS) GO TO 900 ELSE IF (KBDS(14).NE.0) THEN C PRINT *,' SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF CALL FI636 (DATA,MSGA,KBMS, * REFNCE,KPTR,KPDS,KGDS) GO TO 900 END IF ELSE IF (KBDS(17).NE.0) THEN C PRINT *,' MATRIX OF VALS EACH PT' IF (KBDS(14).EQ.0) THEN C PRINT *,' NO SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (KBDS(14).NE.0) THEN C PRINT *,' SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF END IF END IF END IF END IF ELSE IF (IAND(KPTR(14),8).NE.0) THEN C PRINT *,'PROCESSING SPHERICAL HARMONIC COEFFICIENTS' IF (IAND(KPTR(14),4).EQ.0) THEN C PRINT *,' WITH SIMPLE PACKING' IF (IAND(KPTR(14),1).EQ.0) THEN C PRINT *,' WITH NO ADDITIONAL FLAGS' GO TO 5000 ELSE IF (IAND(KPTR(14),1).NE.0) THEN C PRINT *,' WITH ADDITIONAL FLAGS' IF (KBDS(17).EQ.0) THEN C PRINT *,' SINGLE DATUM EACH GRID PT' IF (KBDS(14).EQ.0) THEN C PRINT *,' NO SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (KBDS(14).NE.0) THEN C PRINT *,' SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF END IF ELSE IF (KBDS(17).NE.0) THEN C PRINT *,' MATRIX OF VALS EACH PT' IF (KBDS(14).EQ.0) THEN C PRINT *,' NO SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (KBDS(14).NE.0) THEN C PRINT *,' SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF END IF END IF END IF ELSE IF (IAND(KPTR(14),4).NE.0) THEN C COMPLEX/SECOND ORDER PACKING C PRINT *,' WITH COMPLEX/SECOND ORDER PACKING' IF (IAND(KPTR(14),1).EQ.0) THEN C PRINT *,' WITH NO ADDITIONAL FLAGS' ELSE IF (IAND(KPTR(14),1).NE.0) THEN C PRINT *,' WITH ADDITIONAL FLAGS' IF (KBDS(17).EQ.0) THEN C PRINT *,' SINGLE DATUM EACH GRID PT' IF (KBDS(14).EQ.0) THEN C PRINT *,' NO SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (KBDS(14).NE.0) THEN C PRINT *,' SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF END IF ELSE IF (KBDS(17).NE.0) THEN C PRINT *,' MATRIX OF VALS EACH PT' IF (KBDS(14).EQ.0) THEN C PRINT *,' NO SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF ELSE IF (KBDS(14).NE.0) THEN C PRINT *,' SEC BIT MAP' IF (KBDS(16).EQ.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES CONSTANT WIDTH' ELSE IF (KBDS(16).NE.0) THEN C PRINT *,' SECOND ORDER', C * ' VALUES DIFFERENT WIDTHS' END IF END IF END IF END IF END IF END IF C PRINT *,' NOT PROCESSED - NOT PROCESSED - NOT PROCESSED' KRET = 11 RETURN 4000 CONTINUE C **************************************************************** C C GRID POINT DATA, SIMPLE PACKING, FLOATING POINT, NO ADDN'L FLAGS C SCAL10 = 10.0 ** KPDS(22) SCAL10 = 1.0 / SCAL10 IF (KPDS(3).EQ.23.OR.KPDS(3).EQ.24.OR.KPDS(3).EQ.26. * OR.KPDS(3).EQ.63.OR.KPDS(3).EQ.64) THEN IF (KPDS(3).EQ.26) THEN KADD = 72 ELSE IF (KPDS(3).EQ.63.OR.KPDS(3).EQ.64) THEN KADD = 91 ELSE KADD = 37 END IF CALL GBYTES (MSGA,KSAVE,KPTR(8),KBITS,0,KNR) KPTR(8) = KPTR(8) + KBITS * KNR II = 1 KENTRY = KPTR(10) DO 4001 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 4001 CONTINUE DO 4002 I = 2, KADD DATA(I) = DATA(1) 4002 CONTINUE ELSE IF (KPDS(3).EQ.21.OR.KPDS(3).EQ.22.OR.KPDS(3).EQ.25. * OR.KPDS(3).EQ.61.OR.KPDS(3).EQ.62) THEN CALL GBYTES (MSGA,KSAVE,KPTR(8),KBITS,0,KNR) II = 1 KENTRY = KPTR(10) DO 4011 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 4011 CONTINUE IF (KPDS(3).EQ.25) THEN KADD = 71 ELSE IF (KPDS(3).EQ.61.OR.KPDS(3).EQ.62) THEN KADD = 90 ELSE KADD = 36 END IF LASTP = KENTRY - KADD DO 4012 I = LASTP+1, KENTRY DATA(I) = DATA(LASTP) 4012 CONTINUE ELSE CALL GBYTES (MSGA,KSAVE,KPTR(8),KBITS,0,KNR) II = 1 KENTRY = KPTR(10) DO 500 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 500 CONTINUE END IF GO TO 900 C ------------- PROCESS SPHERICAL HARMONIC COEFFICIENTS, C SIMPLE PACKING, FLOATING POINT, NO ADDN'L FLAGS 5000 CONTINUE C PRINT *,'CHECK POINT SPECTRAL COEFF' KPTR(8) = IBYT12 CALL GBYTE (MSGA,KKK,KPTR(8),32) KPTR(8) = KPTR(8) + 32 C C THE NEXT CODE WILL CONVERT THE IBM370 FOATING POINT C TO THE FLOATING POINT USED ON YOUR MACHINE. C C 1ST TEST TO SEE IN ON 32 OR 64 BIT WORD MACHINE C LW = 4 OR 8; IF 8 MAY BE A CRAY C CALL W3FI01(LW) IF (LW.EQ.4) THEN CALL GBYTE (KK,JSGN,0,1) CALL GBYTE (KK,JEXP,1,7) CALL GBYTE (KK,IFR,8,24) ELSE CALL GBYTE (KK,JSGN,32,1) CALL GBYTE (KK,JEXP,33,7) CALL GBYTE (KK,IFR,40,24) ENDIF C IF (IFR.EQ.0) THEN REALKK = 0.0 ELSE IF (JEXP.EQ.0.AND.IFR.EQ.0) THEN REALKK = 0.0 ELSE REALKK = FLOAT(IFR) * 16.0 ** (JEXP - 64 - 6) IF (JSGN.NE.0) REALKK = -REALKK END IF DATA(1) = REALKK CALL GBYTES (MSGA,KSAVE,KPTR(8),KBITS,0,KNR) C -------------- DO 6000 I = 1, KENTRY DATA(I+1) = REFNCE + FLOAT(KSAVE(I)) * SCALE 6000 CONTINUE 900 CONTINUE C PRINT *,'EXIT FI635' RETURN END SUBROUTINE FI636 (DATA,MSGA,KBMS,REFNCE,KPTR,KPDS,KGDS) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: FI636 PROCESS SECOND ORDER PACKING C PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 92-09-22 C C ABSTRACT: PROCESS SECOND ORDER PACKING FROM THE BINARY DATA SECTION C (BDS) FOR SINGLE DATA ITEMS GRID POINT DATA C C PROGRAM HISTORY LOG: C 93-06-08 CAVANAUGH C 93-12-15 CAVANAUGH MODIFIED SECOND ORDER POINTERS TO FIRST ORDER C VALUES AND SECOND ORDER VALUES CORRECTLY. C 95-04-26 R.E.JONES FI636 CORECTION FOR 2ND ORDER COMPLEX C UNPACKING. C 95-10-31 IREDELL REMOVED SAVES AND PRINTS C C USAGE: CALL FI636 (DATA,MSGA,KBMS,REFNCE,KPTR,KPDS,KGDS) C INPUT ARGUMENT LIST: C C MSGA - ARRAY CONTAINING GRIB MESSAGE C REFNCE - REFERENCE VALUE C KPTR - WORK ARRAY C C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) C DATA - LOCATION OF OUTPUT ARRAY C WORKING ARRAY C KBDS(1) - N1 C KBDS(2) - N2 C KBDS(3) - P1 C KBDS(4) - P2 C KBDS(5) - BIT POINTER TO 2ND ORDER WIDTHS C KBDS(6) - " " " " " BIT MAPS C KBDS(7) - " " " FIRST ORDER VALUES C KBDS(8) - " " " SECOND ORDER VALUES C KBDS(9) - " " START OF BDS C KBDS(10) - " " MAIN BIT MAP C KBDS(11) - BINARY SCALING C KBDS(12) - DECIMAL SCALING C KBDS(13) - BIT WIDTH OF FIRST ORDER VALUES C KBDS(14) - BIT MAP FLAG C 0 = NO SECOND ORDER BIT MAP C 1 = SECOND ORDER BIT MAP PRESENT C KBDS(15) - SECOND ORDER BIT WIDTH C KBDS(16) - CONSTANT / DIFFERENT WIDTHS C 0 = CONSTANT WIDTHS C 1 = DIFFERENT WIDTHS C KBDS(17) - SINGLE DATUM / MATRIX C 0 = SINGLE DATUM AT EACH GRID POINT C 1 = MATRIX OF VALUES AT EACH GRID POINT C (18-20)- UNUSED C C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT. C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C MACHINE: HDS, CRAY C C$$$ REAL DATA(*) REAL REFN REAL REFNCE C INTEGER KBDS(20) INTEGER KPTR(*) INTEGER JREF,BMAP2(12500) INTEGER I,IBDS INTEGER KBIT,IFOVAL,ISOVAL INTEGER KPDS(*),KGDS(*) C LOGICAL*1 KBMS(*) C CHARACTER*1 MSGA(*) C EQUIVALENCE (JREF,REFN) C ******************* SETUP ****************************** C PRINT *,'ENTER FI636' C START OF BMS (BIT POINTER) DO I = 1,20 KBDS(I) = 0 END DO C BYTE START OF BDS IBDS = KPTR(2) + KPTR(3) + KPTR(4) + KPTR(5) C PRINT *,'KPTR(2-5) ',KPTR(2),KPTR(3),KPTR(4),KPTR(5) C BIT START OF BDS JPTR = IBDS * 8 C PRINT *,'JPTR ',JPTR KBDS(9) = JPTR C PRINT *,'START OF BDS ',KBDS(9) C BINARY SCALE VALUE BDS BYTES 5-6 CALL GBYTE (MSGA,ISIGN,JPTR+32,1) CALL GBYTE (MSGA,KBDS(11),JPTR+33,15) IF (ISIGN.GT.0) THEN KBDS(11) = - KBDS(11) END IF C PRINT *,'BINARY SCALE VALUE =',KBDS(11) C EXTRACT REFERENCE VALUE CALL GBYTE(MSGA,JREF,JPTR+48,32) C PRINT *,'DECODED REFERENCE VALUE =',REFN,REFNCE C F O BIT WIDTH CALL GBYTE(MSGA,KBDS(13),JPTR+80,8) JPTR = JPTR + 88 C AT START OF BDS BYTE 12 C EXTRACT N1 CALL GBYTE (MSGA,KBDS(1),JPTR,16) C PRINT *,'N1 = ',KBDS(1) JPTR = JPTR + 16 C EXTENDED FLAGS CALL GBYTE (MSGA,KFLAG,JPTR,8) C ISOLATE BIT MAP FLAG IF (IAND(KFLAG,32).NE.0) THEN KBDS(14) = 1 ELSE KBDS(14) = 0 END IF IF (IAND(KFLAG,16).NE.0) THEN KBDS(16) = 1 ELSE KBDS(16) = 0 END IF IF (IAND(KFLAG,64).NE.0) THEN KBDS(17) = 1 ELSE KBDS(17) = 0 END IF JPTR = JPTR + 8 C EXTRACT N2 CALL GBYTE (MSGA,KBDS(2),JPTR,16) C PRINT *,'N2 = ',KBDS(2) JPTR = JPTR + 16 C EXTRACT P1 CALL GBYTE (MSGA,KBDS(3),JPTR,16) C PRINT *,'P1 = ',KBDS(3) JPTR = JPTR + 16 C EXTRACT P2 CALL GBYTE (MSGA,KBDS(4),JPTR,16) C PRINT *,'P2 = ',KBDS(4) JPTR = JPTR + 16 C SKIP RESERVED BYTE JPTR = JPTR + 8 C START OF SECOND ORDER BIT WIDTHS KBDS(5) = JPTR C COMPUTE START OF SECONDARY BIT MAP IF (KBDS(14).NE.0) THEN C FOR INCLUDED SECONDARY BIT MAP JPTR = JPTR + (KBDS(3) * 8) KBDS(6) = JPTR ELSE C FOR CONSTRUCTED SECONDARY BIT MAP KBDS(6) = 0 END IF C CREATE POINTER TO START OF FIRST ORDER VALUES KBDS(7) = KBDS(9) + KBDS(1) * 8 - 8 C PRINT *,'BIT POINTER TO START OF FOVALS',KBDS(7) C CREATE POINTER TO START OF SECOND ORDER VALUES KBDS(8) = KBDS(9) + KBDS(2) * 8 - 8 C PRINT *,'BIT POINTER TO START OF SOVALS',KBDS(8) C PRINT *,'KBDS( 1) - N1 ',KBDS( 1) C PRINT *,'KBDS( 2) - N2 ',KBDS( 2) C PRINT *,'KBDS( 3) - P1 ',KBDS( 3) C PRINT *,'KBDS( 4) - P2 ',KBDS( 4) C PRINT *,'KBDS( 5) - BIT PTR - 2ND ORDER WIDTHS ',KBDS( 5) C PRINT *,'KBDS( 6) - " " " " BIT MAPS ',KBDS( 6) C PRINT *,'KBDS( 7) - " " F O VALS ',KBDS( 7) C PRINT *,'KBDS( 8) - " " S O VALS ',KBDS( 8) C PRINT *,'KBDS( 9) - " " START OF BDS ',KBDS( 9) C PRINT *,'KBDS(10) - " " MAIN BIT MAP ',KBDS(10) C PRINT *,'KBDS(11) - BINARY SCALING ',KBDS(11) C PRINT *,'KPDS(22) - DECIMAL SCALING ',KPDS(22) C PRINT *,'KBDS(13) - FO BIT WIDTH ',KBDS(13) C PRINT *,'KBDS(14) - 2ND ORDER BIT MAP FLAG ',KBDS(14) C PRINT *,'KBDS(15) - 2ND ORDER BIT WIDTH ',KBDS(15) C PRINT *,'KBDS(16) - CONSTANT/DIFFERENT WIDTHS ',KBDS(16) C PRINT *,'KBDS(17) - SINGLE DATUM/MATRIX ',KBDS(17) C PRINT *,'REFNCE VAL ',REFNCE C ************************* PROCESS DATA ********************** IJ = 0 C ======================================================== IF (KBDS(14).EQ.0) THEN C NO BIT MAP, MUST CONSTRUCT ONE IF (KGDS(2).EQ.65535) THEN IF (KGDS(20).EQ.255) THEN C PRINT *,'CANNOT BE USED HERE' ELSE C POINT TO PL LP = KPTR(9) + KPTR(2)*8 + KPTR(3)*8 + KGDS(20)*8 - 8 C PRINT *,'LP = ',LP JT = 0 DO 2000 JZ = 1, KGDS(3) C GET NUMBER IN CURRENT ROW CALL GBYTE (MSGA,NUMBER,LP,16) C INCREMENT TO NEXT ROW NUMBER LP = LP + 16 C PRINT *,'NUMBER IN ROW',JZ,' = ',NUMBER DO 1500 JQ = 1, NUMBER IF (JQ.EQ.1) THEN CALL SBYTE (BMAP2,1,JT,1) ELSE CALL SBYTE (BMAP2,0,JT,1) END IF JT = JT + 1 1500 CONTINUE 2000 CONTINUE END IF ELSE IF (IAND(KGDS(11),32).EQ.0) THEN C ROW BY ROW C PRINT *,' ROW BY ROW' KOUT = KGDS(3) KIN = KGDS(2) ELSE C COL BY COL C PRINT *,' COL BY COL' KIN = KGDS(3) KOUT = KGDS(2) END IF C PRINT *,'KIN=',KIN,' KOUT= ',KOUT DO 200 I = 1, KOUT DO 150 J = 1, KIN IF (J.EQ.1) THEN CALL SBYTE (BMAP2,1,IJ,1) ELSE CALL SBYTE (BMAP2,0,IJ,1) END IF IJ = IJ + 1 150 CONTINUE 200 CONTINUE END IF END IF C ======================================================== C PRINT 99,(BMAP2(J),J=1,110) C99 FORMAT ( 10(1X,Z8.8)) C CALL BINARY (BMAP2,2) C FOR EACH GRID POINT ENTRY C SCALE2 = 2.0**KBDS(11) SCAL10 = 10.0**KPDS(22) C PRINT *,'SCALE VALUES - ',SCALE2,SCAL10 DO 1000 I = 1, KPTR(10) C GET NEXT MASTER BIT MAP BIT POSITION C IF NEXT MASTER BIT MAP BIT POSITION IS 'ON' (1) IF (KBMS(I)) THEN C WRITE(6,900)I,KBMS(I) C 900 FORMAT (1X,I4,3X,14HMAIN BIT IS ON,3X,L4) IF (KBDS(14).NE.0) THEN CALL GBYTE (MSGA,KBIT,KBDS(6),1) ELSE CALL GBYTE (BMAP2,KBIT,KBDS(6),1) END IF C PRINT *,'KBDS(6) =',KBDS(6),' KBIT =',KBIT KBDS(6) = KBDS(6) + 1 IF (KBIT.NE.0) THEN C PRINT *,' SOB ON' C GET NEXT FIRST ORDER PACKED VALUE CALL GBYTE (MSGA,IFOVAL,KBDS(7),KBDS(13)) KBDS(7) = KBDS(7) + KBDS(13) C PRINT *,'FOVAL =',IFOVAL C GET SECOND ORDER BIT WIDTH CALL GBYTE (MSGA,KBDS(15),KBDS(5),8) KBDS(5) = KBDS(5) + 8 C PRINT *,KBDS(7)-KBDS(13),' FOVAL =',IFOVAL,' KBDS(5)=', C * ,KBDS(5), 'ISOWID =',KBDS(15) ELSE C PRINT *,' SOB NOT ON' END IF ISOVAL = 0 IF (KBDS(15).EQ.0) THEN C IF SECOND ORDER BIT WIDTH = 0 C THEN SECOND ORDER VALUE IS 0 C SO CALCULATE DATA VALUE FOR THIS POINT C DATA(I) = (REFNCE + (FLOAT(IFOVAL) * SCALE2)) / SCAL10 ELSE CALL GBYTE (MSGA,ISOVAL,KBDS(8),KBDS(15)) KBDS(8) = KBDS(8) + KBDS(15) END IF DATA(I) = (REFNCE + (FLOAT(IFOVAL + ISOVAL) * * SCALE2)) / SCAL10 C PRINT *,I,DATA(I),REFNCE,IFOVAL,ISOVAL,SCALE2,SCAL10 ELSE C WRITE(6,901) I,KBMS(I) C 901 FORMAT (1X,I4,3X,15HMAIN BIT NOT ON,3X,L4) DATA(I) = 0.0 END IF C PRINT *,I,DATA(I),IFOVAL,ISOVAL,KBDS(5),KBDS(15) 1000 CONTINUE C ************************************************************** C PRINT *,'EXIT FI636' RETURN END SUBROUTINE FI637(*,J,KPDS,KGDS,KRET) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: FI637 GRIB GRID/SIZE TEST C PRGMMR: CAVANAUGH ORG: W/NMC42 DATE: 91-09-13 C C ABSTRACT: TO TEST WHEN GDS IS AVAILABLE TO SEE IF SIZE MISMATCH C ON EXISTING GRIDS (BY CENTER) IS INDICATED C C PROGRAM HISTORY LOG: C 91-09-13 CAVANAUGH C 95-10-31 IREDELL REMOVED SAVES AND PRINTS C C USAGE: CALL FI637(*,J,KPDS,KGDS,KRET) C INPUT ARGUMENT LIST: C J - SIZE FOR INDICATED GRID C KPDS - C KGDS - C C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) C KRET - ERROR RETURN C C REMARKS: C KRET - C = 9 - GDS INDICATES SIZE MISMATCH WITH STD GRID C C SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT. C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C MACHINE: HDS C C$$$ INTEGER KPDS(*) INTEGER KGDS(*) INTEGER J INTEGER I C --------------------------------------- C --------------------------------------- C IF GDS NOT INDICATED, RETURN C ---------------------------------------- IF (IAND(KPDS(4),128).EQ.0) RETURN C --------------------------------------- C GDS IS INDICATED, PROCEED WITH TESTING C --------------------------------------- IF (KGDS(2).EQ.65535) THEN RETURN END IF I = KGDS(2) * KGDS(3) C --------------------------------------- C INTERNATIONAL SET C --------------------------------------- IF (KPDS(3).GE.21.AND.KPDS(3).LE.26) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.37.AND.KPDS(3).LE.44) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).EQ.50) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.61.AND.KPDS(3).LE.64) THEN IF (I.NE.J) THEN RETURN 1 END IF C --------------------------------------- C TEST ECMWF CONTENT C --------------------------------------- ELSE IF (KPDS(1).EQ.98) THEN KRET = 9 IF (KPDS(3).GE.1.AND.KPDS(3).LE.16) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE KRET = 5 RETURN 1 END IF C --------------------------------------- C U.K. MET OFFICE, BRACKNELL C --------------------------------------- ELSE IF (KPDS(1).EQ.74) THEN KRET = 9 IF (KPDS(3).GE.25.AND.KPDS(3).LE.26) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE KRET = 5 RETURN 1 END IF C --------------------------------------- C CANADA C --------------------------------------- ELSE IF (KPDS(1).EQ.54) THEN C PRINT *,' NO CURRENT LISTING OF CANADIAN GRIDS' RETURN 1 C --------------------------------------- C JAPAN METEOROLOGICAL AGENCY C --------------------------------------- ELSE IF (KPDS(1).EQ.34) THEN C PRINT *,' NO CURRENT LISTING OF JMA GRIDS' RETURN 1 C --------------------------------------- C NAVY - FNOC C --------------------------------------- ELSE IF (KPDS(1).EQ.58) THEN IF (KPDS(3).GE.37.AND.KPDS(3).LE.44) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.220.AND.KPDS(3).LE.221) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).EQ.223) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE KRET = 5 RETURN 1 END IF C --------------------------------------- C U.S. GRIDS C --------------------------------------- ELSE IF (KPDS(1).EQ.7) THEN KRET = 9 IF (KPDS(3).GE.1.AND.KPDS(3).LE.4) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).EQ.5.OR.KPDS(3).EQ.6) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).EQ.8) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.27.AND.KPDS(3).LE.30) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.33.AND.KPDS(3).LE.34) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.37.AND.KPDS(3).LE.44) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).EQ.53) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.55.AND.KPDS(3).LE.56) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.67.AND.KPDS(3).LE.77) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.85.AND.KPDS(3).LE.86) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).EQ.87) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.90.AND.KPDS(3).LE.98) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).EQ.100.OR.KPDS(3).EQ.101) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.103.AND.KPDS(3).LE.107) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).EQ.126) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).EQ.196) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE IF (KPDS(3).GE.201.AND.KPDS(3).LE.216) THEN IF (I.NE.J) THEN RETURN 1 END IF ELSE KRET = 5 RETURN 1 END IF ELSE KRET = 10 RETURN 1 END IF C ------------------------------------ C NORMAL EXIT C ------------------------------------ KRET = 0 RETURN END