From a60c4674697f658a0ad7213f385ae271bb8d64ad Mon Sep 17 00:00:00 2001 From: Jim Bresch Date: Sat, 13 Jan 2018 22:37:55 -0700 Subject: [PATCH] Add capability to ungrib RDA's ECMWF ds113.1 (9km) output. The high-resolution EC output are stored in a non-standard grib1 file. The record-length is larger than can be stored by the 3-byte header record, so EC created a hack to the header which breaks ungrib. Decoding it properly would require reverse-engineering NCEP's wgrib C code and adding it to ungrib. The simplest work-around is to use wgrib to determine the record-length and feed that to ungrib. This is done via the new ec_rec_len namelist variable (default of zero). Most of the code changes are to pass the ec_rec_len down to the low-level grib routines that need the record length. The value for ec113.1 is 26214508 and should be constant. If they create an even higher-resolution dataset then we'll feed that to ungrib using ec_rec_len. --- ungrib/src/g1print.F | 11 +++++++---- ungrib/src/gribcode.F | 46 +++++++++++++++++++++++++++------------------- ungrib/src/rd_grib1.F | 8 ++++---- ungrib/src/read_namelist.F | 7 ++++--- ungrib/src/ungrib.F | 20 +++++++++++++------- 5 files changed, 55 insertions(+), 37 deletions(-) diff --git a/ungrib/src/g1print.F b/ungrib/src/g1print.F index 8f5ba15..d0ceda6 100644 --- a/ungrib/src/g1print.F +++ b/ungrib/src/g1print.F @@ -32,7 +32,10 @@ program g1print logical :: idb = .FALSE. integer :: year character(len=5) :: field + integer :: ec_rec_len + ec_rec_len = 0 +! ec_rec_len = 26214508 flnm = ' ' call parse_args(ierr, a1='v', l1=ivb, a2='V', l2=idb, hlast=flnm) if (ierr.ne.0) then @@ -58,13 +61,13 @@ program g1print endif irec = 0 - call gribget(munit, ierr) + call gribget(munit, ierr, ec_rec_len) do while (ierr.eq.0) irec = irec + 1 - call gribheader(0,igherr) + call gribheader(0,igherr, ec_rec_len) if (igherr /= 0) then call deallogrib - call gribget(munit, ierr) + call gribget(munit, ierr, ec_rec_len) cycle endif @@ -111,7 +114,7 @@ program g1print call deallogrib - call gribget(munit, ierr) + call gribget(munit, ierr, ec_rec_len) enddo if (ierr.eq.1) write(*,'(/,"***** End-Of-File on C unit ", I3,/)') munit call c_close( munit, 0, ierr) diff --git a/ungrib/src/gribcode.F b/ungrib/src/gribcode.F index 40c3e13..26baf20 100644 --- a/ungrib/src/gribcode.F +++ b/ungrib/src/gribcode.F @@ -14,11 +14,11 @@ ! ! ! The main user interfaces are: ! ! ! -! SUBROUTINE GRIBGET(NUNIT, IERR) ! +! SUBROUTINE GRIBGET(NUNIT, IERR, ec_rec_len) ! ! Read a single GRIB record from UNIX file-descriptor NUNIT into array ! ! GREC. No unpacking of any header or data values is performed. ! ! ! -! SUBROUTINE GRIBREAD(NUNIT, DATA, NDATA, IERR) ! +! SUBROUTINE GRIBREAD(NUNIT, DATA, NDATA, IERR, ec_rec_len) ! ! Read a single GRIB record from UNIX file-descriptor NUNIT, and unpack ! ! all header and data values into the appropriate arrays. ! ! ! @@ -251,7 +251,7 @@ contains !=============================================================================! !=============================================================================! ! - integer function gribsize(trec, ilen, ierr) + integer function gribsize(trec, ilen, ierr, ec_rec_len) !-----------------------------------------------------------------------------! ! Return the size of a single GRIB record. ! ! ! @@ -272,7 +272,7 @@ contains ! ! !-----------------------------------------------------------------------------! implicit none - integer :: ilen + integer :: ilen, ec_rec_len integer, dimension(ilen) :: trec integer :: isz0 = 32 integer :: isz1 = 0 @@ -295,7 +295,11 @@ contains if (ied.eq.1) then ! Grib1 ! Find the size of the whole GRIB record - call gbyte_g1(trec, gribsize, 32, 24) + if ( ec_rec_len .ne. 0 ) then + gribsize = ec_rec_len + else + call gbyte_g1(trec, gribsize, 32, 24) + endif ! GRIB Edition 0 does not include the total size, so we have to sum up ! ! the sizes of the individual sections ! @@ -350,7 +354,7 @@ contains !=============================================================================! !=============================================================================! ! - subroutine findgrib(nunit, isize, ierr) + subroutine findgrib(nunit, isize, ierr, ec_rec_len) !-----------------------------------------------------------------------------! ! ! @@ -382,7 +386,7 @@ contains ! ! !-----------------------------------------------------------------------------! implicit none - integer, intent(in) :: nunit + integer, intent(in) :: nunit, ec_rec_len integer, intent(out) :: isize integer, intent(out) :: ierr @@ -446,7 +450,7 @@ contains #ifdef BYTESWAP call swap4(trec, isz) #endif - isize = gribsize(trec, isz, ierr) + isize = gribsize(trec, isz, ierr, ec_rec_len) end subroutine findgrib ! @@ -763,7 +767,7 @@ end module module_grib !=============================================================================! !=============================================================================! ! -subroutine gribget(nunit, ierr) +subroutine gribget(nunit, ierr, ec_rec_len) use module_grib !-----------------------------------------------------------------------------! ! ! @@ -790,12 +794,12 @@ subroutine gribget(nunit, ierr) implicit none - integer :: nunit + integer :: nunit, ec_rec_len integer :: ierr integer :: isz, isize ! Position the file pointer at the beginning of the GRIB record. - call findgrib(nunit, isize, ierr) + call findgrib(nunit, isize, ierr, ec_rec_len) if (ierr.ne.0) return ! Allocate the GREC array to be able to hold the data @@ -821,7 +825,7 @@ end subroutine gribget !=============================================================================! !=============================================================================! ! -subroutine gribread(nunit, data, ndata, debug_level, ierr) +subroutine gribread(nunit, data, ndata, debug_level, ierr, ec_rec_len) !-----------------------------------------------------------------------------! ! Read one grib record, unpack the header and data information. ! ! ! @@ -852,7 +856,7 @@ subroutine gribread(nunit, data, ndata, debug_level, ierr) implicit none - integer :: nunit + integer :: nunit, ec_rec_len integer :: debug_level integer :: ierr real, allocatable, dimension(:) :: datarray @@ -863,12 +867,12 @@ subroutine gribread(nunit, data, ndata, debug_level, ierr) ierr = 0 - call gribget(nunit, ierr) + call gribget(nunit, ierr, ec_rec_len) if (ierr.ne.0) return ! Unpack the header information - call gribheader(debug_level,ierr) + call gribheader(debug_level,ierr, ec_rec_len) ! Determine the size of the data array from the information in the header, ! and allocate the array DATARRAY to hold that data. @@ -1250,14 +1254,14 @@ end subroutine glccone !=============================================================================! !=============================================================================! ! -subroutine gribheader(debug_level,ierr) +subroutine gribheader(debug_level,ierr, ec_rec_len) ! ! IERR non-zero means there was a problem unpacking the grib header. ! use module_grib implicit none integer :: debug_level - integer :: ierr + integer :: ierr, ec_rec_len integer, parameter :: nsec1 = 24 @@ -1284,10 +1288,14 @@ subroutine gribheader(debug_level,ierr) ! Section 0: sec0(2) = ied if (ied.eq.1) then - call gbyte_g1(grec, sec0(1), 32, 24) + if ( ec_rec_len .ne. 0 ) then + sec0(1) = ec_rec_len + else + call gbyte_g1(grec, sec0(1), 32, 24) + endif iskip = 64 elseif (ied.eq.0) then - sec0(1) = gribsize(grec,200, ierr) + sec0(1) = gribsize(grec,200, ierr, 0) iskip = 32 endif diff --git a/ungrib/src/rd_grib1.F b/ungrib/src/rd_grib1.F index b5b92cd..114bf08 100644 --- a/ungrib/src/rd_grib1.F +++ b/ungrib/src/rd_grib1.F @@ -60,7 +60,7 @@ ! ! !*****************************************************************************! SUBROUTINE rd_grib1(IUNIT, gribflnm, level, field, hdate, & - ierr, iuarr, debug_level) + ierr, iuarr, debug_level, ec_rec_len) use table use gridinfo use datarray @@ -68,7 +68,7 @@ SUBROUTINE rd_grib1(IUNIT, gribflnm, level, field, hdate, & implicit none - integer :: debug_level + 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 @@ -124,7 +124,7 @@ SUBROUTINE rd_grib1(IUNIT, gribflnm, level, field, hdate, & ! Read a single GRIB record, but do no unpacking now: - call gribget(iuarr(iunit), ierr) + 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) @@ -134,7 +134,7 @@ SUBROUTINE rd_grib1(IUNIT, gribflnm, level, field, hdate, & ! ! Unpack the header information: ! - call gribheader(debug_level,igherr) + call gribheader(debug_level,igherr,ec_rec_len) if (igherr /= 0) then field = "NULL" call deallogrib diff --git a/ungrib/src/read_namelist.F b/ungrib/src/read_namelist.F index b931c7c..4d8133d 100644 --- a/ungrib/src/read_namelist.F +++ b/ungrib/src/read_namelist.F @@ -1,6 +1,6 @@ subroutine read_namelist(hstart, hend, delta_time, ntimes,& ordered_by_date, debug_level, out_format, prefix, & - add_lvls, new_plvl_in, interp_type) + add_lvls, new_plvl_in, interp_type, ec_rec_len) use misc_definitions_module use module_debug @@ -15,7 +15,7 @@ subroutine read_namelist(hstart, hend, delta_time, ntimes,& integer :: debug_level real, dimension(:) :: new_plvl_in logical :: add_lvls - integer :: interp_type + integer :: interp_type, ec_rec_len integer :: ierr integer :: idts @@ -69,7 +69,7 @@ subroutine read_namelist(hstart, hend, delta_time, ntimes,& namelist /ungrib/ out_format, & ordered_by_date, prefix, & - add_lvls, new_plvl, interp_type + add_lvls, new_plvl, interp_type, ec_rec_len allocate(new_plvl(size(new_plvl_in))) @@ -97,6 +97,7 @@ subroutine read_namelist(hstart, hend, delta_time, ntimes,& add_lvls = .false. new_plvl = -99999. interp_type = 0 + ec_rec_len = 0 ! Start routine: diff --git a/ungrib/src/ungrib.F b/ungrib/src/ungrib.F index 9dd7ed7..46599b5 100644 --- a/ungrib/src/ungrib.F +++ b/ungrib/src/ungrib.F @@ -57,7 +57,7 @@ program ungrib interface subroutine read_namelist(hstart, hend, delta_time, ntimes,& ordered_by_date, debug_level, out_format, prefix, & - add_lvls, new_plvl, interp_type) + add_lvls, new_plvl, interp_type, ec_rec_len) use misc_definitions_module @@ -70,7 +70,7 @@ program ungrib character(len=MAX_FILENAME_LEN) :: prefix logical :: add_lvls real, dimension(:) :: new_plvl - integer :: interp_type + integer :: interp_type, ec_rec_len end subroutine read_namelist end interface @@ -84,7 +84,7 @@ program ungrib real , dimension(maxlvl) :: plvl, new_plvl integer :: iplvl logical :: add_lvls - integer :: interp_type + integer :: interp_type, ec_rec_len integer :: nlvl @@ -116,7 +116,7 @@ program ungrib call read_namelist(hstart, hend, interval, ntimes, & ordered_by_date, debug_level, out_format, prefix, & - add_lvls, new_plvl, interp_type) + add_lvls, new_plvl, interp_type, ec_rec_len) call mprintf(.true.,INFORM,"Interval value: %i seconds or %f hours", & i1=interval, f1=float(interval)/3600.) @@ -127,8 +127,14 @@ program ungrib ! ----------------- ! Determine GRIB Edition number grib_version=0 - call edition_num(nunit1, gribflnm, grib_version, ierr) - call mprintf((ierr.eq.3),ERROR,"GRIB file problem") +! NCAR/RDA's EC 113 data are grib1 with a non-standard record length + if ( ec_rec_len .eq. 0 ) then + call edition_num(nunit1, gribflnm, grib_version, ierr) + call mprintf((ierr.eq.3),ERROR,"GRIB file problem") + else + grib_version=1 + ierr=0 + endif if (grib_version.eq.2) then vtable_columns=12 #if defined (USE_PNG) && (USE_JPEG2000) @@ -226,7 +232,7 @@ program ungrib "flnm = %s",s1=gribflnm) ! Read one record at a time from GRIB1 (and older Editions) call rd_grib1(nunit1, gribflnm, level, field, & - hdate, ierr, iuarr, debug_level) + hdate, ierr, iuarr, debug_level, ec_rec_len) else ! Read one file of records from GRIB2. -- 2.11.4.GIT