module grib2_read_mod
!
! module: functions to read grib2 files
!
! Ming Hu
!
! program history log:
! 2017-04-10 Hu initial build
!
! Subroutines Included:
! sub drefresh_cldsurf - initialize RR related variables to default
!
use netcdf
implicit none
! set default to private
private
! set subroutines to public
public :: read_grib2_head_dim
public :: read_grib2_head_time
public :: read_grib2_sngle
contains
subroutine read_grib2_head_dim(filename,nx,ny,rlonmin,rlatmax,rdx,rdy,kgds)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_grib2 read grib2 file head dimension information
! prgmmr: Ming Hu org: GSD date: 2015-05-20
!
! abstract: read grib2 file head
!
!
! program history log:
! 2015-05-20 Hu, initial documentation
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: Zeus
!
!$$$ end documentation block
! use gridmod, only: idsl5,regional
use grib_mod
implicit none
character*100,intent(in) :: filename
integer, intent(out) :: nx,ny
real, intent(out) :: rlonmin,rlatmax
real*8, intent(out) :: rdx,rdy
integer, intent(out) :: KGDS(200)
!
integer :: nz
!
type(gribfield) :: gfld
logical :: expand=.true.
integer :: ifile
character(len=1),allocatable,dimension(:) :: cgrib
integer,parameter :: msk1=32000
integer :: lskip, lgrib,iseek
integer :: currlen
integer :: icount , lengrib
integer :: listsec0(3)
integer :: listsec1(13)
integer year, month, day, hour, minute, second, fcst
integer :: numfields,numlocal,maxlocal,ierr
integer :: grib_edition
integer :: itot
! real :: dx,dy,lat1,lon1
real :: scale_factor
!
!
integer :: nn,n,j,iret
real :: fldmax,fldmin,sum
!
!
scale_factor=1.0e6
ifile=10
loopfile: do nn=1,1
! write(6,*) 'read in grib2 file head', trim(filename)
lskip=0
lgrib=0
iseek=0
icount=0
itot=0
currlen=0
! Open GRIB2 file
call baopenr(ifile,trim(filename),iret)
if (iret.eq.0) then
VERSION: do
! Search opend file for the next GRIB2 messege (record).
call skgb(ifile,iseek,msk1,lskip,lgrib)
! Check for EOF, or problem
if (lgrib.eq.0) then
exit
endif
! Check size, if needed allocate more memory.
if (lgrib.gt.currlen) then
if (allocated(cgrib)) deallocate(cgrib)
allocate(cgrib(lgrib))
currlen=lgrib
endif
! Read a given number of bytes from unblocked file.
call baread(ifile,lskip,lgrib,lengrib,cgrib)
if(lgrib.ne.lengrib) then
write(*,*) 'ERROR, read_grib2 lgrib ne lengrib', &
lgrib,lengrib
stop 1234
endif
iseek=lskip+lgrib
icount=icount+1
! Unpack GRIB2 field
call gb_info(cgrib,lengrib,listsec0,listsec1, &
numfields,numlocal,maxlocal,ierr)
if(ierr.ne.0) then
write(6,*) 'Error querying GRIB2 message',ierr
stop
endif
itot=itot+numfields
grib_edition=listsec0(2)
if (grib_edition.ne.2) then
exit VERSION
endif
! write(*,*) 'listsec0=',listsec0
! write(*,*) 'listsec1=',listsec1
! write(*,*) 'numfields=',numfields
! get information form grib2 file
n=1
call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
year =gfld%idsect(6) !(FOUR-DIGIT) YEAR OF THE DATA
month =gfld%idsect(7) ! MONTH OF THE DATA
day =gfld%idsect(8) ! DAY OF THE DATA
hour =gfld%idsect(9) ! HOUR OF THE DATA
minute=gfld%idsect(10) ! MINUTE OF THE DATA
second=gfld%idsect(11) ! SECOND OF THE DATA
! write(*,*) 'year,month,day,hour,minute,second='
! write(*,*) year,month,day,hour,minute,second
! write(*,*) 'source center =',gfld%idsect(1)
! write(*,*) 'Indicator of model =',gfld%ipdtmpl(5)
! write(*,*) 'observation level (m)=',gfld%ipdtmpl(12)
! write(*,*) 'map projection=',gfld%igdtnum
if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical
! Equidistant
nx = gfld%igdtmpl(8)
ny = gfld%igdtmpl(9)
nz = 1
rdx = gfld%igdtmpl(17)/scale_factor
rdy = gfld%igdtmpl(18)/scale_factor
rlatmax = gfld%igdtmpl(12)/scale_factor
rlonmin = gfld%igdtmpl(13)/scale_factor
! write(*,*) 'nx,ny=',nx,ny
! write(*,*) 'dx,dy=',rdx,rdy
! write(*,*) 'lat1,lon1=',rlatmax,rlonmin
elseif (gfld%igdtnum.eq.20) then ! Polar stereographic projection
nx = gfld%igdtmpl(8)
ny = gfld%igdtmpl(9)
nz = 1
rdx = gfld%igdtmpl(15)/scale_factor
rdy = gfld%igdtmpl(16)/scale_factor
rlatmax = gfld%igdtmpl(10)/scale_factor
rlonmin = gfld%igdtmpl(11)/scale_factor
! write(*,*) gfld%igdtmpl(1:18)
!
! fill in this one to use gdswiz
! POLAR STEREOGRAPHIC GRIDS
! KGDS(i)
! (2) - N(I) NR POINTS ALONG LAT CIRCLE
! (3) - N(J) NR POINTS ALONG LON CIRCLE
! (4) - LA(1) LATITUDE OF ORIGIN
! (5) - LO(1) LONGITUDE OF ORIGIN
! (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17)
! (7) - LOV GRID ORIENTATION
! (8) - DX - X DIRECTION INCREMENT
! (9) - DY - Y DIRECTION INCREMENT
! (10) - PROJECTION CENTER FLAG
! (11) - SCANNING MODE (RIGHT ADJ COPY OF OCTET 28)
KGDS=0
KGDS(1)=5
KGDS(2) = nx
KGDS(3) = ny
KGDS(4) = gfld%igdtmpl(10)/1000
KGDS(5) = gfld%igdtmpl(11)/1000
KGDS(6) = 72
KGDS(7) = gfld%igdtmpl(14)/1000
KGDS(8) = gfld%igdtmpl(15)/1000
KGDS(9) = gfld%igdtmpl(16)/1000
KGDS(10) = gfld%igdtmpl(17)
KGDS(11) = gfld%igdtmpl(18)
! write(*,*) KGDS(1:11)
else
write(*,*) 'unknown projection'
stop 1235
endif
call gf_free(gfld)
enddo VERSION ! skgb
endif
CALL BACLOSE(ifile,ierr)
nullify(gfld%local)
if (allocated(cgrib)) deallocate(cgrib)
enddo loopfile
return
end subroutine read_grib2_head_dim
subroutine read_grib2_sngle(filename,ntot,ice,snow,lb)
!subroutine read_grib2_sngle(filename,ntot,var)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_grib2 read grib2 file
! prgmmr: Ming Hu org: GSD date: 2015-05-20
!
! abstract: read grib2 file
!
!
! program history log:
! 2015-05-20 parrish, initial documentation
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: Zeu
!
!$$$ end documentation block
! use gridmod, only: idsl5,regional
use grib_mod
implicit none
character*100,intent(in) :: filename
integer, intent(in) :: ntot
real, intent(out) :: ice(ntot)
real, intent(out) :: snow(ntot)
logical*1,intent(out) :: lb(ntot)
!
integer :: height
!
type(gribfield) :: gfld
logical :: expand=.true.
integer :: ifile
character(len=1),allocatable,dimension(:) :: cgrib
integer,parameter :: msk1=32000
integer :: lskip, lgrib,iseek
integer :: currlen
integer :: icount , lengrib
integer :: listsec0(3)
integer :: listsec1(13)
integer year, month, day, hour, minute, second, fcst
integer :: numfields,numlocal,maxlocal,ierr
integer :: grib_edition
integer :: itot
integer :: nx,ny
real :: dx,dy,lat1,lon1
real :: scale_factor
!
!
integer :: paramid
integer :: nn,n,j,iret
real :: fldmax,fldmin,sum
!
!
scale_factor=1.0e6
ifile=12
loopfile: do nn=1,1
! write(6,*) 'read mosaic in grib2 file ', trim(filename)
lskip=0
lgrib=0
iseek=0
icount=0
itot=0
currlen=0
! Open GRIB2 file
call baopenr(ifile,trim(filename),iret)
if (iret.eq.0) then
VERSION: do
! Search opend file for the next GRIB2 messege (record).
call skgb(ifile,iseek,msk1,lskip,lgrib)
! Check for EOF, or problem
if (lgrib.eq.0) then
exit
endif
! Check size, if needed allocate more memory.
if (lgrib.gt.currlen) then
if (allocated(cgrib)) deallocate(cgrib)
allocate(cgrib(lgrib))
currlen=lgrib
endif
! Read a given number of bytes from unblocked file.
call baread(ifile,lskip,lgrib,lengrib,cgrib)
if(lgrib.ne.lengrib) then
write(*,*) 'ERROR, read_grib2 lgrib ne lengrib', &
lgrib,lengrib
stop 1234
endif
iseek=lskip+lgrib
icount=icount+1
! Unpack GRIB2 field
call gb_info(cgrib,lengrib,listsec0,listsec1, &
numfields,numlocal,maxlocal,ierr)
if(ierr.ne.0) then
write(6,*) 'Error querying GRIB2 message',ierr
stop
endif
itot=itot+numfields
grib_edition=listsec0(2)
if (grib_edition.ne.2) then
exit VERSION
endif
! write(*,*) 'listsec0=',listsec0
! write(*,*) 'listsec1=',listsec1
! write(*,*) 'numfields=',numfields
! get information form grib2 file
n=1
call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
year =gfld%idsect(6) !(FOUR-DIGIT) YEAR OF THE DATA
month =gfld%idsect(7) ! MONTH OF THE DATA
day =gfld%idsect(8) ! DAY OF THE DATA
hour =gfld%idsect(9) ! HOUR OF THE DATA
minute=gfld%idsect(10) ! MINUTE OF THE DATA
second=gfld%idsect(11) ! SECOND OF THE DATA
! write(*,*) 'year,month,day,hour,minute,second='
! write(*,*) year,month,day,hour,minute,second
! write(*,*) 'source center =',gfld%idsect(1)
! write(*,*) 'Indicator of model =',gfld%ipdtmpl(5)
! write(*,*) 'observation level (m)=',gfld%ipdtmpl(12)
! write(*,*) 'map projection=',gfld%igdtnum
height=gfld%ipdtmpl(12)
if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical
! Equidistant
nx = gfld%igdtmpl(8)
ny = gfld%igdtmpl(9)
dx = gfld%igdtmpl(17)/scale_factor
dy = gfld%igdtmpl(18)/scale_factor
lat1 = gfld%igdtmpl(12)/scale_factor
lon1 = gfld%igdtmpl(13)/scale_factor
! write(*,*) 'nx,ny=',nx,ny
! write(*,*) 'dx,dy=',dx,dy
! write(*,*) 'lat1,lon1=',lat1,lon1
elseif (gfld%igdtnum.eq.20) then ! Polar stereographic projection
nx = gfld%igdtmpl(8)
ny = gfld%igdtmpl(9)
dx = gfld%igdtmpl(15)/scale_factor
dy = gfld%igdtmpl(16)/scale_factor
lat1 = gfld%igdtmpl(10)/scale_factor
lon1 = gfld%igdtmpl(11)/scale_factor
else
write(*,*) 'unknown projection'
stop 1235
endif
! write(*,*) 'Product Definition Template Number=',gfld%ipdtnum
! write(*,*) 'Number of elements=',gfld%ipdtlen
! write(*,*) 'specified Product Definition Template=', &
! gfld%ipdtmpl(1:gfld%ipdtlen)
paramid=gfld%ipdtmpl(2)
call gf_free(gfld)
! Continue to unpack GRIB2 field.
! write(*,*) 'number of fields =',numfields
NUM_FIELDS: do n = 1, numfields
! e.g. U and V would =2, otherwise its usually =1
call gf_getfld(cgrib,lengrib,n,.true.,expand,gfld,ierr)
if (ierr.ne.0) then
write(*,*) ' ERROR extracting field gf_getfld = ',ierr
cycle
endif
! write(*,*) 'gfld%ndpts=',n,gfld%ndpts
! write(*,*) 'gfld%unpacked=',n,gfld%unpacked
fldmax=gfld%fld(1)
fldmin=gfld%fld(1)
sum=gfld%fld(1)
if(ntot .lt. gfld%ndpts) then
write(*,*) 'Error, wrong dimension ',ntot, gfld%ndpts
stop 1234
endif
! do j=1,gfld%ndpts
if(paramid==0) then
do j=1,ntot
ice(j)=gfld%fld(j)
enddo
elseif(paramid==201) then
do j=1,ntot
snow(j)=gfld%fld(j)
lb(j)=gfld%bmap(j)
enddo
else
write(*,*) 'unknow parameter'
stop
endif
! write(*,*) 'height,max,min',height,maxval(var),minval(var)
call gf_free(gfld)
enddo NUM_FIELDS
enddo VERSION ! skgb
endif
CALL BACLOSE(ifile,ierr)
if (allocated(cgrib)) deallocate(cgrib)
nullify(gfld%local)
enddo loopfile
return
end subroutine read_grib2_sngle
subroutine read_grib2_head_time(filename,year,month,day,hour,minute)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_grib2 read grib2 file head time information
! prgmmr: Ming Hu org: GSD date: 2017-05-20
!
! abstract: read grib2 file head
!
!
! program history log:
! 2015-05-20 Hu, initial documentation
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: Zeus
!
!$$$ end documentation block
! use gridmod, only: idsl5,regional
use grib_mod
implicit none
character*100,intent(in) :: filename
integer, intent(out) :: year, month, day, hour, minute
!
integer :: nz
!
type(gribfield) :: gfld
logical :: expand=.true.
integer :: ifile
character(len=1),allocatable,dimension(:) :: cgrib
integer,parameter :: msk1=32000
integer :: lskip, lgrib,iseek
integer :: currlen
integer :: icount , lengrib
integer :: listsec0(3)
integer :: listsec1(13)
integer :: second, fcst
integer :: numfields,numlocal,maxlocal,ierr
integer :: grib_edition
integer :: itot
! real :: dx,dy,lat1,lon1
real :: scale_factor
!
!
integer :: nn,n,j,iret
real :: fldmax,fldmin,sum
!
!
scale_factor=1.0e6
ifile=10
loopfile: do nn=1,1
! write(6,*) 'read in grib2 file head', trim(filename)
lskip=0
lgrib=0
iseek=0
icount=0
itot=0
currlen=0
! Open GRIB2 file
call baopenr(ifile,trim(filename),iret)
if (iret.eq.0) then
VERSION: do
! Search opend file for the next GRIB2 messege (record).
call skgb(ifile,iseek,msk1,lskip,lgrib)
! Check for EOF, or problem
if (lgrib.eq.0) then
exit
endif
! Check size, if needed allocate more memory.
if (lgrib.gt.currlen) then
if (allocated(cgrib)) deallocate(cgrib)
allocate(cgrib(lgrib))
currlen=lgrib
endif
! Read a given number of bytes from unblocked file.
call baread(ifile,lskip,lgrib,lengrib,cgrib)
if(lgrib.ne.lengrib) then
write(*,*) 'ERROR, read_grib2 lgrib ne lengrib', &
lgrib,lengrib
stop 1234
endif
iseek=lskip+lgrib
icount=icount+1
! Unpack GRIB2 field
call gb_info(cgrib,lengrib,listsec0,listsec1, &
numfields,numlocal,maxlocal,ierr)
if(ierr.ne.0) then
write(6,*) 'Error querying GRIB2 message',ierr
stop
endif
itot=itot+numfields
grib_edition=listsec0(2)
if (grib_edition.ne.2) then
exit VERSION
endif
! write(*,*) 'listsec0=',listsec0
! write(*,*) 'listsec1=',listsec1
! write(*,*) 'numfields=',numfields
! get information form grib2 file
n=1
call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
year =gfld%idsect(6) !(FOUR-DIGIT) YEAR OF THE DATA
month =gfld%idsect(7) ! MONTH OF THE DATA
day =gfld%idsect(8) ! DAY OF THE DATA
hour =gfld%idsect(9) ! HOUR OF THE DATA
minute=gfld%idsect(10) ! MINUTE OF THE DATA
second=gfld%idsect(11) ! SECOND OF THE DATA
! write(*,*) 'year,month,day,hour,minute,second='
! write(*,*) year,month,day,hour,minute,second
call gf_free(gfld)
enddo VERSION ! skgb
endif
CALL BACLOSE(ifile,ierr)
nullify(gfld%local)
if (allocated(cgrib)) deallocate(cgrib)
enddo loopfile
return
end subroutine read_grib2_head_time
end module grib2_read_mod