| blob | 26baf20db3621f36997da34a60299c335effaa27 |
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, ec_rec_len) !
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, ec_rec_len) !
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, ec_rec_len)
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, ec_rec_len
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 if ( ec_rec_len .ne. 0 ) then
299 gribsize = ec_rec_len
300 else
301 call gbyte_g1(trec, gribsize, 32, 24)
302 endif
304 ! GRIB Edition 0 does not include the total size, so we have to sum up !
305 ! the sizes of the individual sections !
307 elseif (ied.eq.0) then
308 ! Grib old edition
309 ! Find the size of section 1.
310 call gbyte_g1(trec, isz1, isz0, 24)
311 isz1 = isz1 * 8
312 call gbyte_g1(trec, iflag, isz0+56, 8)
313 if ((iflag.eq.128).or.(iflag.eq.192)) then ! section 2 is there
314 ! Find the size of section 2.
315 call gbyte_g1(trec, isz2, isz0+isz1, 24)
316 isz2 = isz2 * 8
317 endif
318 if ((iflag.eq.64).or.(iflag.eq.192)) then ! Section 3 is there
319 ! Find the size of section 3.
320 call gbyte_g1(trec, isz3, isz0+isz1+isz2, 24)
321 isz3 = isz3 * 8
322 endif
323 ! Find the size of section 4.
324 call gbyte_g1(trec, isz4, isz0+isz1+isz2+isz3, 24)
325 isz4 = isz4 * 8
327 ! Total the sizes of sections 0 through 5.
328 gribsize = (isz0+isz1+isz2+isz3+isz4+isz5) / 8
330 elseif (ied.eq.2) then
331 ! Grib2
332 CALL getarg ( 0 , pname )
333 write(*,'("*** stopping in gribcode ***\n")')
334 write(*,'("\tI was expecting a Grib1 file, but this is a Grib2 file.")')
335 if ( index(pname,'ungrib.exe') .ne. 0 ) then
336 write(*,'("\tIt is possible this is because your GRIBFILE.XXX files")')
337 write(*,'("\tare not all of the same type.")')
338 write(*,'("\tWPS can handle both file types, but a separate ungrib")')
339 write(*,'("\tjob must be run for each Grib type.\n")')
340 else
341 write(*,'("\tUse g2print on Grib2 files\n")')
342 endif
343 stop 'gribsize in gribcode'
344 else
345 write(*,'("Error trying to read grib edition number in gribsize.")')
346 write(*,'("Possible corrupt grib file.")')
347 write(6,*) 'Incorrect edition number = ',ied
348 write(6,*) 'Skipping the rest of the file and continuing.'
349 ierr = 1
350 endif
351 end function gribsize
352 !
353 !=============================================================================!
354 !=============================================================================!
355 !=============================================================================!
356 !
357 subroutine findgrib(nunit, isize, ierr, ec_rec_len)
359 !-----------------------------------------------------------------------------!
360 ! !
361 ! Find the string "GRIB", which starts off a GRIB record. !
362 ! !
363 ! Input: !
364 ! NUNIT: The C unit to read from. This should already be opened. !
365 ! !
366 ! Output: !
367 ! ISIZE: The size in bytes of one complete GRIB Record !
368 ! IERR: Error flag, !
369 ! 0 : No error or end-of-file on reading !
370 ! 1 : Hit the end of file !
371 ! 2 : Error on reading !
372 ! !
373 ! Side effects: !
374 ! * The pointer to C unit NUNIT is set to the beginning of the next !
375 ! GRIB record. !
376 ! * The first time FINDGRIB is called, the integer GTEST is set to !
377 ! a value equivalent to the string 'GRIB' !
378 ! !
379 ! Modules: !
380 ! MODULE_GRIB !
381 ! !
382 ! Externals: !
383 ! BN_READ !
384 ! BN_SEEK !
385 ! GRIBSIZE !
386 ! !
387 !-----------------------------------------------------------------------------!
388 implicit none
389 integer, intent(in) :: nunit, ec_rec_len
390 integer, intent(out) :: isize
391 integer, intent(out) :: ierr
393 integer, parameter :: LENTMP=100
394 integer, dimension(lentmp) :: trec
396 integer :: isz, itest, icnt
398 integer, save :: gtest = 0
400 ! Set the integer variable GTEST to hold the integer representation of the
401 ! character string 'GRIB'. This integer variable is later compared to
402 ! integers we read from the GRIB file, to find the beginning of a GRIB record.
404 if (gtest.eq.0) then
405 if (mwsize.eq.32) then
406 gtest = transfer('GRIB', gtest)
407 elseif(mwsize.eq.64) then
408 call gbyte_g1(char(0)//char(0)//char(0)//char(0)//'GRIB', gtest, 0, mwsize)
409 endif
410 endif
411 ierr = 0
412 icnt = 0
414 LOOP : DO
415 ! Read LENTMP bytes into holding array TREC.
416 call bn_read(nunit, trec, lentmp, isz, ierr, 0)
417 if (ierr.eq.1) then
418 return
419 elseif (ierr.eq.2) then
420 write(*,'("Error reading GRIB: IERR = ", I2)') ierr
421 return
422 endif
423 ! Reposition the file pointer back to where we started.
424 call bn_seek(nunit, -isz, 0, 0)
426 ! Compare the first four bytes of TREC with the string 'GRIB' stored in
427 ! integer variable GTEST.
428 if (mwsize.eq.32) then
429 if (trec(1) == gtest) exit LOOP
430 elseif (mwsize.eq.64) then
431 call gbyte_g1(trec, itest, 0, 32)
432 if (itest == gtest) exit LOOP
433 endif
435 ! Advance the file pointer one byte.
436 call bn_seek(nunit, 1, 0, 0)
437 icnt = icnt + 1
438 if ( icnt .gt. 100000) then ! stop if we cannot find the GRIB string
439 write(*,'("*** stopping in findgrib in gribcode ***\n")')
440 write(*,'("\tI could not find the GRIB string in the input file")')
441 write(*,'("\tafter testing the first 100,000 bytes.")')
442 write(*,'("\tThe file may be corrupt or it is not a GRIB file.")')
443 write(*,'("\tPerhaps a gzipped GRIB file or netcdf?\n")')
444 stop 'findgrib'
445 endif
447 ENDDO LOOP
449 !#if defined (DEC) || defined (ALPHA) || defined (alpha) || defined (LINUX)
450 #ifdef BYTESWAP
451 call swap4(trec, isz)
452 #endif
453 isize = gribsize(trec, isz, ierr, ec_rec_len)
455 end subroutine findgrib
456 !
457 !=============================================================================!
458 !=============================================================================!
459 !=============================================================================!
460 !
461 subroutine SGUP_NOBITMAP(datarray, ndat)
462 ! Simple grid-point unpacking
463 implicit none
465 integer :: ndat
466 real , dimension(ndat) :: datarray
467 integer, dimension(ndat) :: IX
468 real :: dfac, bfac
469 integer :: iskip
471 DFAC = 10.**(-sec1(24))
472 BFAC = 2.**sec4(7)
473 if (ied.eq.0) then
474 iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
475 elseif (ied.eq.1) then
476 iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
477 endif
478 ! sec4(8) is the number of bits used per datum value.
479 ! If sec4(8) = 255, assume they mean sec4(8) = 0
480 if (sec4(8) == 255) then
481 sec4(8) = 0
482 endif
483 ! If sec4(8) is 0, assume datarray is constant value of xec4(1)
485 if (sec4(8).eq.0) then
486 !!! HERE IS THE ORIGINAL NCAR CODE:
487 ! datarray = xec4(1)
488 !!! HERE IS WHAT FSL CHANGED IT TO:
489 datarray = DFAC*xec4(1)
490 !!! because even though it is a constant value
491 !!! you still need to scale by the decimal scale factor.
492 else
493 !!! FSL developers MOVED THE CALL TO gbytes FROM line 441 ABOVE
494 !!! BECAUSE IF sec4(8)==0 BEFORE gbytes IS CALLED, THE MASKS ARRAY
495 !!! IN gbytes WILL BE INDEXED OUT OF BOUNDS. C HARROP 9/16/04
496 call gbytes_g1(grec, IX, iskip, sec4(8), 0, ndat)
497 datarray = DFAC * (xec4(1) + (IX*BFAC))
498 endif
499 end subroutine SGUP_NOBITMAP
500 !
501 !=============================================================================!
502 !=============================================================================!
503 !=============================================================================!
504 !
506 subroutine SGUP_BITMAP(datarray, ndat)
507 ! Simple grid-point unpacking, with a bitmap.
508 implicit none
510 integer :: ndat ! Number of data points in the final grid.
511 real , dimension(ndat) :: datarray ! Array holding the final unpacked data.
512 real :: dfac, bfac
513 integer :: iskip, nbm, i, nn
515 integer, allocatable, dimension(:) :: bmdat
517 ! SEC4(1) : The number of bytes in the whole of GRIB Section 4.
518 ! SEC4(6) : The number of unused bits at the end of GRIB Section 4.
519 ! SEC4(8) : The number of bits per data value.
521 datarray = -1.E30
523 ! 1) There are fewer than NDAT data values, because a bitmap was used.
524 ! Compute the number of data values (NBM). There are 11 extra bytes
525 ! in the header section 4. NBM equals the total number of data bits (not
526 ! counting the header bits), minus the number of unused buts, and then
527 ! divided by the number of bits per data value.
529 ! Compute the parameters involved with packing
530 DFAC = 10.**(-sec1(24))
531 BFAC = 2.**sec4(7)
533 ! If sec4(8) is 0, assume datarray is constant value of xec4(1) scaled by DFAC
535 if (sec4(8).eq.0) then
536 where(bitmap(1:ndat).eq.1) datarray = xec4(1) * DFAC
537 return
538 endif
539 nbm = ((sec4(1)-11)*8-sec4(6))/sec4(8)
540 allocate(bmdat(nbm))
542 ! Set ISKIP to the beginning of the data.
543 if (ied.eq.0) then
544 iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
545 elseif (ied.eq.1) then
546 iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
547 endif
549 ! Read the data from the GREC array
550 call gbytes_g1(grec, bmdat, iskip, sec4(8), 0, nbm)
551 ! sec4(8) is the number of bits used per datum value.
552 ! If sec4(8) = 255, assume they mean sec4(8) = 0
553 if (sec4(8) == 255) sec4(8) = 0
555 ! Unpack the data according to packing parameters DFAC, BFAC, and XEC4(1),
556 ! and masked by the bitmap BITMAP.
557 nn = 0
558 do i = 1, ndat
559 if (bitmap(i).eq.1) then
560 nn = nn + 1
561 datarray(i) = DFAC * (xec4(1) + (bmdat(nn)*BFAC))
562 endif
563 enddo
565 ! Deallocate the scratch BMDAT array
566 deallocate(bmdat)
568 end subroutine SGUP_BITMAP
569 !
570 !=============================================================================!
571 !=============================================================================!
572 !=============================================================================!
573 !
575 subroutine CSHUP(pdata, ndat)
576 ! ECMWFs unpacking of ECMWFs Complex Spherical Harmonic packing
577 ! Adapted from ECMWFs GRIBEX package.
578 implicit none
580 integer :: ndat
581 real , dimension(ndat) :: pdata
582 integer, dimension(ndat+500) :: IX
584 integer :: iskip, isign
585 integer :: N1, IPOWER, J, K, M, nval
586 real :: zscale, zref
588 integer :: ic, jm, iuc, il2, inum, jn
589 integer :: inext, ilast, ioff, jrow, index, i, jcol
590 real :: bval
591 integer, allocatable, dimension(:) :: iexp, imant
592 real , dimension(0:400) :: factor
593 real :: power
594 integer :: N
596 index = -1
598 if (ied.eq.0) then
599 iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
600 elseif(ied.eq.1) then
601 iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
602 endif
604 call gbyte_g1(grec,N1,iskip,16)
605 iskip = iskip + 16
607 call gbyte_g1(grec,ipower,iskip,16)
608 iskip = iskip + 16
609 if (ipower.ge.32768) ipower = 32768-ipower
611 ! Unpack the resolution parameters for the initial (small) truncation:
612 call gbyte_g1(grec,J,iskip,8)
613 iskip = iskip + 8
614 call gbyte_g1(grec,K,iskip,8)
615 iskip = iskip + 8
616 call gbyte_g1(grec,M,iskip,8)
617 iskip = iskip + 8
619 zscale = 2.**sec4(7)
621 iskip = N1*8
623 nval = NDAT - (J+1)*(J+2)
625 call gbytes_g1(grec, ix, iskip, sec4(8), 0, nval)
626 ! sec4(8) is the number of bits used per datum value.
627 ! If sec4(8) = 255, assume they mean sec4(8) = 0
628 if (sec4(8) == 255) sec4(8) = 0
630 pdata(1:nval) = (float(ix(1:nval))*zscale)+xec4(1)
632 IUC = NDAT+1
633 IC = NVAL+1
634 DO JM=INFOGRID(1),0,-1
635 IL2=MAX(JM,J+1)
636 INUM=2*(INFOGRID(1)-IL2+1)
637 pdata(iuc-inum:iuc-1) = pdata(ic-inum:ic-1)
638 iuc = iuc - inum
639 ic = ic - inum
640 IUC = IUC-MAX((IL2-JM)*2,0)
641 ENDDO
643 if (ied.eq.0) then
644 iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
645 elseif (ied.eq.1) then
646 iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 18*8
647 endif
649 allocate(iexp(802))
650 allocate(imant(802))
651 ilast=j+1
652 do jrow=1,ilast
653 inext = 2*(ilast-jrow+1)
654 ! extract all the exponents
655 call gbytes_g1(grec, iexp, iskip, 8, 24, inext)
656 ! extract all the mantissas
657 ioff = 8
658 call gbytes_g1(grec, imant, iskip+8, 24, 8, inext)
659 iskip = iskip + inext*32
661 ! Build the real values from mantissas and exponents
662 bval = 2.**(-24)
663 i = 0
664 do jcol = jrow, infogrid(1)+1
665 index = index + 2
666 if (ilast.ge.jcol) then
667 i = i + 1
668 if ((iexp(i).eq.128.or.iexp(i).eq.0).and.(imant(i).eq.0)) then
669 pdata(i) = 0
670 else
671 if (iexp(i).ge.128) then
672 iexp(i) = iexp(i) - 128
673 isign = -1
674 else
675 isign = 1
676 endif
677 pdata(index) = isign*bval*IMANT(i)*16.**(IEXP(i)-64)
678 i = i + 1
679 if (iexp(i).ge.128) then
680 iexp(i) = iexp(i) - 128
681 isign = -1
682 else
683 isign = 1
684 endif
685 pdata(index+1) = isign*bval*IMANT(i)*16.**(IEXP(i)-64)
686 endif
687 endif
688 enddo
689 enddo
691 !Apply power scaling:
693 if (ipower.ne.0) then
694 power = float(ipower) / 1000.0
695 factor(0) = 1.0
696 do n = 1 , infogrid(1)
697 if( ipower .ne. 1000 ) then
698 factor(n) = 1.0 / (n * (n+1) )**power
699 else
700 factor(n) = 1.0 / (n * (n + 1))
701 endif
702 enddo
703 INDEX = -1
704 DO M = 0 , J-1
705 DO N = M , INFOGRID(1)
706 INDEX = INDEX + 2
707 IF ( N .GE. J ) THEN
708 PDATA(INDEX:INDEX+1) = PDATA(INDEX:INDEX+1) * FACTOR(N)
709 ENDIF
710 ENDDO
711 ENDDO
712 DO M = J , INFOGRID(1)
713 DO N = M , INFOGRID(1)
714 INDEX = INDEX + 2
715 PDATA(INDEX:INDEX+1) = PDATA(INDEX:INDEX+1) * FACTOR(N)
716 ENDDO
717 ENDDO
718 endif
720 end subroutine CSHUP
721 !
722 !=============================================================================!
723 !=============================================================================!
724 !=============================================================================!
725 !
726 !
727 !
728 ! Trigonometric functions which deal with degrees, rather than radians:
729 !
730 real function sind(theta)
731 real :: theta
732 sind = sin(theta*degran)
733 end function sind
734 real function cosd(theta)
735 real :: theta
736 cosd = cos(theta*degran)
737 end function cosd
738 real function tand(theta)
739 real :: theta
740 tand = tan(theta*degran)
741 end function tand
742 real function atand(x)
743 real :: x
744 atand = atan(x)*raddeg
745 end function atand
746 real function atan2d(x,y)
747 real :: x,y
748 atan2d = atan2(x,y)*raddeg
749 end function atan2d
750 real function asind(x)
751 real :: x
752 asind = asin(x)*raddeg
753 end function asind
754 real function acosd(x)
755 real :: x
756 acosd = acos(x)*raddeg
757 end function acosd
759 !
760 !=============================================================================!
761 !=============================================================================!
762 !=============================================================================!
763 !
764 end module module_grib
765 !
766 !=============================================================================!
767 !=============================================================================!
768 !=============================================================================!
769 !
770 subroutine gribget(nunit, ierr, ec_rec_len)
771 use module_grib
772 !-----------------------------------------------------------------------------!
773 ! !
774 ! Read a single GRIB record, with no unpacking of any header or data fields. !
775 ! !
776 ! Input: !
777 ! NUNIT: C unit number to read from. This should already be open. !
778 ! !
779 ! Output: !
780 ! IERR: Error flag, Non-zero means there was a problem with the read. !
781 ! !
782 ! Side Effects: !
783 ! The array GREC is allocated, and filled with one GRIB record. !
784 ! The C unit pointer is moved to the end of the GRIB record just read. !
785 ! !
786 ! Modules: !
787 ! MODULE_GRIB !
788 ! !
789 ! Externals: !
790 ! FINDGRIB !
791 ! BN_READ !
792 ! !
793 !-----------------------------------------------------------------------------!
795 implicit none
797 integer :: nunit, ec_rec_len
798 integer :: ierr
799 integer :: isz, isize
801 ! Position the file pointer at the beginning of the GRIB record.
802 call findgrib(nunit, isize, ierr, ec_rec_len)
803 if (ierr.ne.0) return
805 ! Allocate the GREC array to be able to hold the data
807 #if defined (CRAY)
808 #else
809 allocate(grec((isize+(mwsize/8-1))/(mwsize/8)))
810 #endif
812 ! Read the full GRIB record.
814 call bn_read(nunit, grec, isize, isz, ierr, 1)
816 !#if defined (DEC) || defined (ALPHA) || defined (alpha) || defined (LINUX)
817 #ifdef BYTESWAP
818 call swap4(grec, isz)
819 #endif
822 end subroutine gribget
823 !
824 !=============================================================================!
825 !=============================================================================!
826 !=============================================================================!
827 !
828 subroutine gribread(nunit, data, ndata, debug_level, ierr, ec_rec_len)
829 !-----------------------------------------------------------------------------!
830 ! Read one grib record, unpack the header and data information. !
831 ! !
832 ! Input: !
833 ! NUNIT: C Unit to read from. !
834 ! NDATA: Size of array DATA (Should be >= NDAT as computed herein.) !
835 ! !
836 ! Output: !
837 ! DATA: The unpacked data array !
838 ! IERR: Error flag, non-zero means there was a problem. !
839 ! !
840 ! Side Effects: !
841 ! * Header arrays SEC0, SEC1, SEC2, SEC3, SEC4, XEC4, INFOGRID and !
842 ! INFOGRID are filled. !
843 ! * The BITMAP array is filled. !
844 ! * The C unit pointer is advanced to the end of the GRIB record. !
845 ! !
846 ! Modules: !
847 ! MODULE_GRIB !
848 ! !
849 ! Externals: !
850 ! GRIBGET !
851 ! GRIBHEADER !
852 ! GRIBDATA !
853 ! !
854 !-----------------------------------------------------------------------------!
855 use module_grib
857 implicit none
859 integer :: nunit, ec_rec_len
860 integer :: debug_level
861 integer :: ierr
862 real, allocatable, dimension(:) :: datarray
863 integer :: ndata
864 real, dimension(ndata) :: data
866 integer :: ni, nj
868 ierr = 0
870 call gribget(nunit, ierr, ec_rec_len)
871 if (ierr.ne.0) return
873 ! Unpack the header information
875 call gribheader(debug_level,ierr, ec_rec_len)
877 ! Determine the size of the data array from the information in the header,
878 ! and allocate the array DATARRAY to hold that data.
880 if (sec2(4).ne.50) then
881 ni = infogrid(1)
882 nj = infogrid(2)
883 allocate(datarray(ni*nj))
884 else
885 ni = (infogrid(1)+1) * (infogrid(1)+2)
886 nj = 1
887 allocate(datarray(ni*nj))
888 endif
890 ! Unpack the data from the GRIB record, and fill the array DATARRAY.
892 call gribdata(datarray, ni*nj)
894 data(1:ni*nj) = datarray(1:ni*nj)
895 #if defined (CRAY)
896 #else
897 deallocate(grec, datarray)
898 #endif
900 end subroutine gribread
901 !
902 !=============================================================================!
903 !=============================================================================!
904 !=============================================================================!
905 !
906 subroutine get_sec1(ksec1)
907 ! Return the GRIB Section 1 header information, which has already been
908 ! unpacked by subroutine GRIBHEADER.
909 use module_grib
910 integer, dimension(100) :: ksec1
911 ksec1 = sec1
912 end subroutine get_sec1
913 !
914 !=============================================================================!
915 !=============================================================================!
916 !=============================================================================!
917 !
918 subroutine get_sec2(ksec2)
919 ! Return the GRIB Section 2 header information, which has already been
920 ! unpacked by subroutine GRIBHEADER.
921 use module_grib
922 integer, dimension(10) :: ksec2
923 ksec2 = sec2
924 end subroutine get_sec2
925 !
926 !=============================================================================!
927 !=============================================================================!
928 !=============================================================================!
929 !
930 subroutine get_gridinfo(iginfo, ginfo)
931 use module_grib
932 integer, dimension(40) :: iginfo
933 real, dimension(40) :: ginfo
934 iginfo = infogrid
935 ginfo = gridinfo
936 end subroutine get_gridinfo
937 !
938 !=============================================================================!
939 !=============================================================================!
940 !=============================================================================!
941 !
942 subroutine gribprint(isec)
943 use module_grib
944 implicit none
945 integer :: isec
946 integer :: ou = 6
947 character(len=12) :: string = ',t45,":",i8)'
948 character(len=15) :: rstring = ',t45,":",f12.5)'
950 if (isec.eq.0) then
951 write(*,'(/,"GRIB SECTION 0:")')
952 write(ou,'(5x,"Grib Length"'//string) sec0(1)
953 write(ou,'(5x,"Grib Edition"'//string) sec0(2)
954 else if (isec.eq.1) then
955 write(*,'(/,"GRIB SECTION 1:")')
956 write(ou,'(5x,"Length of PDS"'//string) sec1(1)
957 write(ou,'(5x,"Parameter Table Version"'//string) sec1(2)
958 write(ou,'(5x,"Center ID"'//string) sec1(3)
959 write(ou,'(5x,"Process ID"'//string) sec1(4)
960 write(ou,'(5x,"Grid ID"'//string) sec1(5)
961 if (sec1(25) == 1) then
962 write(ou,'(5x,"Is there a Grid Desc. Section (GDS)?",t45,": Yes")')
963 else if (sec1(25) == 0) then
964 write(ou,'(5x,"Is there a Grid Desc. Section (GDS)?",t45,": No")')
965 else
966 print*, 'Unrecognized sec1(25): ', sec1(25)
967 endif
968 if (sec1(26) == 1) then
969 write(ou,'(5x,"Is there a Bit Map Section (BMS)?",t45,": Yes")')
970 else if (sec1(26) == 0) then
971 write(ou,'(5x,"Is there a Bit Map Section (BMS)?",t45,": No")')
972 else
973 print*, 'Unrecognized sec1(26): ', sec1(26)
974 endif
975 write(ou,'(5x,"Parameter"'//string) sec1(7)
976 write(ou,'(5x,"Level type"'//string) sec1(8)
977 if ( (sec1(8) == 101) .or. (sec1(8) == 104) .or. (sec1(8) == 106) .or. &
978 (sec1(8) == 108) .or. (sec1(8) == 110) .or. (sec1(8) == 112) .or. &
979 (sec1(8) == 114) .or. (sec1(8) == 116) .or. (sec1(8) == 120) .or. &
980 (sec1(8) == 121) .or. (sec1(8) == 128) .or. (sec1(8) == 141) ) then
981 write(ou,'(5x,"Hgt, pres, etc. of layer top "'//string) sec1(9)
982 write(ou,'(5x,"Hgt, pres, etc. of layer bottom "'//string) sec1(10)
983 else
984 write(ou,'(5x,"Height, pressure, etc "'//string) sec1(9)
985 endif
986 write(ou,'(5x,"Year"'//string) sec1(11)
987 write(ou,'(5x,"Month"'//string) sec1(12)
988 write(ou,'(5x,"Day"'//string) sec1(13)
989 write(ou,'(5x,"Hour"'//string) sec1(14)
990 write(ou,'(5x,"Minute"'//string) sec1(15)
991 write(ou,'(5x,"Forecast time unit"'//string) sec1(16)
992 write(ou,'(5x,"P1"'//string) sec1(17)
993 write(ou,'(5x,"P2"'//string) sec1(18)
994 write(ou,'(5x,"Time Range Indicator"'//string) sec1(19)
995 write(ou,'(5x,"Number in Ave?"'//string) sec1(20)
996 write(ou,'(5x,"Number missing from ave?"'//string) sec1(21)
997 write(ou,'(5x,"Century"'//string) sec1(22)
998 write(ou,'(5x,"Sub-center"'//string) sec1(23)
999 write(ou,'(5x,"Decimal scale factor"'//string) sec1(24)
1000 elseif ((isec.eq.2) .and. ((sec1(6).eq.128).or.(sec1(6).eq.192))) then
1001 write(*,'(/,"GRIB SECTION 2:")')
1002 write(ou,'(5x,"Length of GRID Desc. Section"'//string) sec2(1)
1003 if ((sec2(2) /= 0).or.(sec2(3) /= 0) .or. (sec2(4) /= 0)) then
1004 write(ou,'(5x,"Number of V. Coordinate Parms"'//string) sec2(2)
1005 write(ou,'(5x,"List Starting point"'//string) sec2(3)
1006 write(ou,'(5x,"Data Representation type"'//string) sec2(4)
1007 endif
1009 if (sec2(4).eq.0) then
1010 write(ou,'(5x,"Cylindrical Equidistant Grid")')
1011 write(ou,'(10x,"NI"'//string) infogrid(1)
1012 write(ou,'(10x,"NJ"'//string) infogrid(2)
1013 write(ou,'(10x,"Lat 1"'//rstring) gridinfo(3)
1014 write(ou,'(10x,"Lon 1"'//rstring) gridinfo(4)
1015 write(ou,'(10x,"Resolution and Component:", t45,":",B8.8)') infogrid(5)
1016 write(ou,'(10x,"Lat NI"'//string) infogrid(6)
1017 write(ou,'(10x,"Lon NJ"'//string) infogrid(7)
1018 write(ou,'(10x,"Delta-Lon"'//string) infogrid(8)
1019 write(ou,'(10x,"Delta-Lat"'//string) infogrid(9)
1020 write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
1021 write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
1022 write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
1024 else if (sec2(4).eq.1) then
1025 write(ou,'(5x,"Mercator Grid")')
1026 write(ou,'(10x,"NI"'//string) infogrid(1)
1027 write(ou,'(10x,"NJ"'//string) infogrid(2)
1028 write(ou,'(10x,"Lat 1"'//rstring) gridinfo(3)
1029 write(ou,'(10x,"Lon 1"'//rstring) gridinfo(4)
1030 write(ou,'(10x,"Resolution and Component",t45,":", B8.8)') infogrid(5)
1031 write(ou,'(10x,"Lat NI"'//rstring) gridinfo(6)
1032 write(ou,'(10x,"Lon NJ"'//rstring) gridinfo(7)
1033 write(ou,'(10x,"Dx"'//rstring) gridinfo(8)
1034 write(ou,'(10x,"Dy"'//rstring) gridinfo(9)
1035 write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
1036 write(ou,'(10x,"Latin"'//rstring) gridinfo(11)
1037 write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
1038 write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
1040 else if (sec2(4).eq.4) then
1041 write(ou,'(5x,"Gaussian Grid")')
1042 write(ou,'(10x,"NI"'//string) infogrid(1)
1043 write(ou,'(10x,"NJ"'//string) infogrid(2)
1044 write(ou,'(10x,"Original (stored) Lat 1"'//rstring) gridinfo(18)
1045 write(ou,'(10x,"Lat 1"'//rstring) gridinfo(3)
1046 write(ou,'(10x,"Lon 1"'//rstring) gridinfo(4)
1047 write(ou,'(10x,"Resolution and Component",t45,":", B8.8)') infogrid(5)
1048 write(ou,'(10x,"Original (stored) Lat NI"'//rstring) gridinfo(17)
1049 write(ou,'(10x,"Lat NI"'//rstring) gridinfo(6)
1050 write(ou,'(10x,"Lon NJ"'//rstring) gridinfo(7)
1051 write(ou,'(10x,"Delta-Lon"'//rstring) gridinfo(8)
1052 write(ou,'(10x,"Delta-Lat"'//rstring) gridinfo(19)
1053 write(ou,'(10x,"Number of lats (pole - eq)"'//string) infogrid(9)
1054 write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
1055 write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
1056 write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
1057 elseif (sec2(4).eq.3) then
1058 write(ou,'(5x,"Lambert Conformal Grid")')
1059 write(ou,'(10x,"NI"'//string) infogrid(1)
1060 write(ou,'(10x,"NJ"'//string) infogrid(2)
1061 write(ou,'(10x,"Lat 1"'//string) infogrid(3)
1062 write(ou,'(10x,"Lon 1"'//string) infogrid(4)
1063 write(ou,'(10x,"Resolution and Component",t45,":", B8.8)') infogrid(5)
1064 write(ou,'(10x,"Lov"'//string) infogrid(6)
1065 write(ou,'(10x,"Dx"'//string) infogrid(7)
1066 write(ou,'(10x,"Dy"'//string) infogrid(8)
1067 write(ou,'(10x,"Projection center"'//string) infogrid(9)
1068 write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
1069 write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
1070 write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
1071 write(ou,'(10x,"Latin 1"'//string) infogrid(11)
1072 write(ou,'(10x,"Latin 2"'//string) infogrid(12)
1073 write(ou,'(10x,"Lat of southern pole"'//string) infogrid(13)
1074 write(ou,'(10x,"Lon of southern pole"'//string) infogrid(14)
1075 elseif (sec2(4).eq.5) then
1076 write(ou,'(5x,"Polar Stereographic Grid")')
1077 write(ou,'(10x,"NI"'//string) infogrid(1)
1078 write(ou,'(10x,"NJ"'//string) infogrid(2)
1079 write(ou,'(10x,"Lat 1"'//string) infogrid(3)
1080 write(ou,'(10x,"Lon 1"'//string) infogrid(4)
1081 write(ou,'(10x,"Resolution and Component", t45,":",B8.8)') infogrid(5)
1082 write(ou,'(10x,"Lov"'//string) infogrid(6)
1083 write(ou,'(10x,"Dx"'//string) infogrid(7)
1084 write(ou,'(10x,"Dy"'//string) infogrid(8)
1085 write(ou,'(10x,"Projection center"'//string) infogrid(9)
1086 write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
1087 write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
1088 write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
1089 elseif (sec2(4).eq.50) then
1090 write(ou,'(5x,"Spherical harmonic components")')
1091 write(ou,'(10x,"J-Pentagonal resolution parm:"'//string) infogrid(1)
1092 write(ou,'(10x,"K-Pentagonal resolution parm:"'//string) infogrid(2)
1093 write(ou,'(10x,"M-Pentagonal resolution parm:"'//string) infogrid(3)
1094 write(ou,'(10x,"Representation type"'//string) infogrid(4)
1095 write(ou,'(10x,"Coefficient storage mode"'//string) infogrid(5)
1096 endif
1097 elseif ((isec.eq.3) .and. (sec1(26).eq.1)) then
1098 write(*,'(/,"GRIB SECTION 3:")')
1099 write(ou,'(5x,"Length of bit map section"'//string) sec3(1)
1100 write(ou,'(5x,"Number of unused bits"'//string) sec3(2)
1101 write(ou,'(5x,"Numeric"'//string) sec3(3)
1103 elseif (isec.eq.4) then
1104 write(*,'(/,"GRIB SECTION 4:")')
1105 write(ou,'(5x,"Length of BDS"'//string) sec4(1)
1106 write(ou,'(5x,"0/1: grid-point or sph. harm. data"'//string) sec4(2)
1107 write(ou,'(5x,"0/1: simple or complex packing"'//string) sec4(3)
1108 write(ou,'(5x,"0/1: floating or integer"'//string) sec4(4)
1109 write(ou,'(5x,"0/1: No addl flags or addl flags"'//string) sec4(5)
1110 write(ou,'(5x,"Unused bits"'//string) sec4(6)
1111 write(ou,'(5x,"Binary Scale Factor"'//string) sec4(7)
1112 write(ou,'(5x,"Reference Value", t45, ":", F18.8)') xec4(1)
1113 write(ou,'(5x,"Number of bits for packing"'//string) sec4(8)
1114 endif
1116 end subroutine gribprint
1117 !
1118 !=============================================================================!
1119 !=============================================================================!
1120 !=============================================================================!
1121 !
1122 subroutine get_bitmap(bm8, ndat)
1123 use module_grib
1124 integer, dimension(ndat) :: bm8
1125 if ((sec1(6).eq.64).or.(sec1(6).eq.192)) then
1126 bm8 = bitmap
1127 else
1128 bm8 = 1
1129 endif
1130 end subroutine get_bitmap
1131 !
1132 !=============================================================================!
1133 !=============================================================================!
1134 !=============================================================================!
1135 !
1136 subroutine gribxyll(x, y, xlat, xlon)
1137 use module_grib
1138 implicit none
1140 real , intent(in) :: x, y
1141 real , intent(out) :: xlat, xlon
1143 real :: r, xkm, ykm, y1
1144 integer :: iscan, jscan
1146 if (sec2(4).eq.0) then ! Cylindrical equidistant grid
1148 xlat = gridinfo(3) + gridinfo(9)*(y-1.)
1149 xlon = gridinfo(4) + gridinfo(8)*(x-1.)
1151 elseif (sec2(4) == 1) then ! Mercator grid
1152 r = grrth*cosd(gtrue1)
1153 xkm = (x-1.)*gridinfo(8)
1154 ykm = (y-1.)*gridinfo(9)
1155 xlon = gridinfo(4) + (xkm/r)*(180./pi)
1156 y1 = r*alog((1.+sind(gridinfo(3)))/cosd(gridinfo(3)))/gridinfo(9)
1157 xlat = 90. - 2. * atan(exp(-gridinfo(9)*(y+y1-1.)/r))*180./pi
1159 elseif (sec2(4) == 3) then ! Lambert Conformal grid
1160 gclon = gridinfo(6)
1161 r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
1162 xlat = 90. - 2.*atand(tand(45.-gtrue1/2.)* &
1163 ((r*gkappa*gridinfo(7))/(grrth*sind(90.-gtrue1)))**(1./gkappa))
1164 xlon = atan2d((x-1.+gx1),-(y-1.+gy1))/gkappa + gclon
1166 elseif (sec2(4) == 5) then ! Polar Stereographic grid
1167 gclon = gridinfo(6)
1168 r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
1169 xlat = 90. - 2.*atan2d((r*gridinfo(7)),(grrth*(1.+sind(gtrue1))))
1170 xlon = atan2d((x-1.+gx1),-(y-1.+gy1)) + gclon
1172 elseif (sec2(4) == 4) then ! Gaussian grid
1174 xlon = gridinfo(4) + gridinfo(8)*(x-1.)
1175 xlat = gridinfo(3) + gridinfo(19)*(y-1.)
1177 else
1178 write(*,'("Unrecognized projection:", I10)') sec2(4)
1179 write(*,'("STOP in GRIBXYLL")')
1180 stop
1181 endif
1183 end subroutine gribxyll
1184 !
1185 !=============================================================================!
1186 !=============================================================================!
1187 !=============================================================================!
1188 !
1189 subroutine gribllxy(xlat, xlon, x, y)
1190 use module_grib
1191 implicit none
1192 real , intent(in) :: xlat, xlon
1193 real , intent(out) :: x, y
1195 real :: r, y1
1197 if (sec2(4) == 0) then ! Cylindrical Equidistant grid
1199 x = 1. + (xlon-gridinfo(4)) / gridinfo(9)
1200 y = 1. + (xlat-gridinfo(3)) / gridinfo(8)
1202 else if (sec2(4) == 1) then ! Mercator grid
1204 r = grrth*cosd(gtrue1)
1205 x = 1.+( (r/gridinfo(8)) * (xlon-gridinfo(4)) * (pi/180.) )
1206 y1 = (r/gridinfo(9))*alog((1.+sind(gridinfo(3)))/cosd(gridinfo(3)))
1207 y = 1. + ((r/gridinfo(9))*alog((1.+sind(xlat))/cosd(xlat)))-y1
1209 else if (sec2(4) == 3) then ! Lambert Conformal grid
1210 gclon = gridinfo(6)
1211 r = grrth/(gridinfo(7)*gkappa)*sind(90.-gtrue1) * &
1212 (tand(45.-xlat/2.)/tand(45.-gtrue1/2.)) ** gkappa
1213 x = r*sind(gkappa*(xlon-gclon)) - gx1 + 1.
1214 y = -r*cosd(gkappa*(xlon-gclon)) - gy1 + 1.
1216 elseif (sec2(4) == 5) then ! Polar Stereographic grid
1217 gclon = gridinfo(6)
1218 r = grrth/gridinfo(7) * tand((90.-xlat)/2.) * (1.+sind(gtrue1))
1219 x = ( r * sind(xlon-gclon)) - gx1 + 1.
1220 y = (-r * cosd(xlon-gclon)) - gy1 + 1.
1222 else
1223 write(*,'("Unrecognized projection:", I10)') sec2(4)
1224 write(*,'("STOP in GRIBLLXY")')
1225 stop
1226 endif
1228 end subroutine gribllxy
1229 !
1230 !=============================================================================!
1231 !=============================================================================!
1232 !=============================================================================!
1233 !
1234 subroutine glccone (fsplat,ssplat,sign,confac)
1235 use module_grib
1236 implicit none
1237 real, intent(in) :: fsplat,ssplat
1238 integer, intent(in) :: sign
1239 real, intent(out) :: confac
1240 if (abs(fsplat-ssplat).lt.1.E-3) then
1241 confac = sind(fsplat)
1242 else
1243 confac = log10(cosd(fsplat))-log10(cosd(ssplat))
1244 confac = confac/(log10(tand(45.-float(sign)*fsplat/2.))- &
1245 log10(tand(45.-float(sign)*ssplat/2.)))
1246 endif
1247 end subroutine glccone
1248 !
1249 !=============================================================================!
1250 !=============================================================================!
1251 !=============================================================================!
1252 !
1253 !=============================================================================!
1254 !=============================================================================!
1255 !=============================================================================!
1256 !
1257 subroutine gribheader(debug_level,ierr, ec_rec_len)
1258 !
1259 ! IERR non-zero means there was a problem unpacking the grib header.
1260 !
1261 use module_grib
1262 implicit none
1263 integer :: debug_level
1264 integer :: ierr, ec_rec_len
1266 integer, parameter :: nsec1 = 24
1268 integer, dimension(nsec1) :: &
1269 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/)
1270 integer :: icount, iskip, ibts, nbm, nbz, i9skip, i17skip
1272 integer :: iman, ichar, isign, iscan
1274 integer, allocatable, dimension(:) :: bm8
1276 real :: r
1277 integer :: isvns
1278 integer :: gsvns = 0
1280 if (gsvns.eq.0) then
1281 if (mwsize.eq.32) then
1282 gsvns = transfer('7777', gsvns)
1283 elseif(mwsize.eq.64) then
1284 call gbyte_g1(char(0)//char(0)//char(0)//char(0)//'7777', gsvns, 0, mwsize)
1285 endif
1286 endif
1288 ! Section 0:
1289 sec0(2) = ied
1290 if (ied.eq.1) then
1291 if ( ec_rec_len .ne. 0 ) then
1292 sec0(1) = ec_rec_len
1293 else
1294 call gbyte_g1(grec, sec0(1), 32, 24)
1295 endif
1296 iskip = 64
1297 elseif (ied.eq.0) then
1298 sec0(1) = gribsize(grec,200, ierr, 0)
1299 iskip = 32
1300 endif
1302 ! Section 1:
1303 i9skip = iskip + 80
1304 i17skip = iskip + 144
1305 do icount = 1, nsec1 - ((1-ied)*6)
1306 ibts = iw1(icount)*8
1307 call gbyte_g1(grec, sec1(icount), iskip, ibts)
1308 iskip = iskip + ibts
1309 enddo
1310 if (ied.eq.0) sec1(22) = 20
1311 ! Sec1 indices 9 and 10 might actually be one value, not two.
1312 ! If this is the case, reread sec1(9), and set sec1(10) to zero:
1313 if ( (sec1(8) == 101) .or. (sec1(8) == 104) .or. (sec1(8) == 106) .or. &
1314 (sec1(8) == 108) .or. (sec1(8) == 110) .or. (sec1(8) == 112) .or. &
1315 (sec1(8) == 114) .or. (sec1(8) == 116) .or. (sec1(8) == 120) .or. &
1316 (sec1(8) == 121) .or. (sec1(8) == 128) .or. (sec1(8) == 141) .or. &
1317 (sec1(8) == 236) ) then
1318 ! No action here.
1319 else
1320 call gbyte_g1(grec, sec1(9), i9skip, 16)
1321 sec1(10) = 0.
1322 endif
1324 if (sec1(24).ge.32768) sec1(24) = 32768-sec1(24)
1326 ! If the TIME/RANGE INDICATOR (sec1(19)) indicates that the time P1
1327 ! is spread over two bytes, then recompute sec1(17) and set sec1(18)
1328 ! to zero.
1329 if (sec1(19) == 10) then
1330 call gbyte_g1(grec, sec1(17), i17skip, 16)
1331 sec1(18) = 0
1332 endif
1334 ! Pull out single bits from sec1(6) for the GDS and BMS flags:
1335 sec1(25) = sec1(6)/128
1336 sec1(26) = mod(sec1(6)/64,2)
1338 ! Section 2:
1339 ! if ((sec1(6) == 128) .or. (sec1(6) == 192)) then
1340 if (sec1(25) == 1) then
1342 if (ied.eq.0) then
1343 iskip = 32 + sec1(1)*8
1344 elseif (ied.eq.1) then
1345 iskip = 64 + sec1(1)*8
1346 endif
1347 call gbyte_g1(grec, sec2(1), iskip, 24)
1348 iskip = iskip + 24
1349 call gbytes_g1(grec, sec2(2), iskip, 8, 0, 3)
1350 iskip = iskip + 8*3
1352 if (sec2(4) == 0) then
1353 ! Lat/Lon Grid:
1354 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)
1355 iskip = iskip + 32
1356 call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
1357 iskip = iskip + 48
1358 call gbyte_g1(grec, infogrid(5), iskip, 8)
1359 iskip = iskip + 8
1360 call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 2)
1361 iskip = iskip + 48
1362 call gbytes_g1(grec, infogrid(8), iskip, 16, 0, 2)
1363 iskip = iskip + 32
1364 call gbyte_g1(grec, infogrid(21), iskip, 1)
1365 infogrid(21) = 1-(infogrid(21)*2)
1366 iskip = iskip + 1
1367 call gbyte_g1(grec, infogrid(22), iskip, 1)
1368 infogrid(22) = (infogrid(22)*2)-1
1369 iskip = iskip + 1
1370 call gbyte_g1(grec, infogrid(10), iskip, 1)
1371 iskip = iskip + 1
1372 iskip = iskip + 5
1373 call gbyte_g1(grec, infogrid(11), iskip, 32)
1374 iskip = iskip + 32
1376 !MGD if ( debug_level .gt. 100 ) then
1377 !MGD print *, "lat/lon grib grid info", infogrid(1), infogrid(3), &
1378 !MGD infogrid(5), infogrid(6), infogrid(8), infogrid(21), &
1379 !MGD infogrid(22), infogrid(10), infogrid(11), infogrid(8)
1380 !MGD end if
1382 infogrid(8) = infogrid(8) * infogrid(21)
1383 infogrid(9) = infogrid(9) * infogrid(22)
1385 gridinfo(1) = float(infogrid(1))
1386 gridinfo(2) = float(infogrid(2))
1387 if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
1388 if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
1389 gridinfo(3) = float(infogrid(3))*0.001
1390 gridinfo(4) = infogrid(4) * 0.001
1391 if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)
1392 if (infogrid(7).ge.8388608) infogrid(7) = 8388608 - infogrid(7)
1393 gridinfo(6) = infogrid(6) * 0.001
1394 gridinfo(7) = infogrid(7) * 0.001
1395 gridinfo(8) = infogrid(8) * 0.001
1396 gridinfo(9) = infogrid(9) * 0.001
1397 gridinfo(21) = float(infogrid(21))
1398 gridinfo(22) = float(infogrid(22))
1399 elseif (sec2(4) == 1) then ! Mercator grid
1400 ! Number of points in X and Y
1401 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)
1402 iskip = iskip + 32
1403 ! Starting lat and lon
1404 call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
1405 iskip = iskip + 48
1406 ! "Resolution and component flags"
1407 call gbyte_g1(grec, infogrid(5), iskip, 8)
1408 iskip = iskip + 8
1409 ! Ending lat and lon
1410 call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 2)
1411 iskip = iskip + 48
1412 ! Truelat, 3 bytes
1413 call gbyte_g1(grec, infogrid(11), iskip, 24)
1414 iskip = iskip + 24
1415 ! "Reserved", i.e., skip a byte
1416 iskip = iskip + 8
1417 ! Scanning mode flags, first three bits of the next byte
1418 ! and skip the last five bits.
1419 call gbyte_g1(grec, infogrid(21), iskip, 1)
1420 infogrid(21) = 1-(infogrid(21)*2)
1421 iskip = iskip + 1
1422 call gbyte_g1(grec, infogrid(22), iskip, 1)
1423 infogrid(22) = (infogrid(22)*2)-1
1424 iskip = iskip + 1
1425 call gbyte_g1(grec, infogrid(10), iskip, 1)
1426 iskip = iskip + 1
1427 iskip = iskip + 5
1428 ! Grid increment in X and Y
1429 call gbytes_g1(grec, infogrid(8), iskip, 24, 0, 2)
1430 iskip = iskip + 48
1431 ! Done reading map specifications.
1432 ! Now do various conversions:
1434 gridinfo(1) = float(infogrid(1)) ! ok
1435 gridinfo(2) = float(infogrid(2)) ! ok
1437 if (infogrid(3) .ge.8388608) infogrid(3) = 8388608 - infogrid(3)
1438 if (infogrid(4) .ge.8388608) infogrid(4) = 8388608 - infogrid(4)
1439 if (infogrid(6) .ge.8388608) infogrid(6) = 8388608 - infogrid(6)
1440 if (infogrid(7) .ge.8388608) infogrid(7) = 8388608 - infogrid(7)
1441 if (infogrid(11).ge.8388608) infogrid(11) = 8388608 - infogrid(11)
1442 gridinfo(3) = infogrid(3) * 0.001
1443 gridinfo(4) = infogrid(4) * 0.001
1444 gridinfo(6) = infogrid(6) * 0.001
1445 gridinfo(7) = infogrid(7) * 0.001
1446 gridinfo(8) = infogrid(8) * 0.001
1447 gridinfo(9) = infogrid(9) * 0.001
1448 gridinfo(11) = infogrid(11) * 0.001
1450 gridinfo(21) = infogrid(21)
1451 gridinfo(22) = infogrid(22)
1453 gridinfo(20) = 6370.949
1454 grrth = gridinfo(20)
1455 gtrue1 = gridinfo(11)
1457 elseif (sec2(4) == 3) then
1458 if (ied.eq.0) then
1459 print '(//,"*** Despair ***"//)'
1460 stop
1461 endif
1462 ! Lambert Conformal:
1463 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)
1464 iskip = iskip + 32
1465 call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
1466 iskip = iskip + 48
1467 if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
1468 if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
1469 call gbyte_g1(grec, infogrid(5), iskip, 8)
1470 iskip = iskip + 8
1471 call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 3)
1472 if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)
1473 iskip = iskip + 72
1474 call gbyte_g1(grec, infogrid(9), iskip, 8)
1475 iskip = iskip + 8
1476 call gbyte_g1(grec, infogrid(21), iskip, 1)
1477 infogrid(21) = 1-(infogrid(21)*2)
1478 iskip = iskip + 1
1479 call gbyte_g1(grec, infogrid(22), iskip, 1)
1480 infogrid(22) = (infogrid(22)*2)-1
1481 iskip = iskip + 1
1482 call gbyte_g1(grec, infogrid(10), iskip, 1)
1483 iskip = iskip + 1
1484 iskip = iskip + 5
1485 call gbytes_g1(grec, infogrid(11), iskip, 24, 0, 4)
1486 if (infogrid(11).ge.8388608) infogrid(11) = 8388608 - infogrid(11)
1487 if (infogrid(12).ge.8388608) infogrid(12) = 8388608 - infogrid(12)
1488 if (infogrid(13).ge.8388608) infogrid(13) = 8388608 - infogrid(13)
1489 if (infogrid(14).ge.8388608) infogrid(14) = 8388608 - infogrid(14)
1490 iskip = iskip + 96
1491 call gbyte_g1(grec, infogrid(15), iskip, 16)
1492 iskip = iskip + 16
1494 infogrid(7) = infogrid(7) * infogrid(21)
1495 infogrid(8) = infogrid(8) * infogrid(22)
1498 gridinfo(1) = float(infogrid(1))
1499 gridinfo(2) = float(infogrid(2))
1500 gridinfo(3) = infogrid(3) * 0.001
1501 gridinfo(4) = infogrid(4) * 0.001
1502 gridinfo(6) = infogrid(6) * 0.001
1503 gridinfo(7) = infogrid(7) * 0.001
1504 gridinfo(8) = infogrid(8) * 0.001
1505 gridinfo(9) = infogrid(9) * 0.001
1506 gridinfo(11) = infogrid(11) * 0.001
1507 gridinfo(12) = infogrid(12) * 0.001
1508 gridinfo(13) = infogrid(13) * 0.001
1509 gridinfo(14) = infogrid(14) * 0.001
1511 gridinfo(20) = 6370
1512 ! a priori knowledge:
1513 if (sec1(5).eq.212) then
1514 gridinfo(3) = 12.190
1515 gridinfo(4) = -133.459
1516 gridinfo(7) = 40.63525
1517 gridinfo(8) = 40.63525
1518 gridinfo(20) = 6370
1519 endif
1521 !=============================================================================!
1522 ! More a priori knowledge: !
1523 ! Correct some bad lat/lon numbers coded into some RUC headers. !
1524 ! !
1525 if (sec1(3) == 59) then ! If FSL
1526 if (sec1(4) == 86) then ! and RUC
1527 if (sec1(5) == 255) then
1528 ! Check to correct bad lat/lon numbers.
1529 if (infogrid(3) == 16909) then
1530 infogrid(3) = 16281
1531 gridinfo(3) = 16.281
1532 endif
1533 if (infogrid(4) == 236809) then
1534 infogrid(4) = 2338622
1535 gridinfo(4) = 233.8622
1536 endif
1537 endif
1538 endif
1539 endif
1540 !=============================================================================!
1543 gridinfo(21) = float(infogrid(21))
1544 gridinfo(22) = float(infogrid(22))
1546 ! Map parameters
1547 glat1 = gridinfo(3)
1548 glon1 = gridinfo(4)
1549 gclon = gridinfo(6)
1550 if (gclon.gt.180.) gclon = -(360.-gclon)
1551 if ((gclon<0).and.(glon1>180)) glon1 = glon1-360.
1552 gtrue1 = gridinfo(11)
1553 gtrue2 = gridinfo(12)
1554 grrth = gridinfo(20)
1555 call glccone(gtrue1, gtrue2, 1, gkappa)
1556 r = grrth/(gridinfo(7)*gkappa)*sind(90.-gtrue1) * &
1557 (tand(45.-glat1/2.)/tand(45.-gtrue1/2.)) ** gkappa
1558 gx1 = r*sind(gkappa*(glon1-gclon))
1559 gy1 = -r*cosd(gkappa*(glon1-gclon))
1561 elseif (sec2(4) == 4) then
1562 ! Gaussian Grid:
1563 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)
1564 iskip = iskip + 32
1565 call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
1566 iskip = iskip + 48
1567 call gbyte_g1(grec, infogrid(5), iskip, 8)
1568 iskip = iskip + 8
1569 call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 2)
1570 iskip = iskip + 48
1571 call gbytes_g1(grec, infogrid(8), iskip, 16, 0, 2)
1572 iskip = iskip + 32
1573 call gbyte_g1(grec, infogrid(21), iskip, 1)
1574 infogrid(21) = 1-(infogrid(21)*2)
1575 iskip = iskip + 1
1576 call gbyte_g1(grec, infogrid(22), iskip, 1)
1577 infogrid(22) = (infogrid(22)*2)-1
1578 iskip = iskip + 1
1579 call gbyte_g1(grec, infogrid(10), iskip, 1)
1580 iskip = iskip + 1
1581 iskip = iskip + 5
1582 call gbyte_g1(grec, infogrid(11), iskip, 32)
1583 iskip = iskip + 32
1585 infogrid(8) = infogrid(8) * infogrid(21)
1587 gridinfo(1) = float(infogrid(1))
1588 gridinfo(2) = float(infogrid(2))
1589 if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
1590 if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
1591 gridinfo(3) = float(infogrid(3))*0.001
1592 gridinfo(4) = infogrid(4) * 0.001
1593 if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)
1594 if (infogrid(7).ge.8388608) infogrid(7) = 8388608 - infogrid(7)
1595 gridinfo(6) = infogrid(6) * 0.001
1596 gridinfo(7) = infogrid(7) * 0.001
1597 gridinfo(8) = infogrid(8) * 0.001
1598 gridinfo(21) = float(infogrid(21))
1599 gridinfo(22) = float(infogrid(22))
1601 ! Compute an approximate delta-latitude and starting latitude.
1602 ! Replace the stored value of starting latitude with approximate one.
1603 gridinfo(18) = gridinfo(3)
1604 infogrid(18) = infogrid(3)
1605 gridinfo(17) = gridinfo(6)
1606 infogrid(17) = infogrid(6)
1607 ! call griblgg(infogrid(2), gridinfo(3), gridinfo(19))
1608 ! infogrid(19) = nint(gridinfo(19)*1000.)
1609 ! infogrid(3) = nint(gridinfo(3)*1000.)
1610 gridinfo(6) = -gridinfo(3)
1611 infogrid(6) = -infogrid(3)
1613 elseif (sec2(4) == 5) then
1614 ! Polar Stereographic Grid
1615 if (ied.eq.0) then
1616 print '(//,"*** Despair ***"//)'
1617 stop
1618 endif
1619 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2) ! NX and NY
1620 iskip = iskip + 32
1621 call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2) ! LAT1 and LON1
1622 iskip = iskip + 48
1623 call gbyte_g1(grec, infogrid(5), iskip, 8) ! Resolution and Component
1624 iskip = iskip + 8
1625 call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 3) ! LOV, DX, and DY
1626 iskip = iskip + 72
1627 call gbyte_g1(grec, infogrid(9), iskip, 8) ! Projection center flag
1628 iskip = iskip + 8
1629 call gbyte_g1(grec, infogrid(21), iskip, 1)
1630 infogrid(21) = 1-(infogrid(21)*2)
1631 iskip = iskip + 1
1632 call gbyte_g1(grec, infogrid(22), iskip, 1)
1633 infogrid(22) = (infogrid(22)*2)-1
1634 iskip = iskip + 1
1635 call gbyte_g1(grec, infogrid(10), iskip, 1)
1636 iskip = iskip + 1
1637 iskip = iskip + 5
1638 ! call gbyte_g1(grec, infogrid(11), iskip, 32) ! Set to 0 (reserved)
1639 iskip = iskip + 32
1641 if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
1642 if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
1643 if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)
1646 infogrid(7) = infogrid(7) * infogrid(21)
1647 infogrid(8) = infogrid(8) * infogrid(22)
1649 gridinfo(1) = float(infogrid(1))
1650 gridinfo(2) = float(infogrid(2))
1651 gridinfo(3) = infogrid(3) * 0.001
1652 gridinfo(4) = infogrid(4) * 0.001
1653 gridinfo(6) = infogrid(6) * 0.001
1654 gridinfo(7) = infogrid(7) * 0.001
1655 gridinfo(8) = infogrid(8) * 0.001
1657 gridinfo(20) = 6370
1659 ! a priori knowledge:
1660 if (sec1(5).eq.240) then
1661 gridinfo(3) = 22.7736
1662 gridinfo(4) = -120.376
1663 gridinfo(7) = 4.7625
1664 gridinfo(8) = 4.7625
1665 gridinfo(20) = 6370
1666 endif
1668 ! Map parameters
1669 glat1 = gridinfo(3)
1670 glon1 = gridinfo(4)
1671 gclon = gridinfo(6)
1672 if (gclon.gt.180.) gclon = -(360.-gclon)
1673 ! GRIB edition 1 Polar Stereographic grids are true at 60 degrees
1674 ! Which hemisphere depends on infogrid(9), the "Projection Center Flag"
1675 grrth = gridinfo(20)
1676 if (infogrid(9) > 127) then
1677 gtrue1 = -60.
1678 r = grrth/gridinfo(7) * tand((-90.-glat1)/2.) * (1.+sind(-gtrue1))
1679 gx1 = -r * sind(glon1-gridinfo(6))
1680 gy1 = -r * cosd(glon1-gridinfo(6))
1681 else
1682 gtrue1 = 60.
1683 r = grrth/gridinfo(7) * tand((90.-glat1)/2.) * (1.+sind(gtrue1))
1684 gx1 = r * sind(glon1-gridinfo(6))
1685 gy1 = -r * cosd(glon1-gridinfo(6))
1686 endif
1688 gridinfo(21) = float(infogrid(21))
1689 gridinfo(22) = float(infogrid(22))
1691 elseif (sec2(4) == 50) then
1692 ! Spherical harmonic coefficients
1693 if (ied.eq.0) then
1694 print '(//,"*** Despair ***"//)'
1695 stop
1696 endif
1697 call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 3)
1698 iskip = iskip + 48
1699 call gbytes_g1(grec, infogrid(4), iskip, 8, 0, 2)
1700 iskip = iskip + 16
1702 iskip = iskip + 18*8
1704 else
1705 call gribprint(0)
1706 call gribprint(1)
1707 call gribprint(2)
1708 call gribprint(3)
1709 call gribprint(4)
1710 write(*,'("Unrecognized grid: ", i8)') sec2(4)
1711 write(*,'("This grid is not currently supported.")')
1712 write(*,'("Write your own program to put the data to the intermediate format")')
1713 stop
1714 endif
1716 endif
1718 ! Section 3
1719 if ((sec1(6).eq.64).or.(sec1(6).eq.192)) then
1720 if (ied.eq.0) then
1721 print '(//,"*** Despair ***"//)'
1722 stop
1723 endif
1725 if (ied.eq.0) then
1726 iskip = 32 + sec1(1)*8 + sec2(1)*8
1727 elseif (ied.eq.1) then
1728 iskip = 64 + sec1(1)*8 + sec2(1)*8
1729 endif
1730 call gbyte_g1(grec, sec3(1), iskip, 24)
1731 iskip = iskip + 24
1732 call gbyte_g1(grec, sec3(2), iskip, 8)
1733 iskip = iskip + 8
1734 call gbyte_g1(grec, sec3(3), iskip, 16)
1735 iskip = iskip + 16
1737 #if defined (CRAY)
1738 #else
1739 allocate(bitmap((sec3(1)-6)*8))
1740 #endif
1741 allocate(bm8((sec3(1)-6)*8))
1742 call gbytes_g1(grec, bm8, iskip, 1, 0, (sec3(1)-6)*8)
1743 bitmap(1:size(bm8)) = bm8(1:size(bm8))
1744 deallocate(bm8)
1745 iskip = iskip + sec3(1)-6
1746 else
1747 sec3 = 0
1748 endif
1750 ! Section 4
1751 if ((sec1(6).eq.128).or.(sec1(6).eq.192)) then
1752 if (ied.eq.0) then
1753 iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8
1754 elseif (ied.eq.1) then
1755 iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8
1756 endif
1757 call gbyte_g1(grec, sec4(1), iskip, 24)
1758 if (sec4(1) > (sec0(1) - sec1(1) - sec2(1) - sec3(1) - 4)) then
1759 write(*,'(/,"*** I have good reason to believe that this GRIB record is")')
1760 write(*,'("*** corrupted or miscoded.",/)')
1761 ierr = 1
1762 return
1763 endif
1764 iskip = iskip + 24
1765 call gbytes_g1(grec, sec4(2), iskip, 1,0,4)
1766 iskip = iskip + 4
1767 call gbyte_g1(grec, sec4(6), iskip, 4)
1768 iskip = iskip + 4
1769 ! Get the binary scale factor
1770 call gbyte_g1(grec, isign, iskip, 1)
1771 iskip = iskip + 1
1772 call gbyte_g1(grec, sec4(7), iskip, 15)
1773 iskip = iskip + 15
1774 sec4(7) = sec4(7) * (-2*isign+1)
1775 ! Get the reference value:
1776 call gbyte_g1(grec, isign, iskip, 1)
1777 iskip = iskip + 1
1778 isign = -2*isign+1
1779 call gbyte_g1(grec, ichar, iskip, 7)
1780 iskip = iskip + 7
1781 call gbyte_g1(grec, iman, iskip, 24)
1782 iskip = iskip + 24
1783 if ( iman .ne. 0 ) then
1784 xec4(1) = float(isign) * (2.**(-24)) * float(iman) * &
1785 (16.**(ichar-64))
1786 else
1787 xec4(1) = 0.
1788 endif
1790 call gbyte_g1(grec,sec4(8), iskip, 8)
1791 ! sec4(8) is the number of bits used per datum value.
1792 ! If sec4(8) = 255, assume they mean sec4(8) = 0
1793 if (sec4(8) == 255) sec4(8) = 0
1794 iskip = iskip + 8
1795 endif
1797 ! Section 5
1798 call gbyte_g1(grec, isvns, ((sec0(1)-4)*8), 32)
1799 if (isvns.ne.gsvns) then
1800 write(*, '("End-of-record mark (7777) not found", 2I10)') isvns
1801 write(*, '("Sec0(1) = ", I8, i2)') sec0(1), sevencount
1802 sevencount = sevencount + 1
1803 if (sevencount > 10) then
1804 write(*,'(//," *** Found more than 10 consecutive bad GRIB records")')
1805 write(*,'(" *** Let''s just stop now.",//)')
1806 write(*,'(" Perhaps the analysis file should have been converted",/,&
1807 &" from COS-Blocked format?",//)')
1808 stop
1809 endif
1810 else
1811 sevencount = 0
1812 endif
1814 ierr = 0
1816 end subroutine gribheader
1817 !
1818 !=============================================================================!
1819 !=============================================================================!
1820 !=============================================================================!
1821 !
1822 subroutine gribdata(datarray, ndat)
1824 !-----------------------------------------------------------------------------!
1825 ! !
1826 ! Read and unpack the data from a GRIB record. !
1827 ! !
1828 ! Input: !
1829 ! NDAT: The size of the data array we expect to unpack. !
1830 ! !
1831 ! Output: !
1832 ! DATARRAY: The unpacked data from the GRIB record !
1833 ! !
1834 ! Side Effects: !
1835 ! STOP if it cannot unpack the data. !
1836 ! !
1837 ! Externals: !
1838 ! SGUP_BITMAP !
1839 ! SGUP_NOBITMAP !
1840 ! CSHUP !
1841 ! !
1842 ! Modules: !
1843 ! MODULE_GRIB !
1844 ! !
1845 !-----------------------------------------------------------------------------!
1846 use module_grib
1848 implicit none
1850 integer :: ndat
1851 real , dimension(ndat) :: datarray
1852 integer, dimension(ndat) :: IX
1854 integer :: iskip, nbm
1856 if (sec4(2) == 0) then ! Grid-point data
1857 if (sec4(3).eq.0) then ! Simple unpacking
1858 if ((sec1(6).eq.64).or.(sec1(6).eq.192)) then ! There is a bitmap
1859 call SGUP_BITMAP(datarray, ndat)
1860 else
1861 call SGUP_NOBITMAP(datarray, ndat)
1862 endif
1863 else
1864 write(*,'(//,"***** No complex unpacking of gridpoint data.")')
1865 write(*,'("***** Option not yet available.",//)')
1866 ! write(*,'("***** Complain to mesouser@ucar.edu",//)')
1867 stop
1868 endif
1869 else
1870 if (sec4(3).eq.0) then ! Simple unpacking
1871 write(*,'(//,"***** No simple unpacking of spherical-harmonic coefficients.")')
1872 write(*,'("***** Option not yet available.",//)')
1873 ! write(*,'("***** Complain to mesouser@ucar.edu",//)')
1874 stop
1875 elseif (sec4(3).eq.1) then
1876 call CSHUP(datarray, ndat)
1877 endif
1878 endif
1880 end subroutine gribdata
1882 subroutine deallogrib
1883 ! Deallocates a couple of arrays that may be allocated.
1884 use module_grib
1885 #if defined (CRAY)
1886 #else
1887 if (allocated(grec)) deallocate(grec)
1888 if (allocated(bitmap)) deallocate(bitmap)
1889 #endif
1890 end subroutine deallogrib
1892 SUBROUTINE gribLGG( NLAT, startlat, deltalat )
1895 implicit none
1896 !
1897 ! LGGAUS finds the Gaussian latitudes by finding the roots of the
1898 ! ordinary Legendre polynomial of degree NLAT using Newtons
1899 ! iteration method.
1900 !
1901 ! On entry:
1902 integer NLAT ! the number of latitudes (degree of the polynomial)
1903 !
1904 ! On exit: for each Gaussian latitude
1906 double precision, dimension(NLAT) :: LATG ! Latitude
1908 ! Approximations to a regular latitude grid:
1909 real :: deltalat
1910 real :: startlat
1912 !-----------------------------------------------------------------------
1914 integer :: iskip = 15
1915 double precision :: sum1 = 0.
1916 double precision :: sum2 = 0.
1917 double precision :: sum3 = 0.
1918 double precision :: sum4 = 0.
1919 double precision :: xn
1921 integer, save :: SAVE_NLAT = -99
1922 real, save :: save_deltalat = -99.
1923 real, save :: save_startlat = -99.
1925 double precision, dimension(nlat) :: COSC, SINC
1926 double precision, parameter :: PI = 3.141592653589793
1927 !
1928 ! -convergence criterion for iteration of cos latitude
1929 double precision, parameter :: XLIM = 1.0E-14
1930 integer :: nzero, i, j
1931 double precision :: fi, fi1, a, b, g, gm, gp, gt, delta, c, d
1933 if (nlat == save_nlat) then
1934 deltalat = save_deltalat
1935 startlat = save_startlat
1936 return
1937 endif
1938 !
1939 ! -the number of zeros between pole and equator
1940 NZERO = NLAT/2
1941 !
1942 ! -set first guess for cos(colat)
1943 DO I=1,NZERO
1944 COSC(I) = SIN( (I-0.5)*PI/NLAT + PI*0.5 )
1945 ENDDO
1946 !
1947 ! -constants for determining the derivative of the polynomial
1948 FI = NLAT
1949 FI1 = FI+1.0
1950 A = FI*FI1 / SQRT(4.0*FI1*FI1-1.0)
1951 B = FI1*FI / SQRT(4.0*FI*FI-1.0)
1952 !
1953 ! -loop over latitudes, iterating the search for each root
1954 DO I=1,NZERO
1955 J=0
1956 !
1957 ! -determine the value of the ordinary Legendre polynomial for
1958 ! -the current guess root
1959 LOOP30 : DO
1960 CALL LGORD( G, COSC(I), NLAT )
1961 !
1962 ! -determine the derivative of the polynomial at this point
1963 CALL LGORD( GM, COSC(I), NLAT-1 )
1964 CALL LGORD( GP, COSC(I), NLAT+1 )
1965 GT = (COSC(I)*COSC(I)-1.0) / (A*GP-B*GM)
1966 !
1967 ! -update the estimate of the root
1968 DELTA = G*GT
1969 COSC(I) = COSC(I) - DELTA
1970 !
1971 ! -if convergence criterion has not been met, keep trying
1972 J = J+1
1973 IF( ABS(DELTA).LE.XLIM ) EXIT LOOP30
1974 ENDDO LOOP30
1975 ENDDO
1976 !
1977 ! Determine the sin(colat)
1978 SINC(1:NZERO) = SIN(ACOS(COSC(1:NZERO)))
1979 !
1980 ! -if NLAT is odd, set values at the equator
1981 IF( MOD(NLAT,2) .NE. 0 ) THEN
1982 I = NZERO+1
1983 SINC(I) = 1.0
1984 latg(i) = 0.
1985 END IF
1987 ! Set the latitudes.
1989 latg(1:NZERO) = dacos(sinc(1:NZERO)) * 180. / pi
1991 ! Determine the southern hemisphere values by symmetry
1992 do i = 1, nzero
1993 latg(nlat-nzero+i) = -latg(nzero+1-i)
1994 enddo
1997 ! Now that we have the true values, find some approximate values.
1999 xn = float(nlat-iskip*2)
2000 do i = iskip+1, nlat-iskip
2001 sum1 = sum1 + latg(i)*float(i)
2002 sum2 = sum2 + float(i)
2003 sum3 = sum3 + latg(i)
2004 sum4 = sum4 + float(i)**2
2005 enddo
2007 b = (xn*sum1 - sum2*sum3) / (xn*sum4 - sum2**2)
2008 a = (sum3 - b * sum2) / xn
2010 deltalat = sngl(b)
2011 startlat = sngl(a + b)
2013 save_nlat = nlat
2014 save_deltalat = deltalat
2015 save_startlat = startlat
2017 contains
2018 SUBROUTINE LGORD( F, COSC, N )
2019 implicit none
2020 !
2021 ! LGORD calculates the value of an ordinary Legendre polynomial at a
2022 ! latitude.
2023 !
2024 ! On entry:
2025 ! COSC - cos(colatitude)
2026 ! N - the degree of the polynomial
2027 !
2028 ! On exit:
2029 ! F - the value of the Legendre polynomial of degree N at
2030 ! latitude asin(COSC)
2031 double precision :: s1, c4, a, b, fk, f, cosc, colat, c1, fn, ang
2032 integer :: n, k
2034 !------------------------------------------------------------------------
2036 colat = acos(cosc)
2037 c1 = sqrt(2.0)
2038 do k=1,n
2039 c1 = c1 * sqrt( 1.0 - 1.0/(4*k*k) )
2040 enddo
2041 fn = n
2042 ang= fn * colat
2043 s1 = 0.0
2044 c4 = 1.0
2045 a =-1.0
2046 b = 0.0
2047 do k=0,n,2
2048 if (k.eq.n) c4 = 0.5 * c4
2049 s1 = s1 + c4 * cos(ang)
2050 a = a + 2.0
2051 b = b + 1.0
2052 fk = k
2053 ang= colat * (fn-fk-2.0)
2054 c4 = ( a * (fn-b+1.0) / ( b * (fn+fn-a) ) ) * c4
2055 enddo
2056 f = s1 * c1
2057 end subroutine lgord
2059 END SUBROUTINE GRIBLGG
2061 SUBROUTINE REORDER_IT (a, nx, ny, dx, dy, iorder)
2063 use module_debug
2065 implicit none
2066 integer :: nx, ny, iorder
2067 integer :: i, j, k, m
2068 real :: dx, dy
2069 real, dimension(nx*ny) :: a, z
2071 if (iorder .eq. 0 .and. dx .gt. 0. .and. dy .lt. 0) return
2072 k = 0
2073 call mprintf(.true.,DEBUG, &
2074 "Reordering GRIB array : dx = %f , dy = %f , iorder = %i", &
2075 f1=dx,f2=dy,i1=iorder)
2076 if (iorder .eq. 0 ) then
2077 if ( dx .lt. 0 .and. dy .lt. 0. ) then
2078 do j = 1, ny
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 .lt. 0 .and. dy .gt. 0. ) then
2086 do j = ny, 1, -1
2087 do i = nx, 1, -1
2088 k = k + 1
2089 m = i * j
2090 z(k) = a(m)
2091 enddo
2092 enddo
2093 else if ( dx .gt. 0 .and. dy .gt. 0. ) then
2094 do j = ny, 1, -1
2095 do i = 1, nx
2096 k = k + 1
2097 m = i * j
2098 z(k) = a(m)
2099 enddo
2100 enddo
2101 endif
2102 else
2103 if ( dx .gt. 0 .and. dy .lt. 0. ) then
2104 do i = 1, nx
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 = 1, ny
2114 k = k + 1
2115 m = i * j
2116 z(k) = a(m)
2117 enddo
2118 enddo
2119 else if ( dx .lt. 0 .and. dy .lt. 0. ) then
2120 do i = nx, 1, -1
2121 do j = ny, 1, -1
2122 k = k + 1
2123 m = i * j
2124 z(k) = a(m)
2125 enddo
2126 enddo
2127 else if ( dx .gt. 0 .and. dy .gt. 0. ) then
2128 do i = 1, nx
2129 do j = ny, 1, -1
2130 k = k + 1
2131 m = i * j
2132 z(k) = a(m)
2133 enddo
2134 enddo
2135 endif
2136 endif
2137 ! now put it back in the 1-d array and reset the dx and dy
2138 do k = 1, nx*ny
2139 a(k) = z(k)
2140 enddo
2141 dx = abs ( dx)
2142 dy = -1 * abs(dy)
2143 return
2144 END SUBROUTINE REORDER_IT