| blob | 40c3e1352e10387d0458a68b73a74ad814f07b6f |
1 ! !
2 !*****************************************************************************!
3 ! !
4 ! This is a package of subroutines to read GRIB-formatted data. It is still !
5 ! under continuous development. It will not be able to read every GRIB dataset!
6 ! you give it, but it will read a good many. !
7 ! !
8 ! Kevin W. Manning !
9 ! NCAR/MMM !
10 ! Summer 1998, and continuing !
11 ! SDG !
12 ! !
13 !*****************************************************************************!
14 ! !
15 ! The main user interfaces are: !
16 ! !
17 ! SUBROUTINE GRIBGET(NUNIT, IERR) !
18 ! Read a single GRIB record from UNIX file-descriptor NUNIT into array !
19 ! GREC. No unpacking of any header or data values is performed. !
20 ! !
21 ! SUBROUTINE GRIBREAD(NUNIT, DATA, NDATA, IERR) !
22 ! Read a single GRIB record from UNIX file-descriptor NUNIT, and unpack !
23 ! all header and data values into the appropriate arrays. !
24 ! !
25 ! SUBROUTINE GRIBHEADER !
26 ! Unpack the header of a GRIB record !
27 ! !
28 ! SUBROUTINE GRIBDATA(DATARRAY, NDAT) !
29 ! Unpack the data in a GRIB record into array DATARRAY !
30 ! !
31 ! SUBROUTINE GRIBPRINT(ISEC) !
32 ! Print the header information from GRIB section ISEC. !
33 ! !
34 ! SUBROUTINE GET_SEC1(KSEC1) !
35 ! Return the header information from Section 1. !
36 ! !
37 ! SUBROUTINE GET_SEC2(KSEC2) !
38 ! Return the header information from Section 2. !
39 ! !
40 ! SUBROUTINE GET_GRIDINFO(IGINFO, GINFO) !
41 ! Return the grid information of the previously-unpacked GRIB header. !
42 ! !
43 ! !
44 !*****************************************************************************!
45 ! !
46 ! !
47 ! The following arrays have meanings as follows: !
48 ! !
49 ! !
50 ! SEC0: GRIB Header Section 0 information !
51 ! !
52 ! 1 : Length of a complete GRIB record !
53 ! 2 : GRIB Edition number !
54 ! !
55 ! !
56 ! SEC1: GRIB Header Section 1 information !
57 ! !
58 ! 1 : Length of GRIB section 1 (bytes) !
59 ! 2 : Parameter Table Version number ???? !
60 ! 3 : Center Identifier ???? !
61 ! 4 : Process Identifier ???? !
62 ! 5 : Grid ID number for pre-specified grids. !
63 ! 6 : Binary bitmap flag: !
64 ! 7 : Parameter ID Number (ON388 Table 2) !
65 ! 8 : Level type (ON388 Table 3) !
66 ! 9 : Level value, or top value of a layer !
67 ! 10 : Bottom value of a layer ( 0 if NA ??) !
68 ! 11 : Year (00-99) !
69 ! 12 : Month (01-12) !
70 ! 13 : Day of the month (01-31) !
71 ! 14 : Hour (00-23) !
72 ! 15 : Minute (00-59) !
73 ! 16 : Forecast time unit: (ON388 Table 4) !
74 ! 17 : Time period 1: !
75 ! 18 : Time period 2: !
76 ! 19 : Time range indicator (ON833 Table 5) !
77 ! 20 : Number of ?? in an average ?? !
78 ! 21 : Number of ?? missing from average ?? !
79 ! 22 : Century (Years 1999 and 2000 are century 20, 2001 is century 21)!
80 ! 23 : Sub-center identifier ?? !
81 ! 24 : Decimal scale factor for ??? !
82 ! !
83 ! !
84 ! !
85 ! !
86 ! !
87 ! SEC2: GRIB Header Section 2 information !
88 ! !
89 ! 1 : Length of GRIB Section 2 !
90 ! 2 : Number of vertical-coordinate parameters ??? !
91 ! 3 : Starting-point of the list of vertical-coordinate parameters ?? !
92 ! 4 : Data-representation type (i.e., grid type) Table ??? !
93 ! : 0 = ?? !
94 ! : 3 = Lambert-conformal grid. !
95 ! : 5 = Polar-stereographic grid. !
96 ! !
97 ! if (SEC2(4) == 0) then LATITUDE/LONGITUDE GRID !
98 ! !
99 ! INFOGRID/GRIDINFO: !
100 ! !
101 ! 1 : I Dimension of the grid !
102 ! 2 : J Dimension of the grid !
103 ! 3 : Starting Latitude of the grid. !
104 ! 4 : Starting Longitude of the grid. !
105 ! 5 : Resolution and component flags. !
106 ! 6 : Ending latitude of the grid. !
107 ! 7 : Ending longitude of the grid. !
108 ! 8 : Longitudinal increment. !
109 ! 9 : Latitudinal incriment. !
110 ! 10 : Scanning mode (bit 3 from Table 8) !
111 ! 21 : Iscan sign (+1/-1) (bit 1 from Table 8) !
112 ! 22 : Jscan sign (+1/-1) (bit 2 from Table 8) !
113 ! !
114 ! !
115 ! elseif (SEC2(4) == 3) then LAMBERT CONFORMAL GRID !
116 ! !
117 ! INFOGRID/GRIDINFO: !
118 ! !
119 ! 1 : I Dimension of the grid !
120 ! 2 : J Dimension of the grid !
121 ! 3 : Starting Latitude of the grid. !
122 ! 4 : Starting Longitude of the grid. !
123 ! 5 : Resolution and component flags. !
124 ! 6 : Center longitude of the projection. !
125 ! 7 : Grid-spacing in the I direction !
126 ! 8 : Grid-spacing in the J direction !
127 ! 9 : Projection center !
128 ! 10 : Scanning mode (bit 3 from Table 8) !
129 ! 11 : First TRUELAT value. !
130 ! 12 : Second TRUELAT value. !
131 ! 13 : Latitude of the southern pole ?? !
132 ! 14 : Longitude of the southern pole ?? !
133 ! 21 : Iscan sign (+1/-1) (bit 1 from Table 8) !
134 ! 22 : Jscan sign (+1/-1) (bit 2 from Table 8) !
135 ! !
136 ! if (SEC2(4) == 4) then GAUSSIAN GRID !
137 ! !
138 ! INFOGRID/GRIDINFO: !
139 ! !
140 ! 1 : I Dimension of the grid !
141 ! 2 : J Dimension of the grid !
142 ! 3 : Starting Latitude of the grid. !
143 ! 4 : Starting Longitude of the grid. !
144 ! 5 : Resolution and component flags. !
145 ! 6 : Ending latitude of the grid. !
146 ! 7 : Ending longitude of the grid. !
147 ! 8 : Longitudinal increment. !
148 ! 9 : Number of latitude circles between pole and equator !
149 ! 10 : Scanning mode (bit 3 from Table 8) !
150 ! 17 : Original (stored) ending latitude !
151 ! 18 : Original (stored) starting latitude !
152 ! 19 : Approximate delta-latitude !
153 ! 21 : Iscan sign (+1/-1) (bit 1 from Table 8) !
154 ! 22 : Jscan sign (+1/-1) (bit 2 from Table 8) !
155 ! !
156 ! !
157 ! elseif (SEC2(4) == 5) then POLAR STEREOGRAPHIC GRID !
158 ! !
159 ! INFOGRID/GRIDINFO: !
160 ! !
161 ! 1 : I Dimension of the grid !
162 ! 2 : J Dimension of the grid !
163 ! 3 : Starting Latitude of the grid. !
164 ! 4 : Starting Longitude of the grid. !
165 ! 5 : Resolution and component flags. !
166 ! 6 : Center longitude of the projection. !
167 ! 7 : Grid-spacing in the I direction !
168 ! 8 : Grid-spacing in the J direction !
169 ! 9 : Projection center !
170 ! 10 : Scanning mode (bit 3 from Table 8) !
171 ! 21 : Iscan sign (+1/-1) (bit 1 from Table 8) !
172 ! 22 : Jscan sign (+1/-1) (bit 2 from Table 8) !
173 ! !
174 ! elseif (SEC2(4) == 50) then SPHERICAL HARMONIC COEFFICIENTS !
175 ! !
176 ! INFOGRID/GRIDINFO: !
177 ! !
178 ! 1 : J-pentagonal resolution parameter !
179 ! 2 : K-pentagonal resolution parameter !
180 ! 3 : M-pentagonal resolution parameter !
181 ! 4 : Spectral representation type (ON388 Table 9) !
182 ! 5 : Coefficient storage mode (ON388 Table 10) !
183 ! !
184 ! elseif (SEC2(4) == ?) then ?? !
185 ! !
186 ! !
187 ! SEC3: GRIB Header Section 3 information: !
188 ! SEC4: GRIB Header Section 4 information: !
189 !
190 !
191 !
192 module module_grib
193 !
194 ! Machine wordsize must be known for the various unpacking routines to work.
195 ! Machine wordsize is set through CPP Directives.
196 ! Use options -DBIT32 (for 32-bit word-size) or -DBIT64 (for 64-bit wordsize)
197 ! for the CPP pass of the compiler.
198 !
199 #if defined (BIT32)
200 integer, parameter :: MWSIZE = 32 ! Machine word size in bits
201 #elif defined (BIT64)
202 integer, parameter :: MWSIZE = 64 ! Machine word size in bits
203 #endif
206 ! Array GREC holds a single packed GRIB record (header and all).
207 ! Array BITMAP holds the bitmap (if a bitmap was used).
208 !
209 ! For some reason, the cray does not like grec to be allocatable.
210 !
211 #if defined (CRAY)
212 integer, dimension(100000) :: grec
213 integer, dimension(100000) :: bitmap
214 #else
215 integer, allocatable, save, dimension(:) :: grec
216 integer, allocatable, save, dimension(:) :: bitmap
217 #endif
219 ! SEC0 holds the Section 0 header information
220 ! SEC1 holds the Section 1 header information
221 ! SEC2 holds the Section 2 header information
222 ! SEC3 holds the Section 3 header information
223 ! SEC4 holds the Section 4 header information
224 ! XEC4 holds floating-point Section 4 header information
226 integer, dimension(2) :: sec0
227 integer, dimension(100) :: sec1
228 integer, dimension(10) :: sec2
229 integer, dimension(10) :: sec3
230 integer, dimension(10) :: sec4
231 real, dimension(1) :: xec4
233 integer :: sevencount = 0
235 ! INFOGRID holds integer values defining the grid.
236 ! GRIDINFO holds floating-point values definint the grid
238 integer, dimension(40) :: infogrid
239 real, dimension(40) :: gridinfo
241 integer :: ied
242 real, parameter :: pi = 3.1415926534
243 real, parameter :: degran = pi/180.
244 real, parameter :: raddeg = 1./degran
246 real :: glat1, glon1, gclon, gtrue1, gtrue2, grrth, gx1, gy1, gkappa
248 contains
249 !
250 !=============================================================================!
251 !=============================================================================!
252 !=============================================================================!
253 !
254 integer function gribsize(trec, ilen, ierr)
255 !-----------------------------------------------------------------------------!
256 ! Return the size of a single GRIB record. !
257 ! !
258 ! Input: !
259 ! TREC: At least a portion of the complete GRIB record. !
260 ! ILEN: The size of array TREC. !
261 ! !
262 ! Output !
263 ! GRIBSIZE: The size of the full GRIB record !
264 ! IERR : 0= no errors, 1 = read error
265 ! !
266 ! Side Effects: !
267 ! * Module variable IED is set to the GRIB Edition number. !
268 ! * STOP, if not GRIB Edition 0 or 1 !
269 ! !
270 ! Externals !
271 ! GBYTE !
272 ! !
273 !-----------------------------------------------------------------------------!
274 implicit none
275 integer :: ilen
276 integer, dimension(ilen) :: trec
277 integer :: isz0 = 32
278 integer :: isz1 = 0
279 integer :: isz2 = 0
280 integer :: isz3 = 0
281 integer :: isz4 = 0
282 integer :: isz5 = 32
283 integer :: iflag
284 integer :: ierr
285 character :: pname*132
287 ierr = 0
288 ! Unpack the GRIB Edition number, located in the eighth byte (bits 57-64) !
289 ! of array TREC. !
291 call gbyte_g1(trec, ied, 56, 8)
293 ! GRIB Edition 1 has the size of the whole GRIB record right up front. !
295 if (ied.eq.1) then
296 ! Grib1
297 ! Find the size of the whole GRIB record
298 call gbyte_g1(trec, gribsize, 32, 24)
300 ! GRIB Edition 0 does not include the total size, so we have to sum up !
301 ! the sizes of the individual sections !
303 elseif (ied.eq.0) then
304 ! Grib old edition
305 ! Find the size of section 1.
306 call gbyte_g1(trec, isz1, isz0, 24)
307 isz1 = isz1 * 8
308 call gbyte_g1(trec, iflag, isz0+56, 8)
309 if ((iflag.eq.128).or.(iflag.eq.192)) then ! section 2 is there
310 ! Find the size of section 2.
311 call gbyte_g1(trec, isz2, isz0+isz1, 24)
312 isz2 = isz2 * 8
313 endif
314 if ((iflag.eq.64).or.(iflag.eq.192)) then ! Section 3 is there
315 ! Find the size of section 3.
316 call gbyte_g1(trec, isz3, isz0+isz1+isz2, 24)
317 isz3 = isz3 * 8
318 endif
319 ! Find the size of section 4.
320 call gbyte_g1(trec, isz4, isz0+isz1+isz2+isz3, 24)
321 isz4 = isz4 * 8
323 ! Total the sizes of sections 0 through 5.
324 gribsize = (isz0+isz1+isz2+isz3+isz4+isz5) / 8
326 elseif (ied.eq.2) then
327 ! Grib2
328 CALL getarg ( 0 , pname )
329 write(*,'("*** stopping in gribcode ***\n")')
330 write(*,'("\tI was expecting a Grib1 file, but this is a Grib2 file.")')
331 if ( index(pname,'ungrib.exe') .ne. 0 ) then
332 write(*,'("\tIt is possible this is because your GRIBFILE.XXX files")')
333 write(*,'("\tare not all of the same type.")')
334 write(*,'("\tWPS can handle both file types, but a separate ungrib")')
335 write(*,'("\tjob must be run for each Grib type.\n")')
336 else
337 write(*,'("\tUse g2print on Grib2 files\n")')
338 endif
339 stop 'gribsize in gribcode'
340 else
341 write(*,'("Error trying to read grib edition number in gribsize.")')
342 write(*,'("Possible corrupt grib file.")')
343 write(6,*) 'Incorrect edition number = ',ied
344 write(6,*) 'Skipping the rest of the file and continuing.'
345 ierr = 1
346 endif
347 end function gribsize
348 !
349 !=============================================================================!
350 !=============================================================================!
351 !=============================================================================!
352 !
353 subroutine findgrib(nunit, isize, ierr)
355 !-----------------------------------------------------------------------------!
356 ! !
357 ! Find the string "GRIB", which starts off a GRIB record. !
358 ! !
359 ! Input: !
360 ! NUNIT: The C unit to read from. This should already be opened. !
361 ! !
362 ! Output: !
363 ! ISIZE: The size in bytes of one complete GRIB Record !
364 ! IERR: Error flag, !
365 ! 0 : No error or end-of-file on reading !
366 ! 1 : Hit the end of file !
367 ! 2 : Error on reading !
368 ! !
369 ! Side effects: !
370 ! * The pointer to C unit NUNIT is set to the beginning of the next !
371 ! GRIB record. !
372 ! * The first time FINDGRIB is called, the integer GTEST is set to !
373 ! a value equivalent to the string 'GRIB' !
374 ! !
375 ! Modules: !
376 ! MODULE_GRIB !
377 ! !
378 ! Externals: !
379 ! BN_READ !
380 ! BN_SEEK !
381 ! GRIBSIZE !
382 ! !
383 !-----------------------------------------------------------------------------!
384 implicit none
385 integer, intent(in) :: nunit
386 integer, intent(out) :: isize
387 integer, intent(out) :: ierr
389 integer, parameter :: LENTMP=100
390 integer, dimension(lentmp) :: trec
392 integer :: isz, itest, icnt
394 integer, save :: gtest = 0
396 ! Set the integer variable GTEST to hold the integer representation of the
397 ! character string 'GRIB'. This integer variable is later compared to
398 ! integers we read from the GRIB file, to find the beginning of a GRIB record.
400 if (gtest.eq.0) then
401 if (mwsize.eq.32) then
402 gtest = transfer('GRIB', gtest)
403 elseif(mwsize.eq.64) then
404 call gbyte_g1(char(0)//char(0)//char(0)//char(0)//'GRIB', gtest, 0, mwsize)
405 endif
406 endif
407 ierr = 0
408 icnt = 0
410 LOOP : DO
411 ! Read LENTMP bytes into holding array TREC.
412 call bn_read(nunit, trec, lentmp, isz, ierr, 0)
413 if (ierr.eq.1) then
414 return
415 elseif (ierr.eq.2) then
416 write(*,'("Error reading GRIB: IERR = ", I2)') ierr
417 return
418 endif
419 ! Reposition the file pointer back to where we started.
420 call bn_seek(nunit, -isz, 0, 0)
422 ! Compare the first four bytes of TREC with the string 'GRIB' stored in
423 ! integer variable GTEST.
424 if (mwsize.eq.32) then
425 if (trec(1) == gtest) exit LOOP
426 elseif (mwsize.eq.64) then
427 call gbyte_g1(trec, itest, 0, 32)
428 if (itest == gtest) exit LOOP
429 endif
431 ! Advance the file pointer one byte.
432 call bn_seek(nunit, 1, 0, 0)
433 icnt = icnt + 1
434 if ( icnt .gt. 100000) then ! stop if we cannot find the GRIB string
435 write(*,'("*** stopping in findgrib in gribcode ***\n")')
436 write(*,'("\tI could not find the GRIB string in the input file")')
437 write(*,'("\tafter testing the first 100,000 bytes.")')
438 write(*,'("\tThe file may be corrupt or it is not a GRIB file.")')
439 write(*,'("\tPerhaps a gzipped GRIB file or netcdf?\n")')
440 stop 'findgrib'
441 endif
443 ENDDO LOOP
445 !#if defined (DEC) || defined (ALPHA) || defined (alpha) || defined (LINUX)
446 #ifdef BYTESWAP
447 call swap4(trec, isz)
448 #endif
449 isize = gribsize(trec, isz, ierr)
451 end subroutine findgrib
452 !
453 !=============================================================================!
454 !=============================================================================!
455 !=============================================================================!
456 !
457 subroutine SGUP_NOBITMAP(datarray, ndat)
458 ! Simple grid-point unpacking
459 implicit none
461 integer :: ndat
462 real , dimension(ndat) :: datarray
463 integer, dimension(ndat) :: IX
464 real :: dfac, bfac
465 integer :: iskip
467 DFAC = 10.**(-sec1(24))
468 BFAC = 2.**sec4(7)
469 if (ied.eq.0) then
470 iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
471 elseif (ied.eq.1) then
472 iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
473 endif
474 ! sec4(8) is the number of bits used per datum value.
475 ! If sec4(8) = 255, assume they mean sec4(8) = 0
476 if (sec4(8) == 255) then
477 sec4(8) = 0
478 endif
479 ! If sec4(8) is 0, assume datarray is constant value of xec4(1)
481 if (sec4(8).eq.0) then
482 !!! HERE IS THE ORIGINAL NCAR CODE:
483 ! datarray = xec4(1)
484 !!! HERE IS WHAT FSL CHANGED IT TO:
485 datarray = DFAC*xec4(1)
486 !!! because even though it is a constant value
487 !!! you still need to scale by the decimal scale factor.
488 else
489 !!! FSL developers MOVED THE CALL TO gbytes FROM line 441 ABOVE
490 !!! BECAUSE IF sec4(8)==0 BEFORE gbytes IS CALLED, THE MASKS ARRAY
491 !!! IN gbytes WILL BE INDEXED OUT OF BOUNDS. C HARROP 9/16/04
492 call gbytes_g1(grec, IX, iskip, sec4(8), 0, ndat)
493 datarray = DFAC * (xec4(1) + (IX*BFAC))
494 endif
495 end subroutine SGUP_NOBITMAP
496 !
497 !=============================================================================!
498 !=============================================================================!
499 !=============================================================================!
500 !
502 subroutine SGUP_BITMAP(datarray, ndat)
503 ! Simple grid-point unpacking, with a bitmap.
504 implicit none
506 integer :: ndat ! Number of data points in the final grid.
507 real , dimension(ndat) :: datarray ! Array holding the final unpacked data.
508 real :: dfac, bfac
509 integer :: iskip, nbm, i, nn
511 integer, allocatable, dimension(:) :: bmdat
513 ! SEC4(1) : The number of bytes in the whole of GRIB Section 4.
514 ! SEC4(6) : The number of unused bits at the end of GRIB Section 4.
515 ! SEC4(8) : The number of bits per data value.
517 datarray = -1.E30
519 ! 1) There are fewer than NDAT data values, because a bitmap was used.
520 ! Compute the number of data values (NBM). There are 11 extra bytes
521 ! in the header section 4. NBM equals the total number of data bits (not
522 ! counting the header bits), minus the number of unused buts, and then
523 ! divided by the number of bits per data value.
525 ! Compute the parameters involved with packing
526 DFAC = 10.**(-sec1(24))
527 BFAC = 2.**sec4(7)
529 ! If sec4(8) is 0, assume datarray is constant value of xec4(1) scaled by DFAC
531 if (sec4(8).eq.0) then
532 where(bitmap(1:ndat).eq.1) datarray = xec4(1) * DFAC
533 return
534 endif
535 nbm = ((sec4(1)-11)*8-sec4(6))/sec4(8)
536 allocate(bmdat(nbm))
538 ! Set ISKIP to the beginning of the data.
539 if (ied.eq.0) then
540 iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
541 elseif (ied.eq.1) then
542 iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
543 endif
545 ! Read the data from the GREC array
546 call gbytes_g1(grec, bmdat, iskip, sec4(8), 0, nbm)
547 ! sec4(8) is the number of bits used per datum value.
548 ! If sec4(8) = 255, assume they mean sec4(8) = 0
549 if (sec4(8) == 255) sec4(8) = 0
551 ! Unpack the data according to packing parameters DFAC, BFAC, and XEC4(1),
552 ! and masked by the bitmap BITMAP.
553 nn = 0
554 do i = 1, ndat
555 if (bitmap(i).eq.1) then
556 nn = nn + 1
557 datarray(i) = DFAC * (xec4(1) + (bmdat(nn)*BFAC))
558 endif
559 enddo
561 ! Deallocate the scratch BMDAT array
562 deallocate(bmdat)
564 end subroutine SGUP_BITMAP
565 !
566 !=============================================================================!
567 !=============================================================================!
568 !=============================================================================!
569 !
571 subroutine CSHUP(pdata, ndat)
572 ! ECMWFs unpacking of ECMWFs Complex Spherical Harmonic packing
573 ! Adapted from ECMWFs GRIBEX package.
574 implicit none
576 integer :: ndat
577 real , dimension(ndat) :: pdata
578 integer, dimension(ndat+500) :: IX
580 integer :: iskip, isign
581 integer :: N1, IPOWER, J, K, M, nval
582 real :: zscale, zref
584 integer :: ic, jm, iuc, il2, inum, jn
585 integer :: inext, ilast, ioff, jrow, index, i, jcol
586 real :: bval
587 integer, allocatable, dimension(:) :: iexp, imant
588 real , dimension(0:400) :: factor
589 real :: power
590 integer :: N
592 index = -1
594 if (ied.eq.0) then
595 iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
596 elseif(ied.eq.1) then
597 iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
598 endif
600 call gbyte_g1(grec,N1,iskip,16)
601 iskip = iskip + 16
603 call gbyte_g1(grec,ipower,iskip,16)
604 iskip = iskip + 16
605 if (ipower.ge.32768) ipower = 32768-ipower
607 ! Unpack the resolution parameters for the initial (small) truncation:
608 call gbyte_g1(grec,J,iskip,8)
609 iskip = iskip + 8
610 call gbyte_g1(grec,K,iskip,8)
611 iskip = iskip + 8
612 call gbyte_g1(grec,M,iskip,8)
613 iskip = iskip + 8
615 zscale = 2.**sec4(7)
617 iskip = N1*8
619 nval = NDAT - (J+1)*(J+2)
621 call gbytes_g1(grec, ix, iskip, sec4(8), 0, nval)
622 ! sec4(8) is the number of bits used per datum value.
623 ! If sec4(8) = 255, assume they mean sec4(8) = 0
624 if (sec4(8) == 255) sec4(8) = 0
626 pdata(1:nval) = (float(ix(1:nval))*zscale)+xec4(1)
628 IUC = NDAT+1
629 IC = NVAL+1
630 DO JM=INFOGRID(1),0,-1
631 IL2=MAX(JM,J+1)
632 INUM=2*(INFOGRID(1)-IL2+1)
633 pdata(iuc-inum:iuc-1) = pdata(ic-inum:ic-1)
634 iuc = iuc - inum
635 ic = ic - inum
636 IUC = IUC-MAX((IL2-JM)*2,0)
637 ENDDO
639 if (ied.eq.0) then
640 iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
641 elseif (ied.eq.1) then
642 iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 18*8
643 endif
645 allocate(iexp(802))
646 allocate(imant(802))
647 ilast=j+1
648 do jrow=1,ilast
649 inext = 2*(ilast-jrow+1)
650 ! extract all the exponents
651 call gbytes_g1(grec, iexp, iskip, 8, 24, inext)
652 ! extract all the mantissas
653 ioff = 8
654 call gbytes_g1(grec, imant, iskip+8, 24, 8, inext)
655 iskip = iskip + inext*32
657 ! Build the real values from mantissas and exponents
658 bval = 2.**(-24)
659 i = 0
660 do jcol = jrow, infogrid(1)+1
661 index = index + 2
662 if (ilast.ge.jcol) then
663 i = i + 1
664 if ((iexp(i).eq.128.or.iexp(i).eq.0).and.(imant(i).eq.0)) then
665 pdata(i) = 0
666 else
667 if (iexp(i).ge.128) then
668 iexp(i) = iexp(i) - 128
669 isign = -1
670 else
671 isign = 1
672 endif
673 pdata(index) = isign*bval*IMANT(i)*16.**(IEXP(i)-64)
674 i = i + 1
675 if (iexp(i).ge.128) then
676 iexp(i) = iexp(i) - 128
677 isign = -1
678 else
679 isign = 1
680 endif
681 pdata(index+1) = isign*bval*IMANT(i)*16.**(IEXP(i)-64)
682 endif
683 endif
684 enddo
685 enddo
687 !Apply power scaling:
689 if (ipower.ne.0) then
690 power = float(ipower) / 1000.0
691 factor(0) = 1.0
692 do n = 1 , infogrid(1)
693 if( ipower .ne. 1000 ) then
694 factor(n) = 1.0 / (n * (n+1) )**power
695 else
696 factor(n) = 1.0 / (n * (n + 1))
697 endif
698 enddo
699 INDEX = -1
700 DO M = 0 , J-1
701 DO N = M , INFOGRID(1)
702 INDEX = INDEX + 2
703 IF ( N .GE. J ) THEN
704 PDATA(INDEX:INDEX+1) = PDATA(INDEX:INDEX+1) * FACTOR(N)
705 ENDIF
706 ENDDO
707 ENDDO
708 DO M = J , INFOGRID(1)
709 DO N = M , INFOGRID(1)
710 INDEX = INDEX + 2
711 PDATA(INDEX:INDEX+1) = PDATA(INDEX:INDEX+1) * FACTOR(N)
712 ENDDO
713 ENDDO
714 endif
716 end subroutine CSHUP
717 !
718 !=============================================================================!
719 !=============================================================================!
720 !=============================================================================!
721 !
722 !
723 !
724 ! Trigonometric functions which deal with degrees, rather than radians:
725 !
726 real function sind(theta)
727 real :: theta
728 sind = sin(theta*degran)
729 end function sind
730 real function cosd(theta)
731 real :: theta
732 cosd = cos(theta*degran)
733 end function cosd
734 real function tand(theta)
735 real :: theta
736 tand = tan(theta*degran)
737 end function tand
738 real function atand(x)
739 real :: x
740 atand = atan(x)*raddeg
741 end function atand
742 real function atan2d(x,y)
743 real :: x,y
744 atan2d = atan2(x,y)*raddeg
745 end function atan2d
746 real function asind(x)
747 real :: x
748 asind = asin(x)*raddeg
749 end function asind
750 real function acosd(x)
751 real :: x
752 acosd = acos(x)*raddeg
753 end function acosd
755 !
756 !=============================================================================!
757 !=============================================================================!
758 !=============================================================================!
759 !
760 end module module_grib
761 !
762 !=============================================================================!
763 !=============================================================================!
764 !=============================================================================!
765 !
766 subroutine gribget(nunit, ierr)
767 use module_grib
768 !-----------------------------------------------------------------------------!
769 ! !
770 ! Read a single GRIB record, with no unpacking of any header or data fields. !
771 ! !
772 ! Input: !
773 ! NUNIT: C unit number to read from. This should already be open. !
774 ! !
775 ! Output: !
776 ! IERR: Error flag, Non-zero means there was a problem with the read. !
777 ! !
778 ! Side Effects: !
779 ! The array GREC is allocated, and filled with one GRIB record. !
780 ! The C unit pointer is moved to the end of the GRIB record just read. !
781 ! !
782 ! Modules: !
783 ! MODULE_GRIB !
784 ! !
785 ! Externals: !
786 ! FINDGRIB !
787 ! BN_READ !
788 ! !
789 !-----------------------------------------------------------------------------!
791 implicit none
793 integer :: nunit
794 integer :: ierr
795 integer :: isz, isize
797 ! Position the file pointer at the beginning of the GRIB record.
798 call findgrib(nunit, isize, ierr)
799 if (ierr.ne.0) return
801 ! Allocate the GREC array to be able to hold the data
803 #if defined (CRAY)
804 #else
805 allocate(grec((isize+(mwsize/8-1))/(mwsize/8)))
806 #endif
808 ! Read the full GRIB record.
810 call bn_read(nunit, grec, isize, isz, ierr, 1)
812 !#if defined (DEC) || defined (ALPHA) || defined (alpha) || defined (LINUX)
813 #ifdef BYTESWAP
814 call swap4(grec, isz)
815 #endif
818 end subroutine gribget
819 !
820 !=============================================================================!
821 !=============================================================================!
822 !=============================================================================!
823 !
824 subroutine gribread(nunit, data, ndata, debug_level, ierr)
825 !-----------------------------------------------------------------------------!
826 ! Read one grib record, unpack the header and data information. !
827 ! !
828 ! Input: !
829 ! NUNIT: C Unit to read from. !
830 ! NDATA: Size of array DATA (Should be >= NDAT as computed herein.) !
831 ! !
832 ! Output: !
833 ! DATA: The unpacked data array !
834 ! IERR: Error flag, non-zero means there was a problem. !
835 ! !
836 ! Side Effects: !
837 ! * Header arrays SEC0, SEC1, SEC2, SEC3, SEC4, XEC4, INFOGRID and !
838 ! INFOGRID are filled. !
839 ! * The BITMAP array is filled. !
840 ! * The C unit pointer is advanced to the end of the GRIB record. !
841 ! !
842 ! Modules: !
843 ! MODULE_GRIB !
844 ! !
845 ! Externals: !
846 ! GRIBGET !
847 ! GRIBHEADER !
848 ! GRIBDATA !
849 ! !
850 !-----------------------------------------------------------------------------!
851 use module_grib
853 implicit none
855 integer :: nunit
856 integer :: debug_level
857 integer :: ierr
858 real, allocatable, dimension(:) :: datarray
859 integer :: ndata
860 real, dimension(ndata) :: data
862 integer :: ni, nj
864 ierr = 0
866 call gribget(nunit, ierr)
867 if (ierr.ne.0) return
869 ! Unpack the header information
871 call gribheader(debug_level,ierr)
873 ! Determine the size of the data array from the information in the header,
874 ! and allocate the array DATARRAY to hold that data.
876 if (sec2(4).ne.50) then
877 ni = infogrid(1)
878 nj = infogrid(2)
879 allocate(datarray(ni*nj))
880 else
881 ni = (infogrid(1)+1) * (infogrid(1)+2)
882 nj = 1
883 allocate(datarray(ni*nj))
884 endif
886 ! Unpack the data from the GRIB record, and fill the array DATARRAY.
888 call gribdata(datarray, ni*nj)
890 data(1:ni*nj) = datarray(1:ni*nj)
891 #if defined (CRAY)
892 #else
893 deallocate(grec, datarray)
894 #endif
896 end subroutine gribread
897 !
898 !=============================================================================!
899 !=============================================================================!
900 !=============================================================================!
901 !
902 subroutine get_sec1(ksec1)
903 ! Return the GRIB Section 1 header information, which has already been
904 ! unpacked by subroutine GRIBHEADER.
905 use module_grib
906 integer, dimension(100) :: ksec1
907 ksec1 = sec1
908 end subroutine get_sec1
909 !
910 !=============================================================================!
911 !=============================================================================!
912 !=============================================================================!
913 !
914 subroutine get_sec2(ksec2)
915 ! Return the GRIB Section 2 header information, which has already been
916 ! unpacked by subroutine GRIBHEADER.
917 use module_grib
918 integer, dimension(10) :: ksec2
919 ksec2 = sec2
920 end subroutine get_sec2
921 !
922 !=============================================================================!
923 !=============================================================================!
924 !=============================================================================!
925 !
926 subroutine get_gridinfo(iginfo, ginfo)
927 use module_grib
928 integer, dimension(40) :: iginfo
929 real, dimension(40) :: ginfo
930 iginfo = infogrid
931 ginfo = gridinfo
932 end subroutine get_gridinfo
933 !
934 !=============================================================================!
935 !=============================================================================!
936 !=============================================================================!
937 !
938 subroutine gribprint(isec)
939 use module_grib
940 implicit none
941 integer :: isec
942 integer :: ou = 6
943 character(len=12) :: string = ',t45,":",i8)'
944 character(len=15) :: rstring = ',t45,":",f12.5)'
946 if (isec.eq.0) then
947 write(*,'(/,"GRIB SECTION 0:")')
948 write(ou,'(5x,"Grib Length"'//string) sec0(1)
949 write(ou,'(5x,"Grib Edition"'//string) sec0(2)
950 else if (isec.eq.1) then
951 write(*,'(/,"GRIB SECTION 1:")')
952 write(ou,'(5x,"Length of PDS"'//string) sec1(1)
953 write(ou,'(5x,"Parameter Table Version"'//string) sec1(2)
954 write(ou,'(5x,"Center ID"'//string) sec1(3)
955 write(ou,'(5x,"Process ID"'//string) sec1(4)
956 write(ou,'(5x,"Grid ID"'//string) sec1(5)
957 if (sec1(25) == 1) then
958 write(ou,'(5x,"Is there a Grid Desc. Section (GDS)?",t45,": Yes")')
959 else if (sec1(25) == 0) then
960 write(ou,'(5x,"Is there a Grid Desc. Section (GDS)?",t45,": No")')
961 else
962 print*, 'Unrecognized sec1(25): ', sec1(25)
963 endif
964 if (sec1(26) == 1) then
965 write(ou,'(5x,"Is there a Bit Map Section (BMS)?",t45,": Yes")')
966 else if (sec1(26) == 0) then
967 write(ou,'(5x,"Is there a Bit Map Section (BMS)?",t45,": No")')
968 else
969 print*, 'Unrecognized sec1(26): ', sec1(26)
970 endif
971 write(ou,'(5x,"Parameter"'//string) sec1(7)
972 write(ou,'(5x,"Level type"'//string) sec1(8)
973 if ( (sec1(8) == 101) .or. (sec1(8) == 104) .or. (sec1(8) == 106) .or. &
974 (sec1(8) == 108) .or. (sec1(8) == 110) .or. (sec1(8) == 112) .or. &
975 (sec1(8) == 114) .or. (sec1(8) == 116) .or. (sec1(8) == 120) .or. &
976 (sec1(8) == 121) .or. (sec1(8) == 128) .or. (sec1(8) == 141) ) then
977 write(ou,'(5x,"Hgt, pres, etc. of layer top "'//string) sec1(9)
978 write(ou,'(5x,"Hgt, pres, etc. of layer bottom "'//string) sec1(10)
979 else
980 write(ou,'(5x,"Height, pressure, etc "'//string) sec1(9)
981 endif
982 write(ou,'(5x,"Year"'//string) sec1(11)
983 write(ou,'(5x,"Month"'//string) sec1(12)
984 write(ou,'(5x,"Day"'//string) sec1(13)
985 write(ou,'(5x,"Hour"'//string) sec1(14)
986 write(ou,'(5x,"Minute"'//string) sec1(15)
987 write(ou,'(5x,"Forecast time unit"'//string) sec1(16)
988 write(ou,'(5x,"P1"'//string) sec1(17)
989 write(ou,'(5x,"P2"'//string) sec1(18)
990 write(ou,'(5x,"Time Range Indicator"'//string) sec1(19)
991 write(ou,'(5x,"Number in Ave?"'//string) sec1(20)
992 write(ou,'(5x,"Number missing from ave?"'//string) sec1(21)
993 write(ou,'(5x,"Century"'//string) sec1(22)
994 write(ou,'(5x,"Sub-center"'//string) sec1(23)
995 write(ou,'(5x,"Decimal scale factor"'//string) sec1(24)
996 elseif ((isec.eq.2) .and. ((sec1(6).eq.128).or.(sec1(6).eq.192))) then
997 write(*,'(/,"GRIB SECTION 2:")')
998 write(ou,'(5x,"Length of GRID Desc. Section"'//string) sec2(1)
999 if ((sec2(2) /= 0).or.(sec2(3) /= 0) .or. (sec2(4) /= 0)) then
1000 write(ou,'(5x,"Number of V. Coordinate Parms"'//string) sec2(2)
1001 write(ou,'(5x,"List Starting point"'//string) sec2(3)
1002 write(ou,'(5x,"Data Representation type"'//string) sec2(4)
1003 endif
1005 if (sec2(4).eq.0) then
1006 write(ou,'(5x,"Cylindrical Equidistant Grid")')
1007 write(ou,'(10x,"NI"'//string) infogrid(1)
1008 write(ou,'(10x,"NJ"'//string) infogrid(2)
1009 write(ou,'(10x,"Lat 1"'//rstring) gridinfo(3)
1010 write(ou,'(10x,"Lon 1"'//rstring) gridinfo(4)
1011 write(ou,'(10x,"Resolution and Component:", t45,":",B8.8)') infogrid(5)
1012 write(ou,'(10x,"Lat NI"'//string) infogrid(6)
1013 write(ou,'(10x,"Lon NJ"'//string) infogrid(7)
1014 write(ou,'(10x,"Delta-Lon"'//string) infogrid(8)
1015 write(ou,'(10x,"Delta-Lat"'//string) infogrid(9)
1016 write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
1017 write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
1018 write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
1020 else if (sec2(4).eq.1) then
1021 write(ou,'(5x,"Mercator Grid")')
1022 write(ou,'(10x,"NI"'//string) infogrid(1)
1023 write(ou,'(10x,"NJ"'//string) infogrid(2)
1024 write(ou,'(10x,"Lat 1"'//rstring) gridinfo(3)
1025 write(ou,'(10x,"Lon 1"'//rstring) gridinfo(4)
1026 write(ou,'(10x,"Resolution and Component",t45,":", B8.8)') infogrid(5)
1027 write(ou,'(10x,"Lat NI"'//rstring) gridinfo(6)
1028 write(ou,'(10x,"Lon NJ"'//rstring) gridinfo(7)
1029 write(ou,'(10x,"Dx"'//rstring) gridinfo(8)
1030 write(ou,'(10x,"Dy"'//rstring) gridinfo(9)
1031 write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
1032 write(ou,'(10x,"Latin"'//rstring) gridinfo(11)
1033 write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
1034 write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
1036 else if (sec2(4).eq.4) then
1037 write(ou,'(5x,"Gaussian Grid")')
1038 write(ou,'(10x,"NI"'//string) infogrid(1)
1039 write(ou,'(10x,"NJ"'//string) infogrid(2)
1040 write(ou,'(10x,"Original (stored) Lat 1"'//rstring) gridinfo(18)
1041 write(ou,'(10x,"Lat 1"'//rstring) gridinfo(3)
1042 write(ou,'(10x,"Lon 1"'//rstring) gridinfo(4)
1043 write(ou,'(10x,"Resolution and Component",t45,":", B8.8)') infogrid(5)
1044 write(ou,'(10x,"Original (stored) Lat NI"'//rstring) gridinfo(17)
1045 write(ou,'(10x,"Lat NI"'//rstring) gridinfo(6)
1046 write(ou,'(10x,"Lon NJ"'//rstring) gridinfo(7)
1047 write(ou,'(10x,"Delta-Lon"'//rstring) gridinfo(8)
1048 write(ou,'(10x,"Delta-Lat"'//rstring) gridinfo(19)
1049 write(ou,'(10x,"Number of lats (pole - eq)"'//string) infogrid(9)
1050 write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
1051 write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
1052 write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
1053 elseif (sec2(4).eq.3) then
1054 write(ou,'(5x,"Lambert Conformal Grid")')
1055 write(ou,'(10x,"NI"'//string) infogrid(1)
1056 write(ou,'(10x,"NJ"'//string) infogrid(2)
1057 write(ou,'(10x,"Lat 1"'//string) infogrid(3)
1058 write(ou,'(10x,"Lon 1"'//string) infogrid(4)
1059 write(ou,'(10x,"Resolution and Component",t45,":", B8.8)') infogrid(5)
1060 write(ou,'(10x,"Lov"'//string) infogrid(6)
1061 write(ou,'(10x,"Dx"'//string) infogrid(7)
1062 write(ou,'(10x,"Dy"'//string) infogrid(8)
1063 write(ou,'(10x,"Projection center"'//string) infogrid(9)
1064 write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
1065 write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
1066 write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
1067 write(ou,'(10x,"Latin 1"'//string) infogrid(11)
1068 write(ou,'(10x,"Latin 2"'//string) infogrid(12)
1069 write(ou,'(10x,"Lat of southern pole"'//string) infogrid(13)
1070 write(ou,'(10x,"Lon of southern pole"'//string) infogrid(14)
1071 elseif (sec2(4).eq.5) then
1072 write(ou,'(5x,"Polar Stereographic Grid")')
1073 write(ou,'(10x,"NI"'//string) infogrid(1)
1074 write(ou,'(10x,"NJ"'//string) infogrid(2)
1075 write(ou,'(10x,"Lat 1"'//string) infogrid(3)
1076 write(ou,'(10x,"Lon 1"'//string) infogrid(4)
1077 write(ou,'(10x,"Resolution and Component", t45,":",B8.8)') infogrid(5)
1078 write(ou,'(10x,"Lov"'//string) infogrid(6)
1079 write(ou,'(10x,"Dx"'//string) infogrid(7)
1080 write(ou,'(10x,"Dy"'//string) infogrid(8)
1081 write(ou,'(10x,"Projection center"'//string) infogrid(9)
1082 write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
1083 write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
1084 write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
1085 elseif (sec2(4).eq.50) then
1086 write(ou,'(5x,"Spherical harmonic components")')
1087 write(ou,'(10x,"J-Pentagonal resolution parm:"'//string) infogrid(1)
1088 write(ou,'(10x,"K-Pentagonal resolution parm:"'//string) infogrid(2)
1089 write(ou,'(10x,"M-Pentagonal resolution parm:"'//string) infogrid(3)
1090 write(ou,'(10x,"Representation type"'//string) infogrid(4)
1091 write(ou,'(10x,"Coefficient storage mode"'//string) infogrid(5)
1092 endif
1093 elseif ((isec.eq.3) .and. (sec1(26).eq.1)) then
1094 write(*,'(/,"GRIB SECTION 3:")')
1095 write(ou,'(5x,"Length of bit map section"'//string) sec3(1)
1096 write(ou,'(5x,"Number of unused bits"'//string) sec3(2)
1097 write(ou,'(5x,"Numeric"'//string) sec3(3)
1099 elseif (isec.eq.4) then
1100 write(*,'(/,"GRIB SECTION 4:")')
1101 write(ou,'(5x,"Length of BDS"'//string) sec4(1)
1102 write(ou,'(5x,"0/1: grid-point or sph. harm. data"'//string) sec4(2)
1103 write(ou,'(5x,"0/1: simple or complex packing"'//string) sec4(3)
1104 write(ou,'(5x,"0/1: floating or integer"'//string) sec4(4)
1105 write(ou,'(5x,"0/1: No addl flags or addl flags"'//string) sec4(5)
1106 write(ou,'(5x,"Unused bits"'//string) sec4(6)
1107 write(ou,'(5x,"Binary Scale Factor"'//string) sec4(7)
1108 write(ou,'(5x,"Reference Value", t45, ":", F18.8)') xec4(1)
1109 write(ou,'(5x,"Number of bits for packing"'//string) sec4(8)
1110 endif
1112 end subroutine gribprint
1113 !
1114 !=============================================================================!
1115 !=============================================================================!
1116 !=============================================================================!
1117 !
1118 subroutine get_bitmap(bm8, ndat)
1119 use module_grib
1120 integer, dimension(ndat) :: bm8
1121 if ((sec1(6).eq.64).or.(sec1(6).eq.192)) then
1122 bm8 = bitmap
1123 else
1124 bm8 = 1
1125 endif
1126 end subroutine get_bitmap
1127 !
1128 !=============================================================================!
1129 !=============================================================================!
1130 !=============================================================================!
1131 !
1132 subroutine gribxyll(x, y, xlat, xlon)
1133 use module_grib
1134 implicit none
1136 real , intent(in) :: x, y
1137 real , intent(out) :: xlat, xlon
1139 real :: r, xkm, ykm, y1
1140 integer :: iscan, jscan
1142 if (sec2(4).eq.0) then ! Cylindrical equidistant grid
1144 xlat = gridinfo(3) + gridinfo(9)*(y-1.)
1145 xlon = gridinfo(4) + gridinfo(8)*(x-1.)
1147 elseif (sec2(4) == 1) then ! Mercator grid
1148 r = grrth*cosd(gtrue1)
1149 xkm = (x-1.)*gridinfo(8)
1150 ykm = (y-1.)*gridinfo(9)
1151 xlon = gridinfo(4) + (xkm/r)*(180./pi)
1152 y1 = r*alog((1.+sind(gridinfo(3)))/cosd(gridinfo(3)))/gridinfo(9)
1153 xlat = 90. - 2. * atan(exp(-gridinfo(9)*(y+y1-1.)/r))*180./pi
1155 elseif (sec2(4) == 3) then ! Lambert Conformal grid
1156 gclon = gridinfo(6)
1157 r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
1158 xlat = 90. - 2.*atand(tand(45.-gtrue1/2.)* &
1159 ((r*gkappa*gridinfo(7))/(grrth*sind(90.-gtrue1)))**(1./gkappa))
1160 xlon = atan2d((x-1.+gx1),-(y-1.+gy1))/gkappa + gclon
1162 elseif (sec2(4) == 5) then ! Polar Stereographic grid
1163 gclon = gridinfo(6)
1164 r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
1165 xlat = 90. - 2.*atan2d((r*gridinfo(7)),(grrth*(1.+sind(gtrue1))))
1166 xlon = atan2d((x-1.+gx1),-(y-1.+gy1)) + gclon
1168 elseif (sec2(4) == 4) then ! Gaussian grid
1170 xlon = gridinfo(4) + gridinfo(8)*(x-1.)
1171 xlat = gridinfo(3) + gridinfo(19)*(y-1.)
1173 else
1174 write(*,'("Unrecognized projection:", I10)') sec2(4)
1175 write(*,'("STOP in GRIBXYLL")')
1176 stop
1177 endif
1179 end subroutine gribxyll
1180 !
1181 !=============================================================================!
1182 !=============================================================================!
1183 !=============================================================================!
1184 !
1185 subroutine gribllxy(xlat, xlon, x, y)
1186 use module_grib
1187 implicit none
1188 real , intent(in) :: xlat, xlon
1189 real , intent(out) :: x, y
1191 real :: r, y1
1193 if (sec2(4) == 0) then ! Cylindrical Equidistant grid
1195 x = 1. + (xlon-gridinfo(4)) / gridinfo(9)
1196 y = 1. + (xlat-gridinfo(3)) / gridinfo(8)
1198 else if (sec2(4) == 1) then ! Mercator grid
1200 r = grrth*cosd(gtrue1)
1201 x = 1.+( (r/gridinfo(8)) * (xlon-gridinfo(4)) * (pi/180.) )
1202 y1 = (r/gridinfo(9))*alog((1.+sind(gridinfo(3)))/cosd(gridinfo(3)))
1203 y = 1. + ((r/gridinfo(9))*alog((1.+sind(xlat))/cosd(xlat)))-y1
1205 else if (sec2(4) == 3) then ! Lambert Conformal grid
1206 gclon = gridinfo(6)
1207 r = grrth/(gridinfo(7)*gkappa)*sind(90.-gtrue1) * &
1208 (tand(45.-xlat/2.)/tand(45.-gtrue1/2.)) ** gkappa
1209 x = r*sind(gkappa*(xlon-gclon)) - gx1 + 1.
1210 y = -r*cosd(gkappa*(xlon-gclon)) - gy1 + 1.
1212 elseif (sec2(4) == 5) then ! Polar Stereographic grid
1213 gclon = gridinfo(6)
1214 r = grrth/gridinfo(7) * tand((90.-xlat)/2.) * (1.+sind(gtrue1))
1215 x = ( r * sind(xlon-gclon)) - gx1 + 1.
1216 y = (-r * cosd(xlon-gclon)) - gy1 + 1.
1218 else
1219 write(*,'("Unrecognized projection:", I10)') sec2(4)
1220 write(*,'("STOP in GRIBLLXY")')
1221 stop
1222 endif
1224 end subroutine gribllxy
1225 !
1226 !=============================================================================!
1227 !=============================================================================!
1228 !=============================================================================!
1229 !
1230 subroutine glccone (fsplat,ssplat,sign,confac)
1231 use module_grib
1232 implicit none
1233 real, intent(in) :: fsplat,ssplat
1234 integer, intent(in) :: sign
1235 real, intent(out) :: confac
1236 if (abs(fsplat-ssplat).lt.1.E-3) then
1237 confac = sind(fsplat)
1238 else
1239 confac = log10(cosd(fsplat))-log10(cosd(ssplat))
1240 confac = confac/(log10(tand(45.-float(sign)*fsplat/2.))- &
1241 log10(tand(45.-float(sign)*ssplat/2.)))
1242 endif
1243 end subroutine glccone
1244 !
1245 !=============================================================================!
1246 !=============================================================================!
1247 !=============================================================================!
1248 !
1249 !=============================================================================!
1250 !=============================================================================!
1251 !=============================================================================!
1252 !
1253 subroutine gribheader(debug_level,ierr)
1254 !
1255 ! IERR non-zero means there was a problem unpacking the grib header.
1256 !
1257 use module_grib
1258 implicit none
1259 integer :: debug_level
1260 integer :: ierr
1262 integer, parameter :: nsec1 = 24
1264 integer, dimension(nsec1) :: &
1265 iw1=(/3,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,2/)
1266 integer :: icount, iskip, ibts, nbm, nbz, i9skip, i17skip
1268 integer :: iman, ichar, isign, iscan
1270 integer, allocatable, dimension(:) :: bm8
1272 real :: r
1273 integer :: isvns
1274 integer :: gsvns = 0
1276 if (gsvns.eq.0) then
1277 if (mwsize.eq.32) then
1278 gsvns = transfer('7777', gsvns)
1279 elseif(mwsize.eq.64) then
1280 call gbyte_g1(char(0)//char(0)//char(0)//char(0)//'7777', gsvns, 0, mwsize)
1281 endif
1282 endif
1284 ! Section 0:
1285 sec0(2) = ied
1286 if (ied.eq.1) then
1287 call gbyte_g1(grec, sec0(1), 32, 24)
1288 iskip = 64
1289 elseif (ied.eq.0) then
1290 sec0(1) = gribsize(grec,200, ierr)
1291 iskip = 32
1292 endif
1294 ! Section 1:
1295 i9skip = iskip + 80
1296 i17skip = iskip + 144
1297 do icount = 1, nsec1 - ((1-ied)*6)
1298 ibts = iw1(icount)*8
1299 call gbyte_g1(grec, sec1(icount), iskip, ibts)
1300 iskip = iskip + ibts
1301 enddo
1302 if (ied.eq.0) sec1(22) = 20
1303 ! Sec1 indices 9 and 10 might actually be one value, not two.
1304 ! If this is the case, reread sec1(9), and set sec1(10) to zero:
1305 if ( (sec1(8) == 101) .or. (sec1(8) == 104) .or. (sec1(8) == 106) .or. &
1306 (sec1(8) == 108) .or. (sec1(8) == 110) .or. (sec1(8) == 112) .or. &
1307 (sec1(8) == 114) .or. (sec1(8) == 116) .or. (sec1(8) == 120) .or. &
1308 (sec1(8) == 121) .or. (sec1(8) == 128) .or. (sec1(8) == 141) .or. &
1309 (sec1(8) == 236) ) then
1310 ! No action here.
1311 else
1312 call gbyte_g1(grec, sec1(9), i9skip, 16)
1313 sec1(10) = 0.
1314 endif
1316 if (sec1(24).ge.32768) sec1(24) = 32768-sec1(24)
1318 ! If the TIME/RANGE INDICATOR (sec1(19)) indicates that the time P1
1319 ! is spread over two bytes, then recompute sec1(17) and set sec1(18)
1320 ! to zero.
1321 if (sec1(19) == 10) then
1322 call gbyte_g1(grec, sec1(17), i17skip, 16)
1323 sec1(18) = 0
1324 endif
1326 ! Pull out single bits from sec1(6) for the GDS and BMS flags:
1327 sec1(25) = sec1(6)/128
1328 sec1(26) = mod(sec1(6)/64,2)
1330 ! Section 2:
1331 ! if ((sec1(6) == 128) .or. (sec1(6) == 192)) then
1332 if (sec1(25) == 1) then
1334 if (ied.eq.0) then
1335 iskip = 32 + sec1(1)*8
1336 elseif (ied.eq.1) then
1337 iskip = 64 + sec1(1)*8
1338 endif
1339 call gbyte_g1(grec, sec2(1), iskip, 24)
1340 iskip = iskip + 24
1341 call gbytes_g1(grec, sec2(2), iskip, 8, 0, 3)
1342 iskip = iskip + 8*3
1344 if (sec2(4) == 0) then
1345 ! Lat/Lon Grid:
1346 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)
1347 iskip = iskip + 32
1348 call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
1349 iskip = iskip + 48
1350 call gbyte_g1(grec, infogrid(5), iskip, 8)
1351 iskip = iskip + 8
1352 call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 2)
1353 iskip = iskip + 48
1354 call gbytes_g1(grec, infogrid(8), iskip, 16, 0, 2)
1355 iskip = iskip + 32
1356 call gbyte_g1(grec, infogrid(21), iskip, 1)
1357 infogrid(21) = 1-(infogrid(21)*2)
1358 iskip = iskip + 1
1359 call gbyte_g1(grec, infogrid(22), iskip, 1)
1360 infogrid(22) = (infogrid(22)*2)-1
1361 iskip = iskip + 1
1362 call gbyte_g1(grec, infogrid(10), iskip, 1)
1363 iskip = iskip + 1
1364 iskip = iskip + 5
1365 call gbyte_g1(grec, infogrid(11), iskip, 32)
1366 iskip = iskip + 32
1368 !MGD if ( debug_level .gt. 100 ) then
1369 !MGD print *, "lat/lon grib grid info", infogrid(1), infogrid(3), &
1370 !MGD infogrid(5), infogrid(6), infogrid(8), infogrid(21), &
1371 !MGD infogrid(22), infogrid(10), infogrid(11), infogrid(8)
1372 !MGD end if
1374 infogrid(8) = infogrid(8) * infogrid(21)
1375 infogrid(9) = infogrid(9) * infogrid(22)
1377 gridinfo(1) = float(infogrid(1))
1378 gridinfo(2) = float(infogrid(2))
1379 if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
1380 if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
1381 gridinfo(3) = float(infogrid(3))*0.001
1382 gridinfo(4) = infogrid(4) * 0.001
1383 if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)
1384 if (infogrid(7).ge.8388608) infogrid(7) = 8388608 - infogrid(7)
1385 gridinfo(6) = infogrid(6) * 0.001
1386 gridinfo(7) = infogrid(7) * 0.001
1387 gridinfo(8) = infogrid(8) * 0.001
1388 gridinfo(9) = infogrid(9) * 0.001
1389 gridinfo(21) = float(infogrid(21))
1390 gridinfo(22) = float(infogrid(22))
1391 elseif (sec2(4) == 1) then ! Mercator grid
1392 ! Number of points in X and Y
1393 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)
1394 iskip = iskip + 32
1395 ! Starting lat and lon
1396 call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
1397 iskip = iskip + 48
1398 ! "Resolution and component flags"
1399 call gbyte_g1(grec, infogrid(5), iskip, 8)
1400 iskip = iskip + 8
1401 ! Ending lat and lon
1402 call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 2)
1403 iskip = iskip + 48
1404 ! Truelat, 3 bytes
1405 call gbyte_g1(grec, infogrid(11), iskip, 24)
1406 iskip = iskip + 24
1407 ! "Reserved", i.e., skip a byte
1408 iskip = iskip + 8
1409 ! Scanning mode flags, first three bits of the next byte
1410 ! and skip the last five bits.
1411 call gbyte_g1(grec, infogrid(21), iskip, 1)
1412 infogrid(21) = 1-(infogrid(21)*2)
1413 iskip = iskip + 1
1414 call gbyte_g1(grec, infogrid(22), iskip, 1)
1415 infogrid(22) = (infogrid(22)*2)-1
1416 iskip = iskip + 1
1417 call gbyte_g1(grec, infogrid(10), iskip, 1)
1418 iskip = iskip + 1
1419 iskip = iskip + 5
1420 ! Grid increment in X and Y
1421 call gbytes_g1(grec, infogrid(8), iskip, 24, 0, 2)
1422 iskip = iskip + 48
1423 ! Done reading map specifications.
1424 ! Now do various conversions:
1426 gridinfo(1) = float(infogrid(1)) ! ok
1427 gridinfo(2) = float(infogrid(2)) ! ok
1429 if (infogrid(3) .ge.8388608) infogrid(3) = 8388608 - infogrid(3)
1430 if (infogrid(4) .ge.8388608) infogrid(4) = 8388608 - infogrid(4)
1431 if (infogrid(6) .ge.8388608) infogrid(6) = 8388608 - infogrid(6)
1432 if (infogrid(7) .ge.8388608) infogrid(7) = 8388608 - infogrid(7)
1433 if (infogrid(11).ge.8388608) infogrid(11) = 8388608 - infogrid(11)
1434 gridinfo(3) = infogrid(3) * 0.001
1435 gridinfo(4) = infogrid(4) * 0.001
1436 gridinfo(6) = infogrid(6) * 0.001
1437 gridinfo(7) = infogrid(7) * 0.001
1438 gridinfo(8) = infogrid(8) * 0.001
1439 gridinfo(9) = infogrid(9) * 0.001
1440 gridinfo(11) = infogrid(11) * 0.001
1442 gridinfo(21) = infogrid(21)
1443 gridinfo(22) = infogrid(22)
1445 gridinfo(20) = 6370.949
1446 grrth = gridinfo(20)
1447 gtrue1 = gridinfo(11)
1449 elseif (sec2(4) == 3) then
1450 if (ied.eq.0) then
1451 print '(//,"*** Despair ***"//)'
1452 stop
1453 endif
1454 ! Lambert Conformal:
1455 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)
1456 iskip = iskip + 32
1457 call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
1458 iskip = iskip + 48
1459 if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
1460 if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
1461 call gbyte_g1(grec, infogrid(5), iskip, 8)
1462 iskip = iskip + 8
1463 call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 3)
1464 if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)
1465 iskip = iskip + 72
1466 call gbyte_g1(grec, infogrid(9), iskip, 8)
1467 iskip = iskip + 8
1468 call gbyte_g1(grec, infogrid(21), iskip, 1)
1469 infogrid(21) = 1-(infogrid(21)*2)
1470 iskip = iskip + 1
1471 call gbyte_g1(grec, infogrid(22), iskip, 1)
1472 infogrid(22) = (infogrid(22)*2)-1
1473 iskip = iskip + 1
1474 call gbyte_g1(grec, infogrid(10), iskip, 1)
1475 iskip = iskip + 1
1476 iskip = iskip + 5
1477 call gbytes_g1(grec, infogrid(11), iskip, 24, 0, 4)
1478 if (infogrid(11).ge.8388608) infogrid(11) = 8388608 - infogrid(11)
1479 if (infogrid(12).ge.8388608) infogrid(12) = 8388608 - infogrid(12)
1480 if (infogrid(13).ge.8388608) infogrid(13) = 8388608 - infogrid(13)
1481 if (infogrid(14).ge.8388608) infogrid(14) = 8388608 - infogrid(14)
1482 iskip = iskip + 96
1483 call gbyte_g1(grec, infogrid(15), iskip, 16)
1484 iskip = iskip + 16
1486 infogrid(7) = infogrid(7) * infogrid(21)
1487 infogrid(8) = infogrid(8) * infogrid(22)
1490 gridinfo(1) = float(infogrid(1))
1491 gridinfo(2) = float(infogrid(2))
1492 gridinfo(3) = infogrid(3) * 0.001
1493 gridinfo(4) = infogrid(4) * 0.001
1494 gridinfo(6) = infogrid(6) * 0.001
1495 gridinfo(7) = infogrid(7) * 0.001
1496 gridinfo(8) = infogrid(8) * 0.001
1497 gridinfo(9) = infogrid(9) * 0.001
1498 gridinfo(11) = infogrid(11) * 0.001
1499 gridinfo(12) = infogrid(12) * 0.001
1500 gridinfo(13) = infogrid(13) * 0.001
1501 gridinfo(14) = infogrid(14) * 0.001
1503 gridinfo(20) = 6370
1504 ! a priori knowledge:
1505 if (sec1(5).eq.212) then
1506 gridinfo(3) = 12.190
1507 gridinfo(4) = -133.459
1508 gridinfo(7) = 40.63525
1509 gridinfo(8) = 40.63525
1510 gridinfo(20) = 6370
1511 endif
1513 !=============================================================================!
1514 ! More a priori knowledge: !
1515 ! Correct some bad lat/lon numbers coded into some RUC headers. !
1516 ! !
1517 if (sec1(3) == 59) then ! If FSL
1518 if (sec1(4) == 86) then ! and RUC
1519 if (sec1(5) == 255) then
1520 ! Check to correct bad lat/lon numbers.
1521 if (infogrid(3) == 16909) then
1522 infogrid(3) = 16281
1523 gridinfo(3) = 16.281
1524 endif
1525 if (infogrid(4) == 236809) then
1526 infogrid(4) = 2338622
1527 gridinfo(4) = 233.8622
1528 endif
1529 endif
1530 endif
1531 endif
1532 !=============================================================================!
1535 gridinfo(21) = float(infogrid(21))
1536 gridinfo(22) = float(infogrid(22))
1538 ! Map parameters
1539 glat1 = gridinfo(3)
1540 glon1 = gridinfo(4)
1541 gclon = gridinfo(6)
1542 if (gclon.gt.180.) gclon = -(360.-gclon)
1543 if ((gclon<0).and.(glon1>180)) glon1 = glon1-360.
1544 gtrue1 = gridinfo(11)
1545 gtrue2 = gridinfo(12)
1546 grrth = gridinfo(20)
1547 call glccone(gtrue1, gtrue2, 1, gkappa)
1548 r = grrth/(gridinfo(7)*gkappa)*sind(90.-gtrue1) * &
1549 (tand(45.-glat1/2.)/tand(45.-gtrue1/2.)) ** gkappa
1550 gx1 = r*sind(gkappa*(glon1-gclon))
1551 gy1 = -r*cosd(gkappa*(glon1-gclon))
1553 elseif (sec2(4) == 4) then
1554 ! Gaussian Grid:
1555 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)
1556 iskip = iskip + 32
1557 call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
1558 iskip = iskip + 48
1559 call gbyte_g1(grec, infogrid(5), iskip, 8)
1560 iskip = iskip + 8
1561 call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 2)
1562 iskip = iskip + 48
1563 call gbytes_g1(grec, infogrid(8), iskip, 16, 0, 2)
1564 iskip = iskip + 32
1565 call gbyte_g1(grec, infogrid(21), iskip, 1)
1566 infogrid(21) = 1-(infogrid(21)*2)
1567 iskip = iskip + 1
1568 call gbyte_g1(grec, infogrid(22), iskip, 1)
1569 infogrid(22) = (infogrid(22)*2)-1
1570 iskip = iskip + 1
1571 call gbyte_g1(grec, infogrid(10), iskip, 1)
1572 iskip = iskip + 1
1573 iskip = iskip + 5
1574 call gbyte_g1(grec, infogrid(11), iskip, 32)
1575 iskip = iskip + 32
1577 infogrid(8) = infogrid(8) * infogrid(21)
1579 gridinfo(1) = float(infogrid(1))
1580 gridinfo(2) = float(infogrid(2))
1581 if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
1582 if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
1583 gridinfo(3) = float(infogrid(3))*0.001
1584 gridinfo(4) = infogrid(4) * 0.001
1585 if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)
1586 if (infogrid(7).ge.8388608) infogrid(7) = 8388608 - infogrid(7)
1587 gridinfo(6) = infogrid(6) * 0.001
1588 gridinfo(7) = infogrid(7) * 0.001
1589 gridinfo(8) = infogrid(8) * 0.001
1590 gridinfo(21) = float(infogrid(21))
1591 gridinfo(22) = float(infogrid(22))
1593 ! Compute an approximate delta-latitude and starting latitude.
1594 ! Replace the stored value of starting latitude with approximate one.
1595 gridinfo(18) = gridinfo(3)
1596 infogrid(18) = infogrid(3)
1597 gridinfo(17) = gridinfo(6)
1598 infogrid(17) = infogrid(6)
1599 ! call griblgg(infogrid(2), gridinfo(3), gridinfo(19))
1600 ! infogrid(19) = nint(gridinfo(19)*1000.)
1601 ! infogrid(3) = nint(gridinfo(3)*1000.)
1602 gridinfo(6) = -gridinfo(3)
1603 infogrid(6) = -infogrid(3)
1605 elseif (sec2(4) == 5) then
1606 ! Polar Stereographic Grid
1607 if (ied.eq.0) then
1608 print '(//,"*** Despair ***"//)'
1609 stop
1610 endif
1611 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2) ! NX and NY
1612 iskip = iskip + 32
1613 call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2) ! LAT1 and LON1
1614 iskip = iskip + 48
1615 call gbyte_g1(grec, infogrid(5), iskip, 8) ! Resolution and Component
1616 iskip = iskip + 8
1617 call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 3) ! LOV, DX, and DY
1618 iskip = iskip + 72
1619 call gbyte_g1(grec, infogrid(9), iskip, 8) ! Projection center flag
1620 iskip = iskip + 8
1621 call gbyte_g1(grec, infogrid(21), iskip, 1)
1622 infogrid(21) = 1-(infogrid(21)*2)
1623 iskip = iskip + 1
1624 call gbyte_g1(grec, infogrid(22), iskip, 1)
1625 infogrid(22) = (infogrid(22)*2)-1
1626 iskip = iskip + 1
1627 call gbyte_g1(grec, infogrid(10), iskip, 1)
1628 iskip = iskip + 1
1629 iskip = iskip + 5
1630 ! call gbyte_g1(grec, infogrid(11), iskip, 32) ! Set to 0 (reserved)
1631 iskip = iskip + 32
1633 if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
1634 if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
1635 if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)
1638 infogrid(7) = infogrid(7) * infogrid(21)
1639 infogrid(8) = infogrid(8) * infogrid(22)
1641 gridinfo(1) = float(infogrid(1))
1642 gridinfo(2) = float(infogrid(2))
1643 gridinfo(3) = infogrid(3) * 0.001
1644 gridinfo(4) = infogrid(4) * 0.001
1645 gridinfo(6) = infogrid(6) * 0.001
1646 gridinfo(7) = infogrid(7) * 0.001
1647 gridinfo(8) = infogrid(8) * 0.001
1649 gridinfo(20) = 6370
1651 ! a priori knowledge:
1652 if (sec1(5).eq.240) then
1653 gridinfo(3) = 22.7736
1654 gridinfo(4) = -120.376
1655 gridinfo(7) = 4.7625
1656 gridinfo(8) = 4.7625
1657 gridinfo(20) = 6370
1658 endif
1660 ! Map parameters
1661 glat1 = gridinfo(3)
1662 glon1 = gridinfo(4)
1663 gclon = gridinfo(6)
1664 if (gclon.gt.180.) gclon = -(360.-gclon)
1665 ! GRIB edition 1 Polar Stereographic grids are true at 60 degrees
1666 ! Which hemisphere depends on infogrid(9), the "Projection Center Flag"
1667 grrth = gridinfo(20)
1668 if (infogrid(9) > 127) then
1669 gtrue1 = -60.
1670 r = grrth/gridinfo(7) * tand((-90.-glat1)/2.) * (1.+sind(-gtrue1))
1671 gx1 = -r * sind(glon1-gridinfo(6))
1672 gy1 = -r * cosd(glon1-gridinfo(6))
1673 else
1674 gtrue1 = 60.
1675 r = grrth/gridinfo(7) * tand((90.-glat1)/2.) * (1.+sind(gtrue1))
1676 gx1 = r * sind(glon1-gridinfo(6))
1677 gy1 = -r * cosd(glon1-gridinfo(6))
1678 endif
1680 gridinfo(21) = float(infogrid(21))
1681 gridinfo(22) = float(infogrid(22))
1683 elseif (sec2(4) == 50) then
1684 ! Spherical harmonic coefficients
1685 if (ied.eq.0) then
1686 print '(//,"*** Despair ***"//)'
1687 stop
1688 endif
1689 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 3)
1690 iskip = iskip + 48
1691 call gbytes_g1(grec, infogrid(4), iskip, 8, 0, 2)
1692 iskip = iskip + 16
1694 iskip = iskip + 18*8
1696 else
1697 call gribprint(0)
1698 call gribprint(1)
1699 call gribprint(2)
1700 call gribprint(3)
1701 call gribprint(4)
1702 write(*,'("Unrecognized grid: ", i8)') sec2(4)
1703 write(*,'("This grid is not currently supported.")')
1704 write(*,'("Write your own program to put the data to the intermediate format")')
1705 stop
1706 endif
1708 endif
1710 ! Section 3
1711 if ((sec1(6).eq.64).or.(sec1(6).eq.192)) then
1712 if (ied.eq.0) then
1713 print '(//,"*** Despair ***"//)'
1714 stop
1715 endif
1717 if (ied.eq.0) then
1718 iskip = 32 + sec1(1)*8 + sec2(1)*8
1719 elseif (ied.eq.1) then
1720 iskip = 64 + sec1(1)*8 + sec2(1)*8
1721 endif
1722 call gbyte_g1(grec, sec3(1), iskip, 24)
1723 iskip = iskip + 24
1724 call gbyte_g1(grec, sec3(2), iskip, 8)
1725 iskip = iskip + 8
1726 call gbyte_g1(grec, sec3(3), iskip, 16)
1727 iskip = iskip + 16
1729 #if defined (CRAY)
1730 #else
1731 allocate(bitmap((sec3(1)-6)*8))
1732 #endif
1733 allocate(bm8((sec3(1)-6)*8))
1734 call gbytes_g1(grec, bm8, iskip, 1, 0, (sec3(1)-6)*8)
1735 bitmap(1:size(bm8)) = bm8(1:size(bm8))
1736 deallocate(bm8)
1737 iskip = iskip + sec3(1)-6
1738 else
1739 sec3 = 0
1740 endif
1742 ! Section 4
1743 if ((sec1(6).eq.128).or.(sec1(6).eq.192)) then
1744 if (ied.eq.0) then
1745 iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8
1746 elseif (ied.eq.1) then
1747 iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8
1748 endif
1749 call gbyte_g1(grec, sec4(1), iskip, 24)
1750 if (sec4(1) > (sec0(1) - sec1(1) - sec2(1) - sec3(1) - 4)) then
1751 write(*,'(/,"*** I have good reason to believe that this GRIB record is")')
1752 write(*,'("*** corrupted or miscoded.",/)')
1753 ierr = 1
1754 return
1755 endif
1756 iskip = iskip + 24
1757 call gbytes_g1(grec, sec4(2), iskip, 1,0,4)
1758 iskip = iskip + 4
1759 call gbyte_g1(grec, sec4(6), iskip, 4)
1760 iskip = iskip + 4
1761 ! Get the binary scale factor
1762 call gbyte_g1(grec, isign, iskip, 1)
1763 iskip = iskip + 1
1764 call gbyte_g1(grec, sec4(7), iskip, 15)
1765 iskip = iskip + 15
1766 sec4(7) = sec4(7) * (-2*isign+1)
1767 ! Get the reference value:
1768 call gbyte_g1(grec, isign, iskip, 1)
1769 iskip = iskip + 1
1770 isign = -2*isign+1
1771 call gbyte_g1(grec, ichar, iskip, 7)
1772 iskip = iskip + 7
1773 call gbyte_g1(grec, iman, iskip, 24)
1774 iskip = iskip + 24
1775 if ( iman .ne. 0 ) then
1776 xec4(1) = float(isign) * (2.**(-24)) * float(iman) * &
1777 (16.**(ichar-64))
1778 else
1779 xec4(1) = 0.
1780 endif
1782 call gbyte_g1(grec,sec4(8), iskip, 8)
1783 ! sec4(8) is the number of bits used per datum value.
1784 ! If sec4(8) = 255, assume they mean sec4(8) = 0
1785 if (sec4(8) == 255) sec4(8) = 0
1786 iskip = iskip + 8
1787 endif
1789 ! Section 5
1790 call gbyte_g1(grec, isvns, ((sec0(1)-4)*8), 32)
1791 if (isvns.ne.gsvns) then
1792 write(*, '("End-of-record mark (7777) not found", 2I10)') isvns
1793 write(*, '("Sec0(1) = ", I8, i2)') sec0(1), sevencount
1794 sevencount = sevencount + 1
1795 if (sevencount > 10) then
1796 write(*,'(//," *** Found more than 10 consecutive bad GRIB records")')
1797 write(*,'(" *** Let''s just stop now.",//)')
1798 write(*,'(" Perhaps the analysis file should have been converted",/,&
1799 &" from COS-Blocked format?",//)')
1800 stop
1801 endif
1802 else
1803 sevencount = 0
1804 endif
1806 ierr = 0
1808 end subroutine gribheader
1809 !
1810 !=============================================================================!
1811 !=============================================================================!
1812 !=============================================================================!
1813 !
1814 subroutine gribdata(datarray, ndat)
1816 !-----------------------------------------------------------------------------!
1817 ! !
1818 ! Read and unpack the data from a GRIB record. !
1819 ! !
1820 ! Input: !
1821 ! NDAT: The size of the data array we expect to unpack. !
1822 ! !
1823 ! Output: !
1824 ! DATARRAY: The unpacked data from the GRIB record !
1825 ! !
1826 ! Side Effects: !
1827 ! STOP if it cannot unpack the data. !
1828 ! !
1829 ! Externals: !
1830 ! SGUP_BITMAP !
1831 ! SGUP_NOBITMAP !
1832 ! CSHUP !
1833 ! !
1834 ! Modules: !
1835 ! MODULE_GRIB !
1836 ! !
1837 !-----------------------------------------------------------------------------!
1838 use module_grib
1840 implicit none
1842 integer :: ndat
1843 real , dimension(ndat) :: datarray
1844 integer, dimension(ndat) :: IX
1846 integer :: iskip, nbm
1848 if (sec4(2) == 0) then ! Grid-point data
1849 if (sec4(3).eq.0) then ! Simple unpacking
1850 if ((sec1(6).eq.64).or.(sec1(6).eq.192)) then ! There is a bitmap
1851 call SGUP_BITMAP(datarray, ndat)
1852 else
1853 call SGUP_NOBITMAP(datarray, ndat)
1854 endif
1855 else
1856 write(*,'(//,"***** No complex unpacking of gridpoint data.")')
1857 write(*,'("***** Option not yet available.",//)')
1858 ! write(*,'("***** Complain to mesouser@ucar.edu",//)')
1859 stop
1860 endif
1861 else
1862 if (sec4(3).eq.0) then ! Simple unpacking
1863 write(*,'(//,"***** No simple unpacking of spherical-harmonic coefficients.")')
1864 write(*,'("***** Option not yet available.",//)')
1865 ! write(*,'("***** Complain to mesouser@ucar.edu",//)')
1866 stop
1867 elseif (sec4(3).eq.1) then
1868 call CSHUP(datarray, ndat)
1869 endif
1870 endif
1872 end subroutine gribdata
1874 subroutine deallogrib
1875 ! Deallocates a couple of arrays that may be allocated.
1876 use module_grib
1877 #if defined (CRAY)
1878 #else
1879 if (allocated(grec)) deallocate(grec)
1880 if (allocated(bitmap)) deallocate(bitmap)
1881 #endif
1882 end subroutine deallogrib
1884 SUBROUTINE gribLGG( NLAT, startlat, deltalat )
1887 implicit none
1888 !
1889 ! LGGAUS finds the Gaussian latitudes by finding the roots of the
1890 ! ordinary Legendre polynomial of degree NLAT using Newtons
1891 ! iteration method.
1892 !
1893 ! On entry:
1894 integer NLAT ! the number of latitudes (degree of the polynomial)
1895 !
1896 ! On exit: for each Gaussian latitude
1898 double precision, dimension(NLAT) :: LATG ! Latitude
1900 ! Approximations to a regular latitude grid:
1901 real :: deltalat
1902 real :: startlat
1904 !-----------------------------------------------------------------------
1906 integer :: iskip = 15
1907 double precision :: sum1 = 0.
1908 double precision :: sum2 = 0.
1909 double precision :: sum3 = 0.
1910 double precision :: sum4 = 0.
1911 double precision :: xn
1913 integer, save :: SAVE_NLAT = -99
1914 real, save :: save_deltalat = -99.
1915 real, save :: save_startlat = -99.
1917 double precision, dimension(nlat) :: COSC, SINC
1918 double precision, parameter :: PI = 3.141592653589793
1919 !
1920 ! -convergence criterion for iteration of cos latitude
1921 double precision, parameter :: XLIM = 1.0E-14
1922 integer :: nzero, i, j
1923 double precision :: fi, fi1, a, b, g, gm, gp, gt, delta, c, d
1925 if (nlat == save_nlat) then
1926 deltalat = save_deltalat
1927 startlat = save_startlat
1928 return
1929 endif
1930 !
1931 ! -the number of zeros between pole and equator
1932 NZERO = NLAT/2
1933 !
1934 ! -set first guess for cos(colat)
1935 DO I=1,NZERO
1936 COSC(I) = SIN( (I-0.5)*PI/NLAT + PI*0.5 )
1937 ENDDO
1938 !
1939 ! -constants for determining the derivative of the polynomial
1940 FI = NLAT
1941 FI1 = FI+1.0
1942 A = FI*FI1 / SQRT(4.0*FI1*FI1-1.0)
1943 B = FI1*FI / SQRT(4.0*FI*FI-1.0)
1944 !
1945 ! -loop over latitudes, iterating the search for each root
1946 DO I=1,NZERO
1947 J=0
1948 !
1949 ! -determine the value of the ordinary Legendre polynomial for
1950 ! -the current guess root
1951 LOOP30 : DO
1952 CALL LGORD( G, COSC(I), NLAT )
1953 !
1954 ! -determine the derivative of the polynomial at this point
1955 CALL LGORD( GM, COSC(I), NLAT-1 )
1956 CALL LGORD( GP, COSC(I), NLAT+1 )
1957 GT = (COSC(I)*COSC(I)-1.0) / (A*GP-B*GM)
1958 !
1959 ! -update the estimate of the root
1960 DELTA = G*GT
1961 COSC(I) = COSC(I) - DELTA
1962 !
1963 ! -if convergence criterion has not been met, keep trying
1964 J = J+1
1965 IF( ABS(DELTA).LE.XLIM ) EXIT LOOP30
1966 ENDDO LOOP30
1967 ENDDO
1968 !
1969 ! Determine the sin(colat)
1970 SINC(1:NZERO) = SIN(ACOS(COSC(1:NZERO)))
1971 !
1972 ! -if NLAT is odd, set values at the equator
1973 IF( MOD(NLAT,2) .NE. 0 ) THEN
1974 I = NZERO+1
1975 SINC(I) = 1.0
1976 latg(i) = 0.
1977 END IF
1979 ! Set the latitudes.
1981 latg(1:NZERO) = dacos(sinc(1:NZERO)) * 180. / pi
1983 ! Determine the southern hemisphere values by symmetry
1984 do i = 1, nzero
1985 latg(nlat-nzero+i) = -latg(nzero+1-i)
1986 enddo
1989 ! Now that we have the true values, find some approximate values.
1991 xn = float(nlat-iskip*2)
1992 do i = iskip+1, nlat-iskip
1993 sum1 = sum1 + latg(i)*float(i)
1994 sum2 = sum2 + float(i)
1995 sum3 = sum3 + latg(i)
1996 sum4 = sum4 + float(i)**2
1997 enddo
1999 b = (xn*sum1 - sum2*sum3) / (xn*sum4 - sum2**2)
2000 a = (sum3 - b * sum2) / xn
2002 deltalat = sngl(b)
2003 startlat = sngl(a + b)
2005 save_nlat = nlat
2006 save_deltalat = deltalat
2007 save_startlat = startlat
2009 contains
2010 SUBROUTINE LGORD( F, COSC, N )
2011 implicit none
2012 !
2013 ! LGORD calculates the value of an ordinary Legendre polynomial at a
2014 ! latitude.
2015 !
2016 ! On entry:
2017 ! COSC - cos(colatitude)
2018 ! N - the degree of the polynomial
2019 !
2020 ! On exit:
2021 ! F - the value of the Legendre polynomial of degree N at
2022 ! latitude asin(COSC)
2023 double precision :: s1, c4, a, b, fk, f, cosc, colat, c1, fn, ang
2024 integer :: n, k
2026 !------------------------------------------------------------------------
2028 colat = acos(cosc)
2029 c1 = sqrt(2.0)
2030 do k=1,n
2031 c1 = c1 * sqrt( 1.0 - 1.0/(4*k*k) )
2032 enddo
2033 fn = n
2034 ang= fn * colat
2035 s1 = 0.0
2036 c4 = 1.0
2037 a =-1.0
2038 b = 0.0
2039 do k=0,n,2
2040 if (k.eq.n) c4 = 0.5 * c4
2041 s1 = s1 + c4 * cos(ang)
2042 a = a + 2.0
2043 b = b + 1.0
2044 fk = k
2045 ang= colat * (fn-fk-2.0)
2046 c4 = ( a * (fn-b+1.0) / ( b * (fn+fn-a) ) ) * c4
2047 enddo
2048 f = s1 * c1
2049 end subroutine lgord
2051 END SUBROUTINE GRIBLGG
2053 SUBROUTINE REORDER_IT (a, nx, ny, dx, dy, iorder)
2055 use module_debug
2057 implicit none
2058 integer :: nx, ny, iorder
2059 integer :: i, j, k, m
2060 real :: dx, dy
2061 real, dimension(nx*ny) :: a, z
2063 if (iorder .eq. 0 .and. dx .gt. 0. .and. dy .lt. 0) return
2064 k = 0
2065 call mprintf(.true.,DEBUG, &
2066 "Reordering GRIB array : dx = %f , dy = %f , iorder = %i", &
2067 f1=dx,f2=dy,i1=iorder)
2068 if (iorder .eq. 0 ) then
2069 if ( dx .lt. 0 .and. dy .lt. 0. ) then
2070 do j = 1, ny
2071 do i = nx, 1, -1
2072 k = k + 1
2073 m = i * j
2074 z(k) = a(m)
2075 enddo
2076 enddo
2077 else if ( dx .lt. 0 .and. dy .gt. 0. ) then
2078 do j = ny, 1, -1
2079 do i = nx, 1, -1
2080 k = k + 1
2081 m = i * j
2082 z(k) = a(m)
2083 enddo
2084 enddo
2085 else if ( dx .gt. 0 .and. dy .gt. 0. ) then
2086 do j = ny, 1, -1
2087 do i = 1, nx
2088 k = k + 1
2089 m = i * j
2090 z(k) = a(m)
2091 enddo
2092 enddo
2093 endif
2094 else
2095 if ( dx .gt. 0 .and. dy .lt. 0. ) then
2096 do i = 1, nx
2097 do j = 1, ny
2098 k = k + 1
2099 m = i * j
2100 z(k) = a(m)
2101 enddo
2102 enddo
2103 else if ( dx .lt. 0 .and. dy .lt. 0. ) then
2104 do i = nx, 1, -1
2105 do j = 1, ny
2106 k = k + 1
2107 m = i * j
2108 z(k) = a(m)
2109 enddo
2110 enddo
2111 else if ( dx .lt. 0 .and. dy .lt. 0. ) then
2112 do i = nx, 1, -1
2113 do j = ny, 1, -1
2114 k = k + 1
2115 m = i * j
2116 z(k) = a(m)
2117 enddo
2118 enddo
2119 else if ( dx .gt. 0 .and. dy .gt. 0. ) then
2120 do i = 1, nx
2121 do j = ny, 1, -1
2122 k = k + 1
2123 m = i * j
2124 z(k) = a(m)
2125 enddo
2126 enddo
2127 endif
2128 endif
2129 ! now put it back in the 1-d array and reset the dx and dy
2130 do k = 1, nx*ny
2131 a(k) = z(k)
2132 enddo
2133 dx = abs ( dx)
2134 dy = -1 * abs(dy)
2135 return
2136 END SUBROUTINE REORDER_IT