!*****************************************************************************!
! Subroutine RD_GRIB1 !
! !
! Purpose: !
! Read one record from the input GRIB file. Based on the information in !
! the GRIB header and the user-defined Vtable, decide whether the field in !
! the GRIB record is one to process or to skip. If the field is one we !
! want to keep, extract the data from the GRIB record, and pass the data !
! back to the calling routine. !
! !
! Argument list: !
! Input: !
! IUNIT : "Unit Number" to open and read from. Not really a Fortran !
! unit number, since we do not do Fortran I/O for the GRIB !
! files. Nor is it a UNIX File Descriptor returned from a C !
! OPEN statement. It is really just an array index to the !
! array (IUARR) where the UNIX File Descriptor values are !
! stored. !
! GRIBFLNM: File name to open, if it is not already open. !
! IUARR : Array to hold UNIX File descriptors retured from a C open !
! statement. If the value of IUARR(IUNIT) is zero, then the !
! file GRIBFLNM must be opened, and the value of IUARR(IUNIT) !
! becomes the UNIX File descriptor to read from. !
! DEBUG_LEVEL Integer for various amounts of printout. !
! EC_REC_LEN : Record length for EC ds113.0 files. !
! PMIN : Minimum pressure level (Pa) to process. !
! !
! Output: !
! LEVEL : The pressure-level (Pa) of the field to process. !
! FIELD : The field name of the field to process. NULL is returned !
! if we do not want to process the field we read. !
! HDATE : The 19-character date of the field to process. !
! IERR : Error flag: 0 - no error on read from GRIB file. !
! 1 - Hit the end of the GRIB file. !
! 2 - The file GRIBFLNM we tried to open does !
! not exist. !
! Externals !
! Module TABLE !
! Module GRIDINFO !
! Subroutine C_OPEN !
! Subroutine DEALLOGRIB !
! Subroutine GRIBGET !
! Subroutine GRIBHEADER !
! Subroutine GET_SEC1 !
! Subroutine GET_SEC2 !
! Subroutine GET_GRIDINFO !
! Subroutine BUILD_HDATE !
! Subroutine GETH_NEWDATE !
! Subroutine GRIBDATA !
! !
! Side Effects !
! File GRIBFLNM is opened, as necessary !
! !
! Variable MAP from module GRIDINFO is filled in. !
! !
! Numerous side effects from the GRIB-processing routines. !
! !
! Kevin W. Manning !
! NCAR/MMM !
! Summer, 1998, and continuing !
! SDG !
! !
!*****************************************************************************!
SUBROUTINE rd_grib1(IUNIT, gribflnm, level, field, hdate, &
ierr, iuarr, debug_level, ec_rec_len, pmin)
use table
use gridinfo
use datarray
use module_debug
implicit none
integer :: debug_level, ec_rec_len
integer :: iunit ! Array number in IUARR assigned to the C read pointer.
integer, dimension(100) :: KSEC1
integer, dimension(10) :: KSEC2
integer, dimension(40) :: infogrid
real, dimension(40) :: ginfo
!
!-----------------------------------------------------------------------
integer :: iparm, ktype
logical :: lopen
integer :: icenter, iprocess, iscan, ii, isb
integer year, month, day, hour, minute, second, icc, iyy
integer :: fcst
real :: level
character(LEN=*) :: field
character(LEN=132) :: gribflnm
character(LEN=8) :: tmp8
integer, dimension(255) :: iuarr
integer :: ierr, iostat, nunit
integer :: i, lvl2, lvl1
character(LEN=19) :: hdate
integer :: igherr
real :: pmin
! Variables for thinned grids:
logical :: lthinned = .FALSE.
real, allocatable, dimension(:) :: thinnedDataArray
integer, dimension(74) :: npoints_acc
real :: mj, xmj
integer :: np, ny, nx
real :: Va, Vb, Vc, Vd
real, external :: oned
ierr = 0
! If the file GRIBFLNM has not been opened, then IUARR(IUNIT) should be Zero.
! In this case, open the file GRIBFLNM, and store the UNIX File descriptor
! in to IUARR(IUNIT). This way, we will know what UNIX File descriptor to use
! next time we call this RD_GRIB subroutine.
!
if (iuarr(iunit).eq.0) then
if (debug_level.gt.0) then
call c_open(iunit, nunit, gribflnm, 1, ierr, 1)
else
call c_open(iunit, nunit, gribflnm, 1, ierr, -1)
endif
if (ierr.ne.0) then
call deallogrib
ierr = 2
return
endif
iuarr(iunit) = nunit
endif
! Read a single GRIB record, but do no unpacking now:
call gribget(iuarr(iunit), ierr, ec_rec_len)
if (ierr.ne.0) then
call mprintf(.true.,DEBUG,"RD_GRIB1 gribget read error, ierr = %i",i1=ierr)
call deallogrib
return
endif
!
! Unpack the header information:
!
call gribheader(debug_level,igherr,ec_rec_len)
if (igherr /= 0) then
field = "NULL"
call deallogrib
return
endif
!
! Copy header information to arrays KSEC1, KSEC2, INFOGRID, and GRIDINFO
!
call get_sec1(ksec1)
call get_sec2(ksec2)
call get_gridinfo(infogrid, ginfo)
icenter = KSEC1(3) ! Indicator of the source (center) of the data.
iprocess = KSEC1(4) ! Indicator of model (or whatever) which generated the data.
if (icenter.eq.7) then
if (iprocess.eq.83 .or. iprocess.eq.84) then
map%source = 'NCEP MESO NAM Model'
elseif (iprocess.eq.81) then
map%source = 'NCEP GFS Analysis'
elseif (iprocess.eq.82) then
map%source = 'NCEP GFS GDAS/FNL'
elseif (iprocess.eq.89) then
map%source = 'NCEP NMM '
elseif (iprocess.eq.96) then
map%source = 'NCEP GFS Model'
elseif (iprocess.eq.107) then
map%source = 'NCEP GEFS'
elseif (iprocess.eq.109) then
map%source = 'NCEP RTMA'
elseif (iprocess.eq.86 .or. iprocess.eq.100) then
map%source = 'NCEP RUC Model' ! 60 km
elseif (iprocess.eq.101) then
map%source = 'NCEP RUC Model' ! 40 km
elseif (iprocess.eq.105) then
map%source = 'NCEP RUC Model' ! 20 km
elseif (iprocess.eq.140) then
map%source = 'NCEP NARR'
elseif (iprocess.eq.195) then
map%source = 'NCEP CDAS2'
elseif (iprocess.eq.44) then
map%source = 'NCEP SST Analysis'
elseif (iprocess.eq.70) then
map%source = 'GFDL Hurricane Model'
elseif (iprocess.eq.129) then
map%source = 'NCEP GODAS'
elseif (iprocess.eq.25) then
map%source = 'NCEP SNOW COVER ANALYSIS'
else
map%source = 'unknown model from NCEP'
end if
! grid numbers only set for NCEP and AFWA models
write(tmp8,'("GRID ",i3)') KSEC1(5)
map%source(25:32) = tmp8
else if (icenter .eq. 57) then
if (iprocess .eq. 87) then
map%source = 'AFWA AGRMET'
else
map%source = 'AFWA'
endif
write(tmp8,'("GRID ",i3)') KSEC1(5)
map%source(25:32) = tmp8
else if (icenter .eq. 58) then
map%source = 'US Navy FNOC'
else if (icenter .eq. 59) then
if (iprocess .eq. 125) then
map%source = 'NOAA GSD Rapid Refresh'
else if (iprocess .eq. 105) then
map%source = 'NOAA GSD'
else
print *,'Unknown GSD source'
stop
endif
else if (icenter .eq. 60) then
map%source = 'NCAR'
else if (icenter .eq. 98) then
map%source = 'ECMWF'
else if (icenter .eq. 74 .or. icenter .eq. 75 ) then
map%source = 'UKMO'
else
map%source = 'unknown model and orig center'
end if
IPARM=KSEC1(7) ! Indicator of parameter
KTYPE=KSEC1(8) ! Indicator of type of level
! print *,' IPARM, KTYPE, KSEC1(9)', iparm,ktype,ksec1(9)
IF(KTYPE.EQ.1) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.100) THEN
LVL1=FLOAT(KSEC1(9)) * 100.
LVL2=-99
ELSEIF(KTYPE.EQ.101) THEN
LVL1=KSEC1(9)
LVL2=KSEC1(10)
ELSEIF(KTYPE.EQ.102) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.103) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.105) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.107) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.109) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.111) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.112) THEN ! Layer between two depths below surface
LVL1=KSEC1(9)
LVL2=KSEC1(10)
ELSEIF(KTYPE.EQ.113) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.115) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.117) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.119) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.125) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.160) THEN
LVL1=KSEC1(9)
LVL2=-99
ELSEIF(KTYPE.EQ.200) THEN
LVL1=0
LVL2=-99
ELSEIF(KTYPE.EQ.201) THEN
LVL1=0
LVL2=-99
ELSE
LVL1=KSEC1(9)
LVL2=KSEC1(10)
ENDIF
! Check to see that the combination of iparm, ktype, lvl1, and lvl2
! match what has been requested in the Vtable. If not, set the field
! name to NULL, meaning that we do not want to process this one.
field = 'NULL'
do i = 1, maxvar
if (gcode(i).eq.iparm) then
if (lcode(i).eq.ktype) then
if ((level1(i).eq.lvl1) .or. (level1(i) == splatcode) ) then
if (level2(i).eq.lvl2) then
field=namvar(i)
level = -999.
if (ktype.eq.100) then ! Pressure-level
level=lvl1
if ( level .lt. pmin ) field = 'NULL'
elseif (ktype.eq.102) then
level=201300.
elseif ((ktype.eq.116.and.lvl1.le.50.and.lvl2.eq.0) .or. &
(ktype.eq.105).or.(ktype.eq.1) .or. &
(ktype.eq.111).or.(ktype.eq.112) ) then
! level=200100.
level = float(200000+iprty(i))
elseif (ktype.eq.109 .or. ktype.eq.107) then ! hybrid or sigma levels
level = lvl1
elseif (ktype.eq. 6 ) then ! max wind
level = 6.
elseif (ktype.eq. 7 ) then ! trop
level = 7.
elseif (ktype .eq. 160 ) then ! depth below sea-surface (m)
level = 201500.
elseif (ktype .eq. 237 .or. ktype .eq. 238 ) then ! depth of ocean layer
level = 201600.
elseif (ktype .eq. 200 ) then !column variable (TCDC,PWAT,etc.)
level = lvl1 !
endif
if (level .lt. -998. ) then
write(6,*) 'Could not find a level for this Vtable entry'
write(6,*) 'iparm = ',iparm,' ktype = ',ktype,' lvl1 = ',lvl1,' lvl2 = ',lvl2
write(6,*) 'Fix the Vtable or modify rd_grib1.F'
stop 'rd_grib1'
endif
endif
endif
endif
endif
enddo
if (field .eq. 'NULL') then
call deallogrib
return
endif
if ((field.eq.'WEASD').or.(field.eq.'SNOW')) then
level = level + ksec1(19)+1
endif
! Build the 19-character date string, based on GRIB header date and time
! information, including forecast time information:
ICC=KSEC1(22) ! CENTURY OF THE DATA
IYY=KSEC1(11) ! (TWO-DIGIT) YEAR OF THE DATA
MONTH=KSEC1(12) ! MONTH OF THE DATA
DAY=KSEC1(13) ! DAY OF THE DATA
HOUR=KSEC1(14) ! HOUR OF THE DATA
MINUTE=KSEC1(15) ! MINUTE OF THE DATA
SECOND=0
if (ksec1(19) == 3) then
FCST = (KSEC1(17) + KSEC1(18))/2
! TEMPORARY AFWA FIX
! elseif (ksec1(19) == 4 .or. ksec1(19) == 5) then
elseif (ksec1(19) == 4 .or. ksec1(19) == 5 .or. ksec1(19) == 7) then
FCST = KSEC1(18)
else
FCST = KSEC1(17)
endif
! convert the fcst units to hours if necessary
if (ksec1(16) .eq. 254 ) then
fcst = fcst/3600.
elseif (ksec1(16) .eq. 0 ) then
fcst = fcst/60.
endif
if (IYY.EQ.00) then
YEAR = ICC*100
else
YEAR = (ICC-1)*100 + IYY
endif
hdate(1:19) = ' '
call build_hdate(hdate,year,month,day,hour,minute,second)
call geth_newdate(hdate,hdate,3600*fcst)
! Store information about the grid on which the data is.
! This stuff gets stored in the MAP variable, as defined in module GRIDINFO
map%startloc = 'SWCORNER'
map%grid_wind = .true.
! NCEP's grib1 messages (in GDS Octet 17, the Resolution and Component Flags)
! all have '0' for the earth radius flag which the documentation (written by NCEP)
! says is 6367.47, but they really use 6371.229. Hardcode it.
! It's not clear what ECMWF uses. One place says 6367.47 and another 6371.229.
if ( index(map%source,'NCEP') .ne. 0 ) then
map%r_earth = 6371.229
else
map%r_earth = 6367.47
endif
if (ksec2(4).eq.0) then ! Lat/Lon grid
map%igrid = 0
map%nx = infogrid(1)
map%ny = infogrid(2)
map%dx = ginfo(8)
map%dy = ginfo(9)
map%lat1 = ginfo(3)
map%lon1 = ginfo(4)
! If this is global data, then the dx and dy are more accurately
! computed by the number of points than the 3 digits grib permits.
if ( ABS(map%nx * map%dx - 360.) .lt. 2. ) then ! Check if it's a global grid
if ( ABS ( map%dx - (360./real(map%nx)) ) .gt. 0.00001 ) then
!print *,'old dx = ',ginfo(8)
map%dx = 360./real(map%nx)
!print *,'new dx = ',map%dx
endif
endif
if ( ABS((map%ny-1) * map%dy - 2.*abs(map%lat1)) .lt. 1. ) then
if ( ABS ( map%dy - (2.*abs(map%lat1)/real(map%ny-1)) ) .gt. 0.00001 ) then
!print *,'old dy = ',ginfo(9)
map%dy = 2.*abs(map%lat1)/real(map%ny-1)
!print *,'new dy = ',map%dy
endif
endif
write(tmp8,'(b8.8)') infogrid(5)
if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
if (icenter .eq. 7 .and. KSEC1(5) .eq. 173 ) then ! correction for ncep grid 173
map%lat1 = 89.958333
map%lon1 = 0.041667
map%dx = 0.083333333 * sign(1.0,map%dx)
map%dy = 0.083333333 * sign(1.0,map%dy)
endif
! correction for ncep grid 229 added 5/3/07 JFB
if (icenter .eq. 7 .and. KSEC1(5) .eq. 229 ) then
if (ginfo(3) .gt. 89. .and. ginfo(9) .gt. 0.) then
map%dy = -1. * map%dy
endif
endif
! print *, "CE map stuff", map%igrid, map%nx, map%ny, map%dx, &
! map%dy, map%lat1, map%lon1
elseif (ksec2(4).eq.1) then ! Mercator Grid
map%igrid = 1
map%nx = infogrid(1)
map%ny = infogrid(2)
map%dx = ginfo(8) ! km
map%dy = ginfo(9)
map%truelat1 = ginfo(5)
map%truelat2 = 0.
map%lov = 0.
map%lat1 = ginfo(3)
map%lon1 = ginfo(4)
write(tmp8,'(b8.8)') infogrid(5)
if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
elseif (ksec2(4).eq.3) then ! Lambert Conformal Grid
map%igrid = 3
map%nx = infogrid(1)
map%ny = infogrid(2)
map%lov = ginfo(6)
map%truelat1 = ginfo(11)
map%truelat2 = ginfo(12)
map%dx = ginfo(7)
map%dy = ginfo(8)
map%lat1 = ginfo(3)
map%lon1 = ginfo(4)
write(tmp8,'(b8.8)') infogrid(5)
if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
! if (tmp8(2:2) .eq. '0') map%r_earth = 6367.47
elseif(ksec2(4).eq.4) then ! Gaussian Grid
map%igrid = 4
map%nx = infogrid(1)
map%ny = infogrid(2)
map%dx = ginfo(8)
! map%dy = ginfo(19)
map%dy = real (infogrid(9))
map%lon1 = ginfo(4)
map%lat1 = ginfo(3)
write(tmp8,'(b8.8)') infogrid(5)
if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
! If this is global data, then the dx and dy are more accurately
! computed by the number of points than the 3 digits grib permits.
if ( ABS(map%nx * map%dx - 360.) .lt. 2. ) then ! Check if it's a global grid
if ( ABS ( map%dx - (360./real(map%nx)) ) .gt. 0.00001 ) then
! print *,'old dx = ',ginfo(8)
map%dx = 360./real(map%nx)
! print *,'new dx = ',map%dx
endif
endif
elseif (ksec2(4).eq.5) then ! Polar-Stereographic Grid.
map%igrid = 5
map%nx = infogrid(1)
map%ny = infogrid(2)
map%lov = ginfo(6)
map%truelat1 = 60.
map%truelat2 = 91.
map%dx = ginfo(7)
map%dy = ginfo(8)
map%lat1 = ginfo(3)
map%lon1 = ginfo(4)
write(tmp8,'(b8.8)') infogrid(5)
if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
else
print*, 'Unknown Data Representation Type, ksec2(4)= ', ksec2(4)
stop 'rd_grib1'
endif
111 format(' igrid : ', i3, /, &
' nx, ny : ', 2I4, /, &
' truelat1, 2: ', 2F10.4, /, &
' Center Lon : ', F10.4, /, &
' LatLon(1,1): ', 2F10.4, /, &
' DX, DY : ', F10.4, F10.4)
! Special for NCEP/NCAR Reanalysis Project:
! Throw out PSFC on lat/lon grid (save gaussian version)
if ((icenter.eq.7).and.(iprocess.eq.80)) then ! Careful! This combination may refer
! to other products as well.
if ((field.eq.'PSFC').and.(ksec2(4).eq.0)) then
field='NULL'
call deallogrib
return
endif
endif
if (allocated(rdatarray)) deallocate(rdatarray)
allocate(rdatarray(map%nx * map%ny))
! If nx=65535, assume the grid is a thinned grid.
! Process only the NCEP grid IDs is 37 to 44.
if (map%nx.eq.65535) then
if ( (icenter .ne. 7) .or. (KSEC1(5).lt.37) .or. (KSEC1(5).gt.44) ) then
write(*,*) 'Originating center is ',icenter
write(*,*) 'Grid ID is ',KSEC1(5),' Only WAFS grids 37-44 are supported'
write(*,'(" ***** STOP in Subroutine RD_GRIB1.",//)')
stop
endif
lthinned = .TRUE.
map%nx = 73
map%dx = 1.25
else
lthinned = .FALSE.
endif
! Unpack the 2D slab from the GRIB record, and put it in array rdatarray
if (lthinned) then
if (allocated(thinnedDataArray)) deallocate(thinnedDataArray)
allocate(thinnedDataArray(map%nx * map%ny))
call gribdata(thinnedDataArray,3447)
! Calculate how many points for each latitude, and accumulate into array
if ((KSEC1(5).ge.37).and.(KSEC1(5).le.40)) then
! Northern hemisphere:
npoints_acc(1)=0
npoints_acc(2)=73
do i=1,72
np = int(2.0+(90.0/1.25)*cos(i*1.25*3.1415926/180.0))
npoints_acc(i+2)=npoints_acc(i+1)+np
enddo
else
! Southern Hemisphere:
npoints_acc(1)=0
npoints_acc(2)=2
do i=1,71
ii = 72-i
np = int(2.0+(90.0/1.25)*cos(ii*1.25*3.1415926/180.0))
npoints_acc(i+2)=npoints_acc(i+1)+np
enddo
npoints_acc(74) = npoints_acc(73) + 73
endif
! for row number i (where i=1 is the southern edge of the grid)
! npoints_acc(i+1)-npoints_acc(i) = number of points in this line
! npoints_acc(i)+1 = index into thinned array for first point of line
do ny=1,73
np = npoints_acc(ny+1)-npoints_acc(ny) ! Number of points in this line.
do nx=1,73
! Calulate the x index (mj) of thinned array (real value)
mj = (nx-1.0)*(np-1.0)/(72.0)
if (abs(mj - int(mj)) < 1.E-10) then
rdatarray((ny-1)*73+nx) = thinnedDataArray(npoints_acc(ny)+1+int(mj))
else
! Get the 2 closest values from thinned array
Vb = thinnedDataArray(npoints_acc(ny)+1+int(mj))
Vc = thinnedDataArray(npoints_acc(ny)+1+int(mj)+1)
! Get the next two closest, if available:
Va = -999999.
Vd = -999999.
if (mj > 1.0) then
Va = thinnedDataArray(npoints_acc(ny)+1+int(mj)-1)
endif
if (mj < np-2) then
Vd = thinnedDataArray(npoints_acc(ny)+1+int(mj)+2)
endif
if ((Va < -999998.) .or. (Vd < -999998.)) then
! Use 2-point linear interpolation.
rdatarray((ny-1)*73+nx) = Vb*(int(mj)+1.0-mj) + Vc*(mj-int(mj))
else
! Use 4-point overlapping parabolic interpolation.
xmj = mj - float(int(mj))
rdatarray((ny-1)*73+nx) = oned(xmj,Va,Vb,Vc,Vd)
endif
endif
enddo
enddo
else
call gribdata(rdatarray,map%nx*map%ny)
endif
! Some grids are broken and need to be reordered (e.g. NCEP-II in 1997).
! WPS assumes that the grids are ordered consistently with the start location.
call mprintf(.true.,DEBUG, &
"RD_GRIB1 icenter = %i , iprocess = %i , grid = %i",i1=icenter,i2=iprocess,i3=KSEC1(5))
if (icenter .eq. 7 .and. iprocess .eq. 0 .and. KSEC1(5) .eq. 2 ) then
call mprintf(.true.,DEBUG, &
"resetting NCEP2 dx and dy. If this is not NCEP2 data you must modify rd_grib1.f90")
call mprintf(.true.,DEBUG, &
"field = %s , dx = %f , dy = %f , i10 = %i",s1=field,f1=map%dx,f2=map%dy,i1=infogrid(10))
map%dx = 2.5
map%dy = -2.5
! call reorder_it (rdatarray, map%nx, map%ny, map%dx, map%dy, infogrid(10))
endif
! Deallocate a couple of arrays that may have been allocated by the
! GRIB decoding routines.
call deallogrib
END subroutine rd_grib1
real function oned(x, a, b, c, d) Result (Answer)
implicit none
real :: x ! Proportion of the way between B and C. Between 0.0 and 1.0
real :: a, b, c, d
if (abs(x) < 1.E-10) then
Answer = B
return
endif
IF(abs(x-1.) < 1.E-10) then
Answer = C
return
endif
Answer = (1.0-X)*(B+X*(0.5*(C-A)+X*(0.5*(C+A)-B)))+X*(C+(1.0-X)*(0.5 &
*(B-D)+(1.0-X)*(0.5*(B+D)-C)))
end function oned