| blob | dd27a9f74a0b9c2b5570805583d782ffe123bcaa |
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. 120) then
204 map%source = 'NOAA GSD Rapid Refresh'
205 else
206 map%source = 'NOAA GSD'
207 endif
208 else if (icenter .eq. 60) then
209 map%source = 'NCAR'
210 else if (icenter .eq. 98) then
211 map%source = 'ECMWF'
212 else if (icenter .eq. 74 .or. icenter .eq. 75 ) then
213 map%source = 'UKMO'
214 else
215 map%source = 'unknown model and orig center'
216 end if
218 IPARM=KSEC1(7) ! Indicator of parameter
219 KTYPE=KSEC1(8) ! Indicator of type of level
221 ! print *,' IPARM, KTYPE, KSEC1(9)', iparm,ktype,ksec1(9)
223 IF(KTYPE.EQ.1) THEN
224 LVL1=KSEC1(9)
225 LVL2=-99
226 ELSEIF(KTYPE.EQ.100) THEN
227 LVL1=FLOAT(KSEC1(9)) * 100.
228 LVL2=-99
229 ELSEIF(KTYPE.EQ.101) THEN
230 LVL1=KSEC1(9)
231 LVL2=KSEC1(10)
232 ELSEIF(KTYPE.EQ.102) THEN
233 LVL1=KSEC1(9)
234 LVL2=-99
235 ELSEIF(KTYPE.EQ.103) THEN
236 LVL1=KSEC1(9)
237 LVL2=-99
238 ELSEIF(KTYPE.EQ.105) THEN
239 LVL1=KSEC1(9)
240 LVL2=-99
241 ELSEIF(KTYPE.EQ.107) THEN
242 LVL1=KSEC1(9)
243 LVL2=-99
244 ELSEIF(KTYPE.EQ.109) THEN
245 LVL1=KSEC1(9)
246 LVL2=-99
247 ELSEIF(KTYPE.EQ.111) THEN
248 LVL1=KSEC1(9)
249 LVL2=-99
250 ELSEIF(KTYPE.EQ.112) THEN ! Layer between two depths below surface
251 LVL1=KSEC1(9)
252 LVL2=KSEC1(10)
253 ELSEIF(KTYPE.EQ.113) THEN
254 LVL1=KSEC1(9)
255 LVL2=-99
256 ELSEIF(KTYPE.EQ.115) THEN
257 LVL1=KSEC1(9)
258 LVL2=-99
259 ELSEIF(KTYPE.EQ.117) THEN
260 LVL1=KSEC1(9)
261 LVL2=-99
262 ELSEIF(KTYPE.EQ.119) THEN
263 LVL1=KSEC1(9)
264 LVL2=-99
265 ELSEIF(KTYPE.EQ.125) THEN
266 LVL1=KSEC1(9)
267 LVL2=-99
268 ELSEIF(KTYPE.EQ.160) THEN
269 LVL1=KSEC1(9)
270 LVL2=-99
271 ELSEIF(KTYPE.EQ.200) THEN
272 LVL1=0
273 LVL2=-99
274 ELSEIF(KTYPE.EQ.201) THEN
275 LVL1=0
276 LVL2=-99
277 ELSE
278 LVL1=KSEC1(9)
279 LVL2=KSEC1(10)
280 ENDIF
282 ! Check to see that the combination of iparm, ktype, lvl1, and lvl2
283 ! match what has been requested in the Vtable. If not, set the field
284 ! name to NULL, meaning that we do not want to process this one.
286 field = 'NULL'
287 do i = 1, maxvar
288 if (gcode(i).eq.iparm) then
289 if (lcode(i).eq.ktype) then
290 if ((level1(i).eq.lvl1) .or. (level1(i) == splatcode) ) then
291 if (level2(i).eq.lvl2) then
292 field=namvar(i)
293 level = -999.
294 if (ktype.eq.100) then ! Pressure-level
295 level=lvl1
296 elseif (ktype.eq.102) then
297 level=201300.
298 elseif ((ktype.eq.116.and.lvl1.le.50.and.lvl2.eq.0) .or. &
299 (ktype.eq.105).or.(ktype.eq.1) .or. &
300 (ktype.eq.111).or.(ktype.eq.112) ) then
301 ! level=200100.
302 level = float(200000+iprty(i))
303 elseif (ktype.eq.109 .or. ktype.eq.107) then ! hybrid or sigma levels
304 level = lvl1
305 elseif (ktype.eq. 6 ) then ! max wind
306 level = 6.
307 elseif (ktype.eq. 7 ) then ! trop
308 level = 7.
309 elseif (ktype .eq. 160 ) then ! depth below sea-surface (m)
310 level = 201500.
311 elseif (ktype .eq. 237 .or. ktype .eq. 238 ) then ! depth of ocean layer
312 level = 201600.
313 elseif (ktype .eq. 200 ) then !column variable (TCDC,PWAT,etc.)
314 level = lvl1 !
315 endif
316 if (level .lt. -998. ) then
317 write(6,*) 'Could not find a level for this Vtable entry'
318 write(6,*) 'iparm = ',iparm,' ktype = ',ktype,' lvl1 = ',lvl1,' lvl2 = ',lvl2
319 write(6,*) 'Fix the Vtable or modify rd_grib1.F'
320 stop 'rd_grib1'
321 endif
322 endif
323 endif
324 endif
325 endif
326 enddo
328 if (field .eq. 'NULL') then
329 call deallogrib
330 return
331 endif
333 if ((field.eq.'WEASD').or.(field.eq.'SNOW')) then
334 level = level + ksec1(19)+1
335 endif
337 ! Build the 19-character date string, based on GRIB header date and time
338 ! information, including forecast time information:
340 ICC=KSEC1(22) ! CENTURY OF THE DATA
341 IYY=KSEC1(11) ! (TWO-DIGIT) YEAR OF THE DATA
342 MONTH=KSEC1(12) ! MONTH OF THE DATA
343 DAY=KSEC1(13) ! DAY OF THE DATA
344 HOUR=KSEC1(14) ! HOUR OF THE DATA
345 MINUTE=KSEC1(15) ! MINUTE OF THE DATA
346 SECOND=0
347 if (ksec1(19) == 3) then
348 FCST = (KSEC1(17) + KSEC1(18))/2
349 ! TEMPORARY AFWA FIX
350 ! elseif (ksec1(19) == 4 .or. ksec1(19) == 5) then
351 elseif (ksec1(19) == 4 .or. ksec1(19) == 5 .or. ksec1(19) == 7) then
352 FCST = KSEC1(18)
353 else
354 FCST = KSEC1(17)
355 endif
357 if (IYY.EQ.00) then
358 YEAR = ICC*100
359 else
360 YEAR = (ICC-1)*100 + IYY
361 endif
363 hdate(1:19) = ' '
364 call build_hdate(hdate,year,month,day,hour,minute,second)
366 call geth_newdate(hdate,hdate,3600*fcst)
368 ! Store information about the grid on which the data is.
369 ! This stuff gets stored in the MAP variable, as defined in module GRIDINFO
371 map%startloc = 'SWCORNER'
372 map%grid_wind = .true.
373 ! NCEP's grib1 messages (in GDS Octet 17, the Resolution and Component Flags)
374 ! all have '0' for the earth radius flag which the documentation (written by NCEP)
375 ! says is 6367.47, but they really use 6371.229. Hardcode it.
376 ! It's not clear what ECMWF uses. One place says 6367.47 and another 6371.229.
377 if ( index(map%source,'NCEP') .ne. 0 ) then
378 map%r_earth = 6371.229
379 else
380 map%r_earth = 6367.47
381 endif
383 if (ksec2(4).eq.0) then ! Lat/Lon grid
384 map%igrid = 0
385 map%nx = infogrid(1)
386 map%ny = infogrid(2)
387 map%dx = ginfo(8)
388 map%dy = ginfo(9)
389 map%lat1 = ginfo(3)
390 map%lon1 = ginfo(4)
391 ! If this is global data, then the dx and dy are more accurately
392 ! computed by the number of points than the 3 digits grib permits.
393 if ( ABS(map%nx * map%dx - 360.) .lt. 1 ) then
394 if ( ABS ( map%dx - (360./real(map%nx)) ) .gt. 0.00001 ) then
395 !print *,'old dx = ',ginfo(8)
396 map%dx = 360./real(map%nx)
397 !print *,'new dx = ',map%dx
398 endif
399 endif
400 if ( ABS((map%ny-1) * map%dy - 2.*abs(map%lat1)) .lt. 1 ) then
401 if ( ABS ( map%dy - (2.*abs(map%lat1)/real(map%ny-1)) ) .gt. 0.00001 ) then
402 !print *,'old dy = ',ginfo(9)
403 map%dy = 2.*abs(map%lat1)/real(map%ny-1)
404 !print *,'new dy = ',map%dy
405 endif
406 endif
407 write(tmp8,'(b8.8)') infogrid(5)
408 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
409 if (icenter .eq. 7 .and. KSEC1(5) .eq. 173 ) then ! correction for ncep grid 173
410 map%lat1 = 89.958333
411 map%lon1 = 0.041667
412 map%dx = 0.083333333 * sign(1.0,map%dx)
413 map%dy = 0.083333333 * sign(1.0,map%dy)
414 endif
415 ! correction for ncep grid 229 added 5/3/07 JFB
416 if (icenter .eq. 7 .and. KSEC1(5) .eq. 229 ) then
417 if (ginfo(3) .gt. 89. .and. ginfo(9) .gt. 0.) then
418 map%dy = -1. * map%dy
419 endif
420 endif
422 ! print *, "CE map stuff", map%igrid, map%nx, map%ny, map%dx, &
423 ! map%dy, map%lat1, map%lon1
425 elseif (ksec2(4).eq.1) then ! Mercator Grid
426 map%igrid = 1
427 map%nx = infogrid(1)
428 map%ny = infogrid(2)
429 map%dx = ginfo(8) ! km
430 map%dy = ginfo(9)
431 map%truelat1 = ginfo(5)
432 map%truelat2 = 0.
433 map%lov = 0.
434 map%lat1 = ginfo(3)
435 map%lon1 = ginfo(4)
436 write(tmp8,'(b8.8)') infogrid(5)
437 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
439 elseif (ksec2(4).eq.3) then ! Lambert Conformal Grid
440 map%igrid = 3
441 map%nx = infogrid(1)
442 map%ny = infogrid(2)
443 map%lov = ginfo(6)
444 map%truelat1 = ginfo(11)
445 map%truelat2 = ginfo(12)
446 map%dx = ginfo(7)
447 map%dy = ginfo(8)
448 map%lat1 = ginfo(3)
449 map%lon1 = ginfo(4)
450 write(tmp8,'(b8.8)') infogrid(5)
451 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
452 ! if (tmp8(2:2) .eq. '0') map%r_earth = 6367.47
454 elseif(ksec2(4).eq.4) then ! Gaussian Grid
455 map%igrid = 4
456 map%nx = infogrid(1)
457 map%ny = infogrid(2)
458 map%dx = ginfo(8)
459 ! map%dy = ginfo(19)
460 map%dy = real (infogrid(9))
461 map%lon1 = ginfo(4)
462 map%lat1 = ginfo(3)
463 write(tmp8,'(b8.8)') infogrid(5)
464 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
466 elseif (ksec2(4).eq.5) then ! Polar-Stereographic Grid.
467 map%igrid = 5
468 map%nx = infogrid(1)
469 map%ny = infogrid(2)
470 map%lov = ginfo(6)
471 map%truelat1 = 60.
472 map%truelat2 = 91.
473 map%dx = ginfo(7)
474 map%dy = ginfo(8)
475 map%lat1 = ginfo(3)
476 map%lon1 = ginfo(4)
477 write(tmp8,'(b8.8)') infogrid(5)
478 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
480 else
481 print*, 'Unknown Data Representation Type, ksec2(4)= ', ksec2(4)
482 stop 'rd_grib1'
483 endif
485 111 format(' igrid : ', i3, /, &
486 ' nx, ny : ', 2I4, /, &
487 ' truelat1, 2: ', 2F10.4, /, &
488 ' Center Lon : ', F10.4, /, &
489 ' LatLon(1,1): ', 2F10.4, /, &
490 ' DX, DY : ', F10.4, F10.4)
492 ! Special for NCEP/NCAR Reanalysis Project:
493 ! Throw out PSFC on lat/lon grid (save gaussian version)
494 if ((icenter.eq.7).and.(iprocess.eq.80)) then ! Careful! This combination may refer
495 ! to other products as well.
496 if ((field.eq.'PSFC').and.(ksec2(4).eq.0)) then
497 field='NULL'
498 call deallogrib
499 return
500 endif
501 endif
503 if (allocated(rdatarray)) deallocate(rdatarray)
504 allocate(rdatarray(map%nx * map%ny))
506 ! If nx=65535, assume the grid is a thinned grid.
507 ! Process only the NCEP grid IDs is 37 to 44.
508 if (map%nx.eq.65535) then
509 if ( (icenter .ne. 7) .or. (KSEC1(5).lt.37) .or. (KSEC1(5).gt.44) ) then
510 write(*,*) 'Originating center is ',icenter
511 write(*,*) 'Grid ID is ',KSEC1(5),' Only WAFS grids 37-44 are supported'
512 write(*,'(" ***** STOP in Subroutine RD_GRIB1.",//)')
513 stop
514 endif
515 lthinned = .TRUE.
516 map%nx = 73
517 map%dx = 1.25
518 else
519 lthinned = .FALSE.
520 endif
522 ! Unpack the 2D slab from the GRIB record, and put it in array rdatarray
524 if (lthinned) then
525 if (allocated(thinnedDataArray)) deallocate(thinnedDataArray)
526 allocate(thinnedDataArray(map%nx * map%ny))
527 call gribdata(thinnedDataArray,3447)
529 ! Calculate how many points for each latitude, and accumulate into array
530 if ((KSEC1(5).ge.37).and.(KSEC1(5).le.40)) then
531 ! Northern hemisphere:
532 npoints_acc(1)=0
533 npoints_acc(2)=73
534 do i=1,72
535 np = int(2.0+(90.0/1.25)*cos(i*1.25*3.1415926/180.0))
536 npoints_acc(i+2)=npoints_acc(i+1)+np
537 enddo
538 else
539 ! Southern Hemisphere:
540 npoints_acc(1)=0
541 npoints_acc(2)=2
542 do i=1,71
543 ii = 72-i
544 np = int(2.0+(90.0/1.25)*cos(ii*1.25*3.1415926/180.0))
545 npoints_acc(i+2)=npoints_acc(i+1)+np
546 enddo
547 npoints_acc(74) = npoints_acc(73) + 73
548 endif
550 ! for row number i (where i=1 is the southern edge of the grid)
551 ! npoints_acc(i+1)-npoints_acc(i) = number of points in this line
552 ! npoints_acc(i)+1 = index into thinned array for first point of line
554 do ny=1,73
555 np = npoints_acc(ny+1)-npoints_acc(ny) ! Number of points in this line.
556 do nx=1,73
557 ! Calulate the x index (mj) of thinned array (real value)
558 mj = (nx-1.0)*(np-1.0)/(72.0)
560 if (abs(mj - int(mj)) < 1.E-10) then
561 rdatarray((ny-1)*73+nx) = thinnedDataArray(npoints_acc(ny)+1+int(mj))
562 else
563 ! Get the 2 closest values from thinned array
564 Vb = thinnedDataArray(npoints_acc(ny)+1+int(mj))
565 Vc = thinnedDataArray(npoints_acc(ny)+1+int(mj)+1)
566 ! Get the next two closest, if available:
567 Va = -999999.
568 Vd = -999999.
569 if (mj > 1.0) then
570 Va = thinnedDataArray(npoints_acc(ny)+1+int(mj)-1)
571 endif
572 if (mj < np-2) then
573 Vd = thinnedDataArray(npoints_acc(ny)+1+int(mj)+2)
574 endif
576 if ((Va < -999998.) .or. (Vd < -999998.)) then
577 ! Use 2-point linear interpolation.
578 rdatarray((ny-1)*73+nx) = Vb*(int(mj)+1.0-mj) + Vc*(mj-int(mj))
579 else
580 ! Use 4-point overlapping parabolic interpolation.
581 xmj = mj - float(int(mj))
582 rdatarray((ny-1)*73+nx) = oned(xmj,Va,Vb,Vc,Vd)
583 endif
584 endif
585 enddo
586 enddo
587 else
588 call gribdata(rdatarray,map%nx*map%ny)
589 endif
591 ! Some grids are broken and need to be reordered (e.g. NCEP-II in 1997).
592 ! WPS assumes that the grids are ordered consistently with the start location.
594 call mprintf(.true.,DEBUG, &
595 "RD_GRIB1 icenter = %i , iprocess = %i , grid = %i",i1=icenter,i2=iprocess,i3=KSEC1(5))
596 if (icenter .eq. 7 .and. iprocess .eq. 0 .and. KSEC1(5) .eq. 2 ) then
597 call mprintf(.true.,DEBUG, &
598 "resetting NCEP2 dx and dy. If this is not NCEP2 data you must modify rd_grib1.f90")
599 call mprintf(.true.,DEBUG, &
600 "field = %s , dx = %f , dy = %f , i10 = %i",s1=field,f1=map%dx,f2=map%dy,i1=infogrid(10))
601 map%dx = 2.5
602 map%dy = -2.5
603 ! call reorder_it (rdatarray, map%nx, map%ny, map%dx, map%dy, infogrid(10))
604 endif
606 ! Deallocate a couple of arrays that may have been allocated by the
607 ! GRIB decoding routines.
609 call deallogrib
611 END subroutine rd_grib1
613 real function oned(x, a, b, c, d) Result (Answer)
614 implicit none
615 real :: x ! Proportion of the way between B and C. Between 0.0 and 1.0
616 real :: a, b, c, d
618 if (abs(x) < 1.E-10) then
619 Answer = B
620 return
621 endif
622 IF(abs(x-1.) < 1.E-10) then
623 Answer = C
624 return
625 endif
626 Answer = (1.0-X)*(B+X*(0.5*(C-A)+X*(0.5*(C+A)-B)))+X*(C+(1.0-X)*(0.5 &
627 *(B-D)+(1.0-X)*(0.5*(B+D)-C)))
628 end function oned