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