ungrib/src/rd_grib1.F

   1 !*****************************************************************************!
   2 ! Subroutine RD_GRIB1                                                         !
   3 !                                                                             !
   4 ! Purpose:                                                                    !
   5 !    Read one record from the input GRIB file.  Based on the information in   !
   6 !    the GRIB header and the user-defined Vtable, decide whether the field in !
   7 !    the GRIB record is one to process or to skip.  If the field is one we    !
   8 !    want to keep, extract the data from the GRIB record, and pass the data   !
   9 !    back to the calling routine.                                             !
  10 !                                                                             !
  11 ! Argument list:                                                              !
  12 !    Input:                                                                   !
  13 !       IUNIT   : "Unit Number" to open and read from.  Not really a Fortran  !
  14 !                 unit number, since we do not do Fortran I/O for the GRIB     !
  15 !                 files.  Nor is it a UNIX File Descriptor returned from a C  !
  16 !                 OPEN statement.  It is really just an array index to the     !
  17 !                 array (IUARR) where the UNIX File Descriptor values are     !
  18 !                 stored.                                                     !
  19 !       GRIBFLNM: File name to open, if it is not already open.               !
  20 !       IUARR   : Array to hold UNIX File descriptors retured from a C open   !
  21 !                 statement.  If the value of IUARR(IUNIT) is zero, then the  !
  22 !                 file GRIBFLNM must be opened, and the value of IUARR(IUNIT) !
  23 !                 becomes the UNIX File descriptor to read from.              !
  24 !       DEBUG_LEVEL Integer for various amounts of printout.                    !
  25 !                                                                             !
  26 !    Output:                                                                  !
  27 !       LEVEL    : The pressure-level (Pa) of the field to process.           !
  28 !       FIELD    : The field name of the field to process.  NULL is returned  !
  29 !                  if we do not want to process the field we read.            !
  30 !       HDATE    : The 19-character date of the field to process.             !
  31 !       IERR     : Error flag: 0 - no error on read from GRIB file.           !
  32 !                              1 - Hit the end of the GRIB file.              !
  33 !                              2 - The file GRIBFLNM we tried to open does    !
  34 !                                  not exist.                                 !
  35 ! Externals                                                                   !
  36 !     Module TABLE                                                            !
  37 !     Module GRIDINFO                                                         !
  38 !     Subroutine C_OPEN                                                       !
  39 !     Subroutine DEALLOGRIB                                                   !
  40 !     Subroutine GRIBGET                                                      !
  41 !     Subroutine GRIBHEADER                                                   !
  42 !     Subroutine GET_SEC1                                                     !
  43 !     Subroutine GET_SEC2                                                     !
  44 !     Subroutine GET_GRIDINFO                                                 !
  45 !     Subroutine BUILD_HDATE                                                  !
  46 !     Subroutine GETH_NEWDATE                                                 !
  47 !     Subroutine GRIBDATA                                                     !
  48 !                                                                             !
  49 ! Side Effects                                                                !
  50 !     File GRIBFLNM is opened, as necessary                                   !
  51 !                                                                             !
  52 !     Variable MAP from module GRIDINFO is filled in.                         !
  53 !                                                                             !
  54 !     Numerous side effects from the GRIB-processing routines.                !
  55 !                                                                             !
  56 ! Kevin W. Manning                                                            !
  57 ! NCAR/MMM                                                                    !
  58 ! Summer, 1998, and continuing                                                !
  59 ! SDG                                                                         !
  60 !                                                                             !
  61 !*****************************************************************************!
  62 SUBROUTINE rd_grib1(IUNIT, gribflnm, level, field, hdate,  &
  63      ierr, iuarr, debug_level, ec_rec_len)
  64   use table
  65   use gridinfo
  66   use datarray
  67   use module_debug
  68
  69   implicit none
  70
  71   integer :: debug_level, ec_rec_len
  72   integer :: iunit ! Array number in IUARR assigned to the C read pointer.
  73   integer, dimension(100) :: KSEC1
  74   integer, dimension(10) :: KSEC2
  75   integer, dimension(40) :: infogrid
  76   real, dimension(40) :: ginfo
  77 !
  78 !-----------------------------------------------------------------------
  79   integer :: iparm, ktype
  80   logical :: lopen
  81
  82   integer :: icenter, iprocess, iscan, ii, isb
  83   integer year, month, day, hour, minute, second, icc, iyy
  84   integer :: fcst
  85   real :: level
  86   character(LEN=*) :: field
  87   character(LEN=132) :: gribflnm
  88   character(LEN=8) :: tmp8
  89   integer, dimension(255) :: iuarr
  90   integer :: ierr, iostat, nunit
  91   integer :: i, lvl2, lvl1
  92   character(LEN=19) :: hdate
  93   integer :: igherr
  94
  95 ! Variables for thinned grids:
  96   logical :: lthinned = .FALSE.
  97   real, allocatable, dimension(:) :: thinnedDataArray
  98   integer, dimension(74) :: npoints_acc
  99   real :: mj, xmj
 100   integer :: np, ny, nx
 101   real :: Va, Vb, Vc, Vd
 102   real, external :: oned
 103
 104   ierr = 0
 105
 106 ! If the file GRIBFLNM has not been opened, then IUARR(IUNIT) should be Zero.
 107 ! In this case, open the file GRIBFLNM, and store the UNIX File descriptor
 108 ! in to IUARR(IUNIT).  This way, we will know what UNIX File descriptor to use
 109 ! next time we call this RD_GRIB subroutine.
 110 !
 111   if (iuarr(iunit).eq.0) then
 112      if (debug_level.gt.0) then
 113         call c_open(iunit, nunit, gribflnm, 1, ierr,  1)
 114      else
 115         call c_open(iunit, nunit, gribflnm, 1, ierr, -1)
 116      endif
 117      if (ierr.ne.0) then
 118         call deallogrib
 119         ierr = 2
 120         return
 121      endif
 122      iuarr(iunit) = nunit
 123   endif
 124
 125 ! Read a single GRIB record, but do no unpacking now:
 126
 127   call gribget(iuarr(iunit), ierr, ec_rec_len)
 128
 129   if (ierr.ne.0) then
 130      call mprintf(.true.,DEBUG,"RD_GRIB1 gribget read error, ierr = %i",i1=ierr)
 131      call deallogrib
 132      return
 133   endif
 134 !
 135 ! Unpack the header information:
 136 !
 137   call gribheader(debug_level,igherr,ec_rec_len)
 138   if (igherr /= 0) then
 139      field = "NULL"
 140      call deallogrib
 141      return
 142   endif
 143 !
 144 ! Copy header information to arrays KSEC1, KSEC2, INFOGRID, and GRIDINFO
 145 !
 146   call get_sec1(ksec1)
 147   call get_sec2(ksec2)
 148   call get_gridinfo(infogrid, ginfo)
 149
 150   icenter = KSEC1(3)        ! Indicator of the source (center) of the data.
 151   iprocess = KSEC1(4)       ! Indicator of model (or whatever) which generated the data.
 152
 153   if (icenter.eq.7) then
 154     if (iprocess.eq.83 .or. iprocess.eq.84) then
 155       map%source = 'NCEP MESO NAM Model'
 156     elseif (iprocess.eq.81) then
 157       map%source = 'NCEP GFS Analysis'
 158     elseif (iprocess.eq.82) then
 159       map%source = 'NCEP GFS GDAS/FNL'
 160     elseif (iprocess.eq.89) then
 161       map%source = 'NCEP NMM '
 162     elseif (iprocess.eq.96) then
 163       map%source = 'NCEP GFS Model'
 164     elseif (iprocess.eq.107) then
 165       map%source = 'NCEP GEFS'
 166     elseif (iprocess.eq.109) then
 167       map%source = 'NCEP RTMA'
 168     elseif (iprocess.eq.86 .or. iprocess.eq.100) then
 169       map%source = 'NCEP RUC Model'    ! 60 km
 170     elseif (iprocess.eq.101) then
 171       map%source = 'NCEP RUC Model'    ! 40 km
 172     elseif (iprocess.eq.105) then
 173       map%source = 'NCEP RUC Model'    ! 20 km
 174     elseif (iprocess.eq.140) then
 175       map%source = 'NCEP NARR'
 176     elseif (iprocess.eq.195) then
 177       map%source = 'NCEP CDAS2'
 178     elseif (iprocess.eq.44) then
 179       map%source = 'NCEP SST Analysis'
 180     elseif (iprocess.eq.70) then
 181       map%source = 'GFDL Hurricane Model'
 182     elseif (iprocess.eq.129) then
 183       map%source = 'NCEP GODAS'
 184     elseif (iprocess.eq.25) then
 185       map%source = 'NCEP SNOW COVER ANALYSIS'
 186     else
 187       map%source = 'unknown model from NCEP'
 188     end if
 189 !  grid numbers only set for NCEP and AFWA models
 190     write(tmp8,'("GRID ",i3)') KSEC1(5)
 191     map%source(25:32) = tmp8
 192   else if (icenter .eq. 57) then
 193     if (iprocess .eq. 87) then
 194       map%source = 'AFWA AGRMET'
 195     else
 196       map%source = 'AFWA'
 197     endif
 198     write(tmp8,'("GRID ",i3)') KSEC1(5)
 199     map%source(25:32) = tmp8
 200   else if (icenter .eq. 58) then
 201     map%source = 'US Navy FNOC'
 202   else if (icenter .eq. 59) then
 203       if (iprocess .eq. 125) then
 204         map%source = 'NOAA GSD Rapid Refresh'
 205       else if (iprocess .eq. 105) then
 206         map%source = 'NOAA GSD'
 207       else
 208         print *,'Unknown GSD source'
 209         stop
 210       endif
 211   else if (icenter .eq. 60) then
 212         map%source = 'NCAR'
 213   else if (icenter .eq. 98) then
 214       map%source = 'ECMWF'
 215   else if (icenter .eq. 74 .or. icenter .eq. 75 ) then
 216       map%source = 'UKMO'
 217   else
 218     map%source = 'unknown model and orig center'
 219   end if
 220
 221   IPARM=KSEC1(7)            ! Indicator of parameter
 222   KTYPE=KSEC1(8)            ! Indicator of type of level
 223
 224 !   print *,' IPARM, KTYPE, KSEC1(9)', iparm,ktype,ksec1(9)
 225
 226   IF(KTYPE.EQ.1) THEN
 227      LVL1=KSEC1(9)
 228      LVL2=-99
 229   ELSEIF(KTYPE.EQ.100) THEN
 230      LVL1=FLOAT(KSEC1(9))  * 100.
 231      LVL2=-99
 232   ELSEIF(KTYPE.EQ.101) THEN
 233      LVL1=KSEC1(9)
 234      LVL2=KSEC1(10)
 235   ELSEIF(KTYPE.EQ.102) THEN
 236      LVL1=KSEC1(9)
 237      LVL2=-99
 238   ELSEIF(KTYPE.EQ.103) THEN
 239      LVL1=KSEC1(9)
 240      LVL2=-99
 241   ELSEIF(KTYPE.EQ.105) THEN
 242      LVL1=KSEC1(9)
 243      LVL2=-99
 244   ELSEIF(KTYPE.EQ.107) THEN
 245      LVL1=KSEC1(9)
 246      LVL2=-99
 247   ELSEIF(KTYPE.EQ.109) THEN
 248      LVL1=KSEC1(9)
 249      LVL2=-99
 250   ELSEIF(KTYPE.EQ.111) THEN
 251      LVL1=KSEC1(9)
 252      LVL2=-99
 253   ELSEIF(KTYPE.EQ.112) THEN ! Layer between two depths below surface
 254      LVL1=KSEC1(9)
 255      LVL2=KSEC1(10)
 256   ELSEIF(KTYPE.EQ.113) THEN
 257      LVL1=KSEC1(9)
 258      LVL2=-99
 259   ELSEIF(KTYPE.EQ.115) THEN
 260      LVL1=KSEC1(9)
 261      LVL2=-99
 262   ELSEIF(KTYPE.EQ.117) THEN
 263      LVL1=KSEC1(9)
 264      LVL2=-99
 265   ELSEIF(KTYPE.EQ.119) THEN
 266      LVL1=KSEC1(9)
 267      LVL2=-99
 268   ELSEIF(KTYPE.EQ.125) THEN
 269      LVL1=KSEC1(9)
 270      LVL2=-99
 271   ELSEIF(KTYPE.EQ.160) THEN
 272      LVL1=KSEC1(9)
 273      LVL2=-99
 274   ELSEIF(KTYPE.EQ.200) THEN
 275      LVL1=0
 276      LVL2=-99
 277   ELSEIF(KTYPE.EQ.201) THEN
 278      LVL1=0
 279      LVL2=-99
 280   ELSE
 281      LVL1=KSEC1(9)
 282      LVL2=KSEC1(10)
 283   ENDIF
 284
 285 ! Check to see that the combination of iparm, ktype, lvl1, and lvl2
 286 ! match what has been requested in the Vtable.  If not, set the field
 287 ! name to NULL, meaning that we do not want to process this one.
 288
 289   field = 'NULL'
 290   do i = 1, maxvar
 291      if (gcode(i).eq.iparm) then
 292         if (lcode(i).eq.ktype) then
 293            if ((level1(i).eq.lvl1) .or. (level1(i) == splatcode) ) then
 294               if (level2(i).eq.lvl2) then
 295                  field=namvar(i)
 296                  level = -999.
 297                  if (ktype.eq.100) then ! Pressure-level
 298                     level=lvl1
 299                  elseif (ktype.eq.102) then
 300                     level=201300.
 301                  elseif ((ktype.eq.116.and.lvl1.le.50.and.lvl2.eq.0) .or. &
 302                       (ktype.eq.105).or.(ktype.eq.1) .or. &
 303                       (ktype.eq.111).or.(ktype.eq.112) ) then
 304                     ! level=200100.
 305                     level = float(200000+iprty(i))
 306                  elseif (ktype.eq.109 .or. ktype.eq.107) then   ! hybrid or sigma levels
 307                     level = lvl1
 308                  elseif (ktype.eq. 6 ) then    ! max wind
 309                     level = 6.
 310                  elseif (ktype.eq. 7 ) then    ! trop
 311                     level = 7.
 312                  elseif (ktype .eq. 160 ) then    ! depth below sea-surface (m)
 313                     level = 201500.
 314                  elseif (ktype .eq. 237 .or. ktype .eq. 238 ) then  ! depth of ocean layer
 315                     level = 201600.
 316                  elseif (ktype .eq. 200 ) then  !column variable (TCDC,PWAT,etc.)
 317                     level = lvl1                !
 318                  endif
 319                  if (level .lt. -998. ) then
 320                    write(6,*) 'Could not find a level for this Vtable entry'
 321                    write(6,*) 'iparm = ',iparm,' ktype = ',ktype,' lvl1 = ',lvl1,' lvl2 = ',lvl2
 322                    write(6,*) 'Fix the Vtable or modify rd_grib1.F'
 323                    stop 'rd_grib1'
 324                  endif
 325               endif
 326            endif
 327         endif
 328      endif
 329   enddo
 330
 331   if (field .eq. 'NULL') then
 332      call deallogrib
 333      return
 334   endif
 335
 336   if ((field.eq.'WEASD').or.(field.eq.'SNOW')) then
 337      level = level + ksec1(19)+1
 338   endif
 339
 340 ! Build the 19-character date string, based on GRIB header date and time
 341 ! information, including forecast time information:
 342
 343   ICC=KSEC1(22)             ! CENTURY OF THE DATA
 344   IYY=KSEC1(11)             ! (TWO-DIGIT) YEAR OF THE DATA
 345   MONTH=KSEC1(12)           ! MONTH OF THE DATA
 346   DAY=KSEC1(13)             ! DAY OF THE DATA
 347   HOUR=KSEC1(14)            ! HOUR OF THE DATA
 348   MINUTE=KSEC1(15)          ! MINUTE OF THE DATA
 349   SECOND=0
 350   if (ksec1(19) == 3) then
 351      FCST = (KSEC1(17) + KSEC1(18))/2
 352 !  TEMPORARY AFWA FIX
 353 !  elseif (ksec1(19) == 4 .or. ksec1(19) == 5) then
 354    elseif (ksec1(19) == 4 .or. ksec1(19) == 5 .or. ksec1(19) == 7) then
 355      FCST = KSEC1(18)
 356   else
 357      FCST = KSEC1(17)
 358   endif
 359 ! convert the fcst units to hours if necessary
 360   if (ksec1(16) .eq. 254 ) then
 361     fcst = fcst/3600.
 362   elseif (ksec1(16) .eq. 0 ) then
 363     fcst = fcst/60.
 364   endif
 365
 366   if (IYY.EQ.00) then
 367      YEAR = ICC*100
 368   else
 369      YEAR = (ICC-1)*100 + IYY
 370   endif
 371
 372   hdate(1:19) = '                   '
 373   call build_hdate(hdate,year,month,day,hour,minute,second)
 374
 375   call geth_newdate(hdate,hdate,3600*fcst)
 376
 377 ! Store information about the grid on which the data is.
 378 ! This stuff gets stored in the MAP variable, as defined in module GRIDINFO
 379
 380   map%startloc = 'SWCORNER'
 381   map%grid_wind = .true.
 382 ! NCEP's grib1 messages (in GDS Octet 17, the Resolution and Component Flags)
 383 ! all have '0' for the earth radius flag which the documentation (written by NCEP)
 384 ! says is 6367.47, but they really use 6371.229. Hardcode it.
 385 ! It's not clear what ECMWF uses. One place says 6367.47 and another 6371.229.
 386   if ( index(map%source,'NCEP') .ne. 0 ) then
 387     map%r_earth = 6371.229
 388   else
 389     map%r_earth = 6367.47
 390   endif
 391
 392   if (ksec2(4).eq.0) then ! Lat/Lon grid
 393      map%igrid = 0
 394      map%nx = infogrid(1)
 395      map%ny = infogrid(2)
 396      map%dx = ginfo(8)
 397      map%dy = ginfo(9)
 398      map%lat1 = ginfo(3)
 399      map%lon1 = ginfo(4)
 400      !  If this is global data, then the dx and dy are more accurately
 401      !  computed by the number of points than the 3 digits grib permits.
 402      if ( ABS(map%nx * map%dx - 360.) .lt. 1 ) then
 403         if ( ABS ( map%dx - (360./real(map%nx)) ) .gt. 0.00001 ) then
 404            !print *,'old dx = ',ginfo(8)
 405            map%dx = 360./real(map%nx)
 406            !print *,'new dx = ',map%dx
 407         endif
 408      endif
 409      if ( ABS((map%ny-1) * map%dy - 2.*abs(map%lat1)) .lt. 1. ) then
 410         if ( ABS ( map%dy - (2.*abs(map%lat1)/real(map%ny-1)) ) .gt. 0.00001 ) then
 411            !print *,'old dy = ',ginfo(9)
 412            map%dy = 2.*abs(map%lat1)/real(map%ny-1)
 413            !print *,'new dy = ',map%dy
 414         endif
 415      endif
 416      write(tmp8,'(b8.8)') infogrid(5)
 417      if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
 418      if (icenter .eq. 7 .and. KSEC1(5) .eq. 173 ) then  ! correction for ncep grid 173
 419        map%lat1 = 89.958333
 420        map%lon1 = 0.041667
 421        map%dx = 0.083333333 * sign(1.0,map%dx)
 422        map%dy = 0.083333333 * sign(1.0,map%dy)
 423      endif
 424 ! correction for ncep grid 229   added 5/3/07   JFB
 425      if (icenter .eq. 7 .and. KSEC1(5) .eq. 229 ) then
 426        if (ginfo(3) .gt. 89. .and. ginfo(9) .gt. 0.) then
 427          map%dy = -1. * map%dy
 428        endif
 429      endif
 430
 431 !    print *, "CE map stuff", map%igrid, map%nx, map%ny, map%dx, &
 432 !    map%dy, map%lat1, map%lon1
 433
 434   elseif (ksec2(4).eq.1) then ! Mercator Grid
 435      map%igrid = 1
 436      map%nx = infogrid(1)
 437      map%ny = infogrid(2)
 438      map%dx = ginfo(8)  ! km
 439      map%dy = ginfo(9)
 440      map%truelat1 = ginfo(5)
 441      map%truelat2 = 0.
 442      map%lov = 0.
 443      map%lat1 = ginfo(3)
 444      map%lon1 = ginfo(4)
 445      write(tmp8,'(b8.8)') infogrid(5)
 446      if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
 447
 448   elseif (ksec2(4).eq.3) then ! Lambert Conformal Grid
 449      map%igrid = 3
 450      map%nx = infogrid(1)
 451      map%ny = infogrid(2)
 452      map%lov = ginfo(6)
 453      map%truelat1 = ginfo(11)
 454      map%truelat2 = ginfo(12)
 455      map%dx = ginfo(7)
 456      map%dy = ginfo(8)
 457      map%lat1 = ginfo(3)
 458      map%lon1 = ginfo(4)
 459      write(tmp8,'(b8.8)') infogrid(5)
 460      if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
 461 !    if (tmp8(2:2) .eq. '0') map%r_earth = 6367.47
 462
 463   elseif(ksec2(4).eq.4) then ! Gaussian Grid
 464      map%igrid = 4
 465      map%nx = infogrid(1)
 466      map%ny = infogrid(2)
 467      map%dx = ginfo(8)
 468 !    map%dy = ginfo(19)
 469      map%dy = real (infogrid(9))
 470      map%lon1 = ginfo(4)
 471      map%lat1 = ginfo(3)
 472      write(tmp8,'(b8.8)') infogrid(5)
 473      if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
 474 !  If this is global data, then the dx and dy are more accurately
 475 !  computed by the number of points than the 3 digits grib permits.
 476      if ( ABS(map%nx * map%dx - 360.) .lt. 1. ) then
 477         if ( ABS ( map%dx - (360./real(map%nx)) ) .gt. 0.00001 ) then
 478          ! print *,'old dx = ',ginfo(8)
 479            map%dx = 360./real(map%nx)
 480          ! print *,'new dx = ',map%dx
 481         endif
 482      endif
 483
 484
 485   elseif (ksec2(4).eq.5) then ! Polar-Stereographic Grid.
 486      map%igrid = 5
 487      map%nx = infogrid(1)
 488      map%ny = infogrid(2)
 489      map%lov = ginfo(6)
 490      map%truelat1 = 60.
 491      map%truelat2 = 91.
 492      map%dx = ginfo(7)
 493      map%dy = ginfo(8)
 494      map%lat1 = ginfo(3)
 495      map%lon1 = ginfo(4)
 496      write(tmp8,'(b8.8)') infogrid(5)
 497      if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
 498
 499   else
 500      print*, 'Unknown Data Representation Type, ksec2(4)= ', ksec2(4)
 501      stop 'rd_grib1'
 502   endif
 503
 504 111  format(' igrid      : ', i3, /, &
 505           ' nx, ny     : ', 2I4, /, &
 506           ' truelat1, 2: ', 2F10.4, /, &
 507           ' Center Lon : ', F10.4, /, &
 508           ' LatLon(1,1): ', 2F10.4, /, &
 509           ' DX, DY     : ', F10.4, F10.4)
 510
 511 ! Special for NCEP/NCAR Reanalysis Project:
 512 !      Throw out PSFC on lat/lon grid (save gaussian version)
 513   if ((icenter.eq.7).and.(iprocess.eq.80)) then   ! Careful! This combination may refer
 514                                                       ! to other products as well.
 515      if ((field.eq.'PSFC').and.(ksec2(4).eq.0)) then
 516         field='NULL'
 517         call deallogrib
 518         return
 519      endif
 520   endif
 521
 522   if (allocated(rdatarray)) deallocate(rdatarray)
 523   allocate(rdatarray(map%nx * map%ny))
 524
 525 ! If nx=65535, assume the grid is a thinned grid.
 526 ! Process only the NCEP grid IDs is 37 to 44.
 527   if (map%nx.eq.65535) then
 528      if ( (icenter .ne. 7) .or. (KSEC1(5).lt.37) .or. (KSEC1(5).gt.44) ) then
 529         write(*,*) 'Originating center is ',icenter
 530         write(*,*) 'Grid ID is ',KSEC1(5),' Only WAFS grids 37-44 are supported'
 531         write(*,'(" ***** STOP in Subroutine RD_GRIB1.",//)')
 532         stop
 533      endif
 534      lthinned = .TRUE.
 535      map%nx = 73
 536      map%dx = 1.25
 537   else
 538      lthinned = .FALSE.
 539   endif
 540
 541 ! Unpack the 2D slab from the GRIB record, and put it in array rdatarray
 542
 543   if (lthinned) then
 544     if (allocated(thinnedDataArray)) deallocate(thinnedDataArray)
 545     allocate(thinnedDataArray(map%nx * map%ny))
 546     call gribdata(thinnedDataArray,3447)
 547
 548     ! Calculate how many points for each latitude, and accumulate into array
 549     if ((KSEC1(5).ge.37).and.(KSEC1(5).le.40)) then
 550        ! Northern hemisphere:
 551        npoints_acc(1)=0
 552        npoints_acc(2)=73
 553        do i=1,72
 554           np = int(2.0+(90.0/1.25)*cos(i*1.25*3.1415926/180.0))
 555           npoints_acc(i+2)=npoints_acc(i+1)+np
 556        enddo
 557     else
 558        ! Southern Hemisphere:
 559        npoints_acc(1)=0
 560        npoints_acc(2)=2
 561        do i=1,71
 562           ii = 72-i
 563           np = int(2.0+(90.0/1.25)*cos(ii*1.25*3.1415926/180.0))
 564           npoints_acc(i+2)=npoints_acc(i+1)+np
 565        enddo
 566        npoints_acc(74) = npoints_acc(73) + 73
 567     endif
 568
 569     ! for row number i (where i=1 is the southern edge of the grid)
 570     !   npoints_acc(i+1)-npoints_acc(i) = number of points in this line
 571     !   npoints_acc(i)+1 = index into thinned array for first point of line
 572
 573     do ny=1,73
 574        np = npoints_acc(ny+1)-npoints_acc(ny) ! Number of points in this line.
 575        do nx=1,73
 576           ! Calulate the x index (mj) of thinned array (real value)
 577           mj = (nx-1.0)*(np-1.0)/(72.0)
 578
 579           if (abs(mj - int(mj)) < 1.E-10) then
 580              rdatarray((ny-1)*73+nx) = thinnedDataArray(npoints_acc(ny)+1+int(mj))
 581           else
 582              ! Get the 2 closest values from thinned array
 583              Vb = thinnedDataArray(npoints_acc(ny)+1+int(mj))
 584              Vc = thinnedDataArray(npoints_acc(ny)+1+int(mj)+1)
 585              ! Get the next two closest, if available:
 586              Va = -999999.
 587              Vd = -999999.
 588              if (mj > 1.0) then
 589                 Va = thinnedDataArray(npoints_acc(ny)+1+int(mj)-1)
 590              endif
 591              if (mj < np-2) then
 592                 Vd = thinnedDataArray(npoints_acc(ny)+1+int(mj)+2)
 593              endif
 594
 595              if ((Va < -999998.) .or. (Vd < -999998.)) then
 596                 ! Use 2-point linear interpolation.
 597                 rdatarray((ny-1)*73+nx) = Vb*(int(mj)+1.0-mj) + Vc*(mj-int(mj))
 598              else
 599                 ! Use 4-point overlapping parabolic interpolation.
 600                 xmj = mj - float(int(mj))
 601                 rdatarray((ny-1)*73+nx) = oned(xmj,Va,Vb,Vc,Vd)
 602              endif
 603           endif
 604        enddo
 605     enddo
 606 else
 607   call gribdata(rdatarray,map%nx*map%ny)
 608 endif
 609
 610 ! Some grids are broken and need to be reordered (e.g. NCEP-II in 1997).
 611 ! WPS assumes that the grids are ordered consistently with the start location.
 612
 613   call mprintf(.true.,DEBUG, &
 614   "RD_GRIB1 icenter = %i , iprocess = %i , grid = %i",i1=icenter,i2=iprocess,i3=KSEC1(5))
 615   if (icenter .eq. 7 .and. iprocess .eq. 0 .and. KSEC1(5) .eq. 2 ) then
 616   call mprintf(.true.,DEBUG, &
 617   "resetting NCEP2 dx and dy. If this is not NCEP2 data you must modify rd_grib1.f90")
 618   call mprintf(.true.,DEBUG, &
 619   "field = %s , dx = %f , dy = %f , i10 = %i",s1=field,f1=map%dx,f2=map%dy,i1=infogrid(10))
 620      map%dx = 2.5
 621      map%dy = -2.5
 622 !   call reorder_it (rdatarray, map%nx, map%ny, map%dx, map%dy, infogrid(10))
 623   endif
 624
 625 ! Deallocate a couple of arrays that may have been allocated by the
 626 ! GRIB decoding routines.
 627
 628   call deallogrib
 629
 630 END subroutine rd_grib1
 631
 632 real function oned(x, a, b, c, d) Result (Answer)
 633   implicit none
 634   real :: x ! Proportion of the way between B and C.  Between 0.0 and 1.0
 635   real :: a, b, c, d
 636
 637   if (abs(x) < 1.E-10) then
 638      Answer = B
 639      return
 640   endif
 641   IF(abs(x-1.) < 1.E-10) then
 642      Answer = C
 643      return
 644   endif
 645   Answer = (1.0-X)*(B+X*(0.5*(C-A)+X*(0.5*(C+A)-B)))+X*(C+(1.0-X)*(0.5 &
 646        *(B-D)+(1.0-X)*(0.5*(B+D)-C)))
 647 end function oned