| blob | d3a56b0fed55e6c46a4a860b42e8b804172a1b85 |
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)
64 use table
65 use gridinfo
66 use datarray
67 use module_debug
69 implicit none
71 integer :: debug_level
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
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
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
104 ierr = 0
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
125 ! Read a single GRIB record, but do no unpacking now:
127 call gribget(iuarr(iunit), ierr)
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)
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)
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.
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
221 IPARM=KSEC1(7) ! Indicator of parameter
222 KTYPE=KSEC1(8) ! Indicator of type of level
224 ! print *,' IPARM, KTYPE, KSEC1(9)', iparm,ktype,ksec1(9)
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
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.
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
331 if (field .eq. 'NULL') then
332 call deallogrib
333 return
334 endif
336 if ((field.eq.'WEASD').or.(field.eq.'SNOW')) then
337 level = level + ksec1(19)+1
338 endif
340 ! Build the 19-character date string, based on GRIB header date and time
341 ! information, including forecast time information:
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
366 if (IYY.EQ.00) then
367 YEAR = ICC*100
368 else
369 YEAR = (ICC-1)*100 + IYY
370 endif
372 hdate(1:19) = ' '
373 call build_hdate(hdate,year,month,day,hour,minute,second)
375 call geth_newdate(hdate,hdate,3600*fcst)
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
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
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
431 ! print *, "CE map stuff", map%igrid, map%nx, map%ny, map%dx, &
432 ! map%dy, map%lat1, map%lon1
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.
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
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.
475 elseif (ksec2(4).eq.5) then ! Polar-Stereographic Grid.
476 map%igrid = 5
477 map%nx = infogrid(1)
478 map%ny = infogrid(2)
479 map%lov = ginfo(6)
480 map%truelat1 = 60.
481 map%truelat2 = 91.
482 map%dx = ginfo(7)
483 map%dy = ginfo(8)
484 map%lat1 = ginfo(3)
485 map%lon1 = ginfo(4)
486 write(tmp8,'(b8.8)') infogrid(5)
487 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
489 else
490 print*, 'Unknown Data Representation Type, ksec2(4)= ', ksec2(4)
491 stop 'rd_grib1'
492 endif
494 111 format(' igrid : ', i3, /, &
495 ' nx, ny : ', 2I4, /, &
496 ' truelat1, 2: ', 2F10.4, /, &
497 ' Center Lon : ', F10.4, /, &
498 ' LatLon(1,1): ', 2F10.4, /, &
499 ' DX, DY : ', F10.4, F10.4)
501 ! Special for NCEP/NCAR Reanalysis Project:
502 ! Throw out PSFC on lat/lon grid (save gaussian version)
503 if ((icenter.eq.7).and.(iprocess.eq.80)) then ! Careful! This combination may refer
504 ! to other products as well.
505 if ((field.eq.'PSFC').and.(ksec2(4).eq.0)) then
506 field='NULL'
507 call deallogrib
508 return
509 endif
510 endif
512 if (allocated(rdatarray)) deallocate(rdatarray)
513 allocate(rdatarray(map%nx * map%ny))
515 ! If nx=65535, assume the grid is a thinned grid.
516 ! Process only the NCEP grid IDs is 37 to 44.
517 if (map%nx.eq.65535) then
518 if ( (icenter .ne. 7) .or. (KSEC1(5).lt.37) .or. (KSEC1(5).gt.44) ) then
519 write(*,*) 'Originating center is ',icenter
520 write(*,*) 'Grid ID is ',KSEC1(5),' Only WAFS grids 37-44 are supported'
521 write(*,'(" ***** STOP in Subroutine RD_GRIB1.",//)')
522 stop
523 endif
524 lthinned = .TRUE.
525 map%nx = 73
526 map%dx = 1.25
527 else
528 lthinned = .FALSE.
529 endif
531 ! Unpack the 2D slab from the GRIB record, and put it in array rdatarray
533 if (lthinned) then
534 if (allocated(thinnedDataArray)) deallocate(thinnedDataArray)
535 allocate(thinnedDataArray(map%nx * map%ny))
536 call gribdata(thinnedDataArray,3447)
538 ! Calculate how many points for each latitude, and accumulate into array
539 if ((KSEC1(5).ge.37).and.(KSEC1(5).le.40)) then
540 ! Northern hemisphere:
541 npoints_acc(1)=0
542 npoints_acc(2)=73
543 do i=1,72
544 np = int(2.0+(90.0/1.25)*cos(i*1.25*3.1415926/180.0))
545 npoints_acc(i+2)=npoints_acc(i+1)+np
546 enddo
547 else
548 ! Southern Hemisphere:
549 npoints_acc(1)=0
550 npoints_acc(2)=2
551 do i=1,71
552 ii = 72-i
553 np = int(2.0+(90.0/1.25)*cos(ii*1.25*3.1415926/180.0))
554 npoints_acc(i+2)=npoints_acc(i+1)+np
555 enddo
556 npoints_acc(74) = npoints_acc(73) + 73
557 endif
559 ! for row number i (where i=1 is the southern edge of the grid)
560 ! npoints_acc(i+1)-npoints_acc(i) = number of points in this line
561 ! npoints_acc(i)+1 = index into thinned array for first point of line
563 do ny=1,73
564 np = npoints_acc(ny+1)-npoints_acc(ny) ! Number of points in this line.
565 do nx=1,73
566 ! Calulate the x index (mj) of thinned array (real value)
567 mj = (nx-1.0)*(np-1.0)/(72.0)
569 if (abs(mj - int(mj)) < 1.E-10) then
570 rdatarray((ny-1)*73+nx) = thinnedDataArray(npoints_acc(ny)+1+int(mj))
571 else
572 ! Get the 2 closest values from thinned array
573 Vb = thinnedDataArray(npoints_acc(ny)+1+int(mj))
574 Vc = thinnedDataArray(npoints_acc(ny)+1+int(mj)+1)
575 ! Get the next two closest, if available:
576 Va = -999999.
577 Vd = -999999.
578 if (mj > 1.0) then
579 Va = thinnedDataArray(npoints_acc(ny)+1+int(mj)-1)
580 endif
581 if (mj < np-2) then
582 Vd = thinnedDataArray(npoints_acc(ny)+1+int(mj)+2)
583 endif
585 if ((Va < -999998.) .or. (Vd < -999998.)) then
586 ! Use 2-point linear interpolation.
587 rdatarray((ny-1)*73+nx) = Vb*(int(mj)+1.0-mj) + Vc*(mj-int(mj))
588 else
589 ! Use 4-point overlapping parabolic interpolation.
590 xmj = mj - float(int(mj))
591 rdatarray((ny-1)*73+nx) = oned(xmj,Va,Vb,Vc,Vd)
592 endif
593 endif
594 enddo
595 enddo
596 else
597 call gribdata(rdatarray,map%nx*map%ny)
598 endif
600 ! Some grids are broken and need to be reordered (e.g. NCEP-II in 1997).
601 ! WPS assumes that the grids are ordered consistently with the start location.
603 call mprintf(.true.,DEBUG, &
604 "RD_GRIB1 icenter = %i , iprocess = %i , grid = %i",i1=icenter,i2=iprocess,i3=KSEC1(5))
605 if (icenter .eq. 7 .and. iprocess .eq. 0 .and. KSEC1(5) .eq. 2 ) then
606 call mprintf(.true.,DEBUG, &
607 "resetting NCEP2 dx and dy. If this is not NCEP2 data you must modify rd_grib1.f90")
608 call mprintf(.true.,DEBUG, &
609 "field = %s , dx = %f , dy = %f , i10 = %i",s1=field,f1=map%dx,f2=map%dy,i1=infogrid(10))
610 map%dx = 2.5
611 map%dy = -2.5
612 ! call reorder_it (rdatarray, map%nx, map%ny, map%dx, map%dy, infogrid(10))
613 endif
615 ! Deallocate a couple of arrays that may have been allocated by the
616 ! GRIB decoding routines.
618 call deallogrib
620 END subroutine rd_grib1
622 real function oned(x, a, b, c, d) Result (Answer)
623 implicit none
624 real :: x ! Proportion of the way between B and C. Between 0.0 and 1.0
625 real :: a, b, c, d
627 if (abs(x) < 1.E-10) then
628 Answer = B
629 return
630 endif
631 IF(abs(x-1.) < 1.E-10) then
632 Answer = C
633 return
634 endif
635 Answer = (1.0-X)*(B+X*(0.5*(C-A)+X*(0.5*(C+A)-B)))+X*(C+(1.0-X)*(0.5 &
636 *(B-D)+(1.0-X)*(0.5*(B+D)-C)))
637 end function oned