source: trunk/mesoscale/LMD_MM_MARS/SRC/WPS/ungrib/src/gribcode.F90 @ 87

Last change on this file since 87 was 11, checked in by aslmd, 14 years ago

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