!*****************************************************************************!
program g2print !
! !
use table
use gridinfo
use storage_module
use filelist
use datarray
implicit none
interface
subroutine parse_args(err, a1, h1, i1, l1, a2, h2, i2, l2,&
a3, h3, i3, l3, hlast)
integer :: err
character(len=*) , optional :: a1, a2, a3
character(len=*), optional :: h1, h2, h3
integer , optional :: i1, i2, i3
logical, optional :: l1, l2, l3
character(len=*), optional :: hlast
end subroutine parse_args
end interface
integer :: nunit1 = 12
character(LEN=120) :: gribflnm
integer :: iprint
integer , parameter :: maxlvl = 100
real :: startlat, startlon, deltalat, deltalon
real :: level
character (LEN=9) :: field
character (LEN=3) :: out_format
logical :: readit
integer, dimension(255) :: iuarr = 0
character (LEN=19) :: HSTART, HEND, HDATE
character(LEN=19) :: hsave = '0000-00-00_00:00:00'
integer :: itime
integer :: ntimes
integer :: interval
integer :: ierr
logical :: ordered_by_date
integer :: debug_level
integer :: grib_version
integer :: vtable_columns
character(len=30) :: hopt
logical :: ivb = .FALSE.
logical :: idb = .FALSE.
! -----------------
gribflnm = ' '
call parse_args(ierr, a1='v', l1=ivb, a2='V', l2=idb, hlast=gribflnm)
if (ierr.ne.0) then
call getarg(0, hopt)
write(*,'(//,"Usage: ", A, " [-v] [-V] file",/)') trim(hopt)
write(*,'(" -v : Print more information about the GRIB records")')
write(*,'(" -V : Print way too much information about the GRIB&
& records")')
write(*,'(" file : GRIB file to read"//)')
stop
endif
! -----------------
! Determine GRIB Edition number
grib_version=0
call edition_num(nunit1, trim(gribflnm), grib_version, ierr)
if (ierr.eq.3) STOP 'GRIB file problem'
debug_level = 0
if (ivb) debug_level = 51
if (idb) debug_level = 101
write(6,*) 'reading from grib file = ',gribflnm
LOOP1 : DO
! At the beginning of LOOP1, we are at a new time period.
! Clear the storage arrays and associated level information.
! If we need to read a new grib record, then read one.
if (grib_version.ne.2) then
write(6,*) 'calling r_grib1 with iunit ', nunit1
write(6,*) 'flnm = ',gribflnm
write(6,*) 'GRIB 1 not yet supported in this code'
stop
! Read one record at a time from GRIB1 (and older Editions)
! call r_grib1(nunit1, gribflnm, level, field, &
! hdate, debug_level, ierr, iuarr, iprint)
else
! Read one file of records from GRIB2.
if (debug_level .gt. 100) write(6,*) 'calling r_grib2'
call r_grib2(nunit1, gribflnm, hdate, &
grib_version, debug_level, ierr)
endif
if (ierr.eq.1) then
! We have hit the end of a file. Exit LOOP1.
exit LOOP1
endif
enddo LOOP1
if (grib_version.ne.2) then
call c_close(iuarr(nunit1), iprint, ierr)
iuarr(nunit1) = 0
endif
! And Now we are done:
print*,' '
print*,' '
print*,' Successful completion of g2print '
contains
subroutine sort_filedates
implicit none
integer :: n
logical :: done
if (nfiles > 1) then
done = .FALSE.
do while ( .not. done)
done = .TRUE.
do n = 1, nfiles-1
if (filedates(n) > filedates(n+1)) then
filedates(size(filedates)) = filedates(n)
filedates(n) = filedates(n+1)
filedates(n+1) = filedates(size(filedates))
filedates(size(filedates)) = '0000-00-00 00:00:00.0000'
done = .FALSE.
endif
enddo
enddo
endif
end subroutine sort_filedates
end program g2print
!*****************************************************************************!
SUBROUTINE r_grib2(junit, gribflnm, hdate, &
grib_edition, debug_level, ireaderr)
use grib_mod
use params
use table ! Included to define cg2code
use gridinfo ! Included to define map%
use storage_module ! Included sub put_storage
real, allocatable, dimension(:) :: hold_array
parameter(msk1=32000,msk2=4000)
character(len=1),allocatable,dimension(:) :: cgrib
integer :: listsec0(3)
integer :: listsec1(13)
integer year, month, day, hour, minute, second, fcst
character(len=*) :: gribflnm
character(len=*) :: hdate
character(len=8) :: pabbrev
character(len=20) :: labbrev
character(len=80) :: tabbrev
integer :: lskip, lgrib
integer :: junit, itot, icount, iseek
integer :: grib_edition
integer :: i, j, ireaderr, ith
integer :: currlen
logical :: unpack, expand
type(gribfield) :: gfld
! For subroutine put_storage
real :: level
real :: scale_factor
integer :: iplvl
character (len=9) :: my_field
! For subroutine outout
integer , parameter :: maxlvl = 100
real , dimension(maxlvl) :: plvl
integer :: nlvl
integer , dimension(maxlvl) :: level_array
logical :: verbose=.false.
integer :: debug_level
character(len=4) :: tmp4
character(len=40) :: string
character(len=13) :: pstring = ',t50,":",i14)'
character(len=15) :: rstring = ',t50,":",f14.5)'
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! SET ARGUMENTS
if (debug_level .gt. 50 ) then
unpack=.true.
else
unpack=.false.
endif
expand=.true.
hdate = '0000-00-00_00:00:00'
ierr=0
itot=0
icount=0
iseek=0
lskip=0
lgrib=0
currlen=0
ith=1
scale_factor = 1e6
! do j = 1,10
! write(6,'("j = ",i4," level1 = ",i8," level2 = ",i8)') j, &
! level1(j), level2(j)
! enddo
!/* IOS Return Codes from BACIO: */
!/* 0 All was well */
!/* -1 Tried to open read only _and_ write only */
!/* -2 Tried to read and write in the same call */
!/* -3 Internal failure in name processing */
!/* -4 Failure in opening file */
!/* -5 Tried to read on a write-only file */
!/* -6 Failed in read to find the 'start' location */
!/* -7 Tried to write to a read only file */
!/* -8 Failed in write to find the 'start' location */
!/* -9 Error in close */
!/* -10 Read or wrote fewer data than requested */
!if ireaderr =1 we have hit the end of a file.
!if ireaderr =2 we have hit the end of all the files.
if ( debug_level .gt. 100 ) verbose = .true.
if (verbose) write(6,*) 'begin r_grib2, flnm = ',gribflnm
! Open a byte-addressable file.
CALL BAOPENR(junit,gribflnm,IOS)
if (verbose) write(6,*) 'back from baopenr, ios = ',ios
if (ios.eq.0) then
VERSION: do
! Search opend file for the next GRIB2 messege (record).
if (verbose) write(6,*) 'calling skgb'
call skgb(junit,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),stat=is)
!print *,'G2 allocate(cgrib(lgrib)) status: ',IS
currlen=lgrib
endif
! Read a given number of bytes from unblocked file.
call baread(junit,lskip,lgrib,lengrib,cgrib)
if (lgrib.ne.lengrib) then
print *,'G2 r_grib2: IO Error.',lgrib,".ne.",lengrib
stop 9
endif
iseek=lskip+lgrib
icount=icount+1
if (verbose) PRINT *,'G2 GRIB MESSAGE ',icount,' starts at',lskip+1
! Unpack GRIB2 field
call gb_info(cgrib,lengrib,listsec0,listsec1, &
numfields,numlocal,maxlocal,ierr)
if (ierr.ne.0) then
write(*,*) ' ERROR querying GRIB2 message = ',ierr
stop 10
endif
itot=itot+numfields
grib_edition=listsec0(2)
if (grib_edition.ne.2) then
exit VERSION
endif
! Additional print statments for developer.
if (verbose) then
print *,'G2 SECTION 0: ',(listsec0(j),j=1,3)
print *,'G2 SECTION 1: ',(listsec1(j),j=1,13)
print *,'G2 Contains ',numlocal,' Local Sections ', &
' and ',numfields,' data fields.'
endif
! ----
! Once per file fill in date, model and projection values.
if (lskip.lt.10) then
! Build the 19-character date string, based on GRIB2 header date
! and time information, including forecast time information:
n=1
call gf_getfld(cgrib,lengrib,n,unpack,expand,gfld,ierr)
if (debug_level .gt. 100 ) then
write(6,*) 'gfld%version = ',gfld%version
if (gfld%discipline .eq. 0) then
string = 'Meteorological products'
else if (gfld%discipline .eq. 1) then
string = 'Hydrological products'
else if (gfld%discipline .eq. 2) then
string = 'Land Surface products'
else
string = 'See code table 0.0'
endif
write(6,*) 'Discipline = ',gfld%discipline,' ',string
write(6,*) 'gfld%idsect(1) = ',gfld%idsect(1)
write(6,*) 'gfld%idsect(2) = ',gfld%idsect(2)
write(6,*) 'gfld%idsect(3) = ',gfld%idsect(3)
write(6,*) 'gfld%idsect(4) = ',gfld%idsect(4)
write(6,*) 'gfld%idsect(5) = ',gfld%idsect(5)
write(6,*) 'gfld%idsect(6) = ',gfld%idsect(6)
write(6,*) 'gfld%idsect(7) = ',gfld%idsect(7)
write(6,*) 'gfld%idsect(8) = ',gfld%idsect(8)
write(6,*) 'gfld%idsect(9) = ',gfld%idsect(9)
write(6,*) 'gfld%idsect(10) = ',gfld%idsect(10)
write(6,*) 'gfld%idsect(11) = ',gfld%idsect(11)
write(6,*) 'gfld%idsect(12) = ',gfld%idsect(12)
write(6,*) 'gfld%idsect(13) = ',gfld%idsect(13)
write(6,*) 'gfld%idsectlen = ',gfld%idsectlen
write(6,*) 'gfld%locallen = ',gfld%locallen
write(6,*) 'gfld%ifldnum = ',gfld%ifldnum
write(6,*) 'gfld%ngrdpts = ',gfld%ngrdpts
write(6,*) 'gfld%numoct_opt = ',gfld%numoct_opt
write(6,*) 'gfld%interp_opt = ',gfld%interp_opt
write(6,*) 'gfld%griddef = ',gfld%griddef
if (gfld%igdtnum .eq. 0) then
string = 'Lat/Lon cylindrical equidistant'
else if (gfld%igdtnum .eq. 1) then
string = 'Rotated Lat/Lon'
else if (gfld%igdtnum .eq. 2) then
string = 'Stretched Lat/Lon'
else if (gfld%igdtnum .eq. 20) then
string = 'Polar Stereographic'
else if (gfld%igdtnum .eq. 30) then
string = 'Lambert Conformal'
else if (gfld%igdtnum .eq. 40) then
string = 'Gaussian Lat/Lon'
else if (gfld%igdtnum .eq. 50) then
string = 'Spherical harmonic coefficients'
else
string = 'see code table 3.1'
endif
write(6,*) 'Grid Template number = ',gfld%igdtnum,' ',string
write(6,*) 'gfld%igdtlen = ',gfld%igdtlen
do i = 1, gfld%igdtlen
write(6,*) 'gfld%igdtmpl(',i,') = ',gfld%igdtmpl(i)
enddo
write(6,*) 'gfld%ipdtnum = ',gfld%ipdtnum
write(6,*) 'gfld%ipdtlen = ',gfld%ipdtlen
if ( gfld%ipdtnum .eq. 0 ) then
do i = 1, gfld%ipdtlen
write(6,*) 'gfld%ipdtmpl(',i,') = ',gfld%ipdtmpl(i)
enddo
endif
write(6,*) 'gfld%num_coord = ',gfld%num_coord
write(6,*) 'gfld%ndpts = ',gfld%ndpts
write(6,*) 'gfld%idrtnum = ',gfld%idrtnum
write(6,*) 'gfld%idrtlen = ',gfld%idrtlen
write(6,*) 'gfld%expanded = ',gfld%expanded
write(6,*) 'gfld%ibmap = ',gfld%ibmap
endif
if (debug_level .le. 50) then
write(6,*) '---------------------------------------------------------------------------'
write(6,*) ' rec Prod Cat Param Lvl Lvl Lvl Name Time Fcst'
write(6,*) ' num Disc num code one two hour'
write(6,*) '---------------------------------------------------------------------------'
endif
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
fcst = 0
! Parse the forecast time info from Sect 4.
if (gfld%ipdtnum.eq.0) then ! Product Definition Template 4.0
! Extract forecast time.
fcst = gfld%ipdtmpl(9)
endif
! Compute valid time.
if (verbose) then
print *, 'ymd',gfld%idsect(6),gfld%idsect(7),gfld%idsect(8)
print *, 'hhmm ',gfld%idsect(9),gfld%idsect(10)
endif
call build_hdate(hdate,year,month,day,hour,minute,second)
if (verbose) print *, 'G2 hdate = ',hdate
! call geth_newdate(hdate,hdate,3600*fcst) ! no need for thin in print
! print *, 'G2 hdate (fcst?) = ',hdate
!--
! Indicator of the source (center) of the data.
icenter = gfld%idsect(1)
! Indicator of model (or whatever) which generated the data.
iprocess = gfld%ipdtmpl(5)
if (icenter.eq.7) then
if (iprocess.eq.83 .or. iprocess.eq.84) then
map%source = 'NCEP NAM Model'
elseif (iprocess.eq.81) then
map%source = 'NCEP GFS Model'
elseif (iprocess.eq.96) then
map%source = 'NCEP GFS Model'
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.109) then
map%source = 'NCEP RTMA'
elseif (iprocess.eq.140) then
map%source = 'NCEP NARR'
else
map%source = 'unknown model from NCEP'
write (6,*) 'unknown NCEP model, iprocess = ',iprocess
end if
else if (icenter .eq. 57) then
if (iprocess .eq. 87) then
map%source = 'AFWA AGRMET'
else
map%source = 'AFWA'
endif
else if (icenter .eq. 98) then
map%source = 'ECMWF'
else
map%source = 'unknown model and orig center'
end if
!--
! 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'
if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical Equidistant
map%igrid = 0
map%nx = gfld%igdtmpl(8)
map%ny = gfld%igdtmpl(9)
map%dx = gfld%igdtmpl(17)
map%dy = gfld%igdtmpl(18)
map%lat1 = gfld%igdtmpl(12)
map%lon1 = gfld%igdtmpl(13)
! Scale dx/dy values to degrees, default range is 1e6.
if (map%dx.gt.10000) then
map%dx = map%dx/scale_factor
endif
if (map%dy.gt.10000) then
map%dy = (map%dy/scale_factor)*(-1)
endif
! Scale lat/lon values to 0-180, default range is 1e6.
if (map%lat1.ge.scale_factor) then
map%lat1 = map%lat1/scale_factor
endif
if (map%lon1.ge.scale_factor) then
map%lon1 = map%lon1/scale_factor
endif
! if (iprocess.eq.81 .and. map%dy.eq.0) then !Hack for AFWA.
! map%dx = 0.5
! map%dy = -0.5
! endif
elseif (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid
map%igrid = 3
map%nx = gfld%igdtmpl(8)
map%ny = gfld%igdtmpl(9)
map%lov = gfld%igdtmpl(14)
map%truelat1 = gfld%igdtmpl(19)
map%truelat2 = gfld%igdtmpl(20)
map%dx = gfld%igdtmpl(15)
map%dy = gfld%igdtmpl(16)
map%lat1 = gfld%igdtmpl(10)
map%lon1 = gfld%igdtmpl(11)
elseif(gfld%igdtnum.eq.40) then ! Gaussian Grid (we will call it lat/lon)
map%igrid = 0
map%nx = gfld%igdtmpl(8)
map%ny = gfld%igdtmpl(9)
map%dx = gfld%igdtmpl(17)
map%dy = gfld%igdtmpl(18) ! ?not in Grid Definition Template 3.40 doc
map%lat1 = gfld%igdtmpl(12)
map%lon1 = gfld%igdtmpl(13)
! Scale dx/dy values to degrees, default range is 1e6.
if (map%dx.gt.10000) then
map%dx = map%dx/scale_factor
endif
if (map%dy.gt.10000) then
map%dy = (map%dy/scale_factor)*(-1)
endif
! Scale lat/lon values to 0-180, default range is 1e6.
if (map%lat1.ge.scale_factor) then
map%lat1 = map%lat1/scale_factor
endif
if (map%lon1.ge.scale_factor) then
map%lon1 = map%lon1/scale_factor
endif
print *,'Gaussian Grid: Dx,Dy,lat,lon',map%dx,map%dy, &
map%lat1,map%lon1
elseif (gfld%igdtnum.eq.20) then ! Polar-Stereographic Grid.
map%igrid = 5
map%nx = gfld%igdtmpl(8)
map%ny = gfld%igdtmpl(9)
!map%lov = gfld%igdtmpl(x) ! ?not in Grid Definition Template 3.20 doc
map%truelat1 = 60.
map%truelat2 = 91.
!map%dx = gfld%igdtmpl(x)
!map%dy = gfld%igdtmpl(x)
map%lat1 = gfld%igdtmpl(10)
map%lon1 = gfld%igdtmpl(11)
else
print*, 'GRIB2 Unknown Projection: ',gfld%igdtnum
print*, 'see Code Table 3.1: Grid Definition Template No'
endif
call gf_free(gfld)
endif
! ----
! Continue to unpack GRIB2 field.
if (debug_level .gt. 100) write(6,*) 'numfields = ',numfields
do n=1,numfields ! e.g. U and V would =2, otherwise its usually =1
call gf_getfld(cgrib,lengrib,n,unpack,expand,gfld,ierr)
if (ierr.ne.0) then
write(*,*) ' ERROR extracting field gf_getfld = ',ierr
cycle
endif
! ------------------------------------
! Additional print information for developer.
pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1), &
gfld%ipdtmpl(2))
if (debug_level .gt. 50 ) then
print *
! print *,'G2 FIELD ',n
if (n==1) then
write(*,'(/,"GRIB2 SECTION 0 - INDICATOR SECTION:")')
write(*,'(5x,"Discipline"'//pstring) gfld%discipline
write(*,'(5x,"GRIB Edition Number"'//pstring) gfld%version
write(*,'(5x,"GRIB length"'//pstring) lengrib
write(*,'(/,"GRIB2 SECTION 1 - IDENTIFICATION SECTION:")')
write(*,'(5x,"Length of Section"'//pstring) gfld%idsectlen
write(*,'(5x,"Originating Center ID"'//pstring) &
gfld%idsect(1)
write(*,'(5x,"Subcenter ID"'//pstring) gfld%idsect(2)
write(*,'(5x,"GRIB Master Table Version"'//pstring) &
gfld%idsect(3)
write(*,'(5x,"GRIB Local Table Version"'//pstring) &
gfld%idsect(4)
write(*,'(5x,"Significance of Reference Time"'//pstring) &
gfld%idsect(5)
write(*,'(5x,"Year"'//pstring) gfld%idsect(6)
write(*,'(5x,"Month"'//pstring) gfld%idsect(7)
write(*,'(5x,"Day"'//pstring) gfld%idsect(8)
write(*,'(5x,"Hour"'//pstring) gfld%idsect(9)
write(*,'(5x,"Minute"'//pstring) gfld%idsect(10)
write(*,'(5x,"Second"'//pstring) gfld%idsect(11)
write(*,'(5x,"Production Status of data"'//pstring) &
gfld%idsect(12)
write(*,'(5x,"Type of processed data"'//pstring) &
gfld%idsect(13)
! print *,'G2 SECTION 1: ',(gfld%idsect(j),j=1,gfld%idsectlen)
endif
write(*,'(/,"GRIB2 SECTION 2 - LOCAL SECTION:")')
write(*,'(5x,"Length of Section 2"'//pstring) gfld%locallen
if ( associated(gfld%local).AND.gfld%locallen.gt.0 ) then
do j = 1, gfld%locallen
write(*,'(5x,"Local value "'//pstring) gfld%local(j)
enddo
! print *,'G2 SECTION 2: ',(gfld%local(j),j=1,gfld%locallen)
endif
write(*,'(/,"GRIB2 SECTION 3 - GRID DEFINITION SECTION:")')
! write(*,'(5x,"Length of Section 3"'//pstring) gfld%unknown
write(*,'(5x,"Source of grid definition"'&
//pstring) gfld%griddef
write(*,'(5x,"Number of grid points"'//pstring) gfld%ngrdpts
write(*,'(5x,"Number of octets for addnl points"'//pstring) &
gfld%numoct_opt
write(*,'(5x,"Interpretation list"'//pstring) &
gfld%interp_opt
write(*,'(5x,"Grid Definition Template Number"'//pstring) &
gfld%igdtnum
if (gfld%igdtnum .eq. 0 .or. gfld%igdtnum .eq. 1 .or. &
gfld%igdtnum .eq. 2 .or. gfld%igdtnum .eq. 3 ) then
if (gfld%igdtnum .eq. 0 ) then
write(*,'(5x,"Lat/Lon or Cylindrical Equidistant Grid")')
else if (gfld%igdtnum .eq. 1 ) then
write(*,'(5x,"Rotated Lat/Lon or Cylind. Equi. Grid")')
else if (gfld%igdtnum .eq. 2 ) then
write(*,'(5x,"Stretched Lat/Lon or Cylind. Equi. Grid")')
else if (gfld%igdtnum .eq. 3 ) then
write(*,'(5x,"Stretched and Rotated Lat/Lon Grid")')
endif
write(*,'(10x,"Shape of the Earth"'//pstring) &
gfld%igdtmpl(1)
write(*,'(10x,"Scale factor of spher. Earth"'//pstring) &
gfld%igdtmpl(2)
write(*,'(10x,"Scaled value of spher. Earth"'//pstring) &
gfld%igdtmpl(3)
write(*,'(10x,"Scale factor of major axis"'//pstring) &
gfld%igdtmpl(4)
write(*,'(10x,"Scaled value of major axis"'//pstring) &
gfld%igdtmpl(5)
write(*,'(10x,"Scale factor of minor axis"'//pstring) &
gfld%igdtmpl(6)
write(*,'(10x,"Scaled value of minor axis"'//pstring) &
gfld%igdtmpl(7)
write(*,'(10x,"Ni - points along a parallel"'//pstring) &
gfld%igdtmpl(8)
write(*,'(10x,"Nj - points along a meridian"'//pstring) &
gfld%igdtmpl(9)
write(*,'(10x,"Basic angle of initial domain"'//pstring)&
gfld%igdtmpl(10)
write(*,'(10x,"Subdivisions of basic angle"'//pstring) &
gfld%igdtmpl(11)
write(*,'(10x,"La1"'//pstring) gfld%igdtmpl(12)
write(*,'(10x,"Lo1"'//pstring) gfld%igdtmpl(13)
write(*,'(10x,"Resolution and Component",t50,":",B14.8)')&
gfld%igdtmpl(14)
write(*,'(10x,"La2"'//pstring) gfld%igdtmpl(15)
write(*,'(10x,"Lo2"'//pstring) gfld%igdtmpl(16)
write(*,'(10x,"Di - i-dir increment"'//pstring) &
gfld%igdtmpl(17)
write(*,'(10x,"Dj - j-dir increment"'//pstring) &
gfld%igdtmpl(18)
write(*,'(10x,"Scanning mode"'//pstring) &
gfld%igdtmpl(19)
if (gfld%igdtnum .eq. 1) then
write(*,'(10x,"Lat of southern pole of project"'//pstring)&
gfld%igdtmpl(20)
write(*,'(10x,"Lon of southern pole of project"'//pstring)&
gfld%igdtmpl(21)
write(*,'(10x,"Angle of rotation of projection"'//pstring)&
gfld%igdtmpl(22)
else if (gfld%igdtnum .eq. 2) then
write(*,'(10x,"Lat of the pole of stretching "'//pstring)&
gfld%igdtmpl(20)
write(*,'(10x,"Lon of the pole of stretching "'//pstring)&
gfld%igdtmpl(21)
write(*,'(10x,"Stretching factor"'//pstring) &
gfld%igdtmpl(22)
else if (gfld%igdtnum .eq. 3) then
write(*,'(10x,"Lat of southern pole of project"'//pstring)&
gfld%igdtmpl(20)
write(*,'(10x,"Lon of southern pole of project"'//pstring)&
gfld%igdtmpl(21)
write(*,'(10x,"Angle of rotation of projection"'//pstring)&
gfld%igdtmpl(22)
write(*,'(10x,"Lat of the pole of stretching "'//pstring)&
gfld%igdtmpl(23)
write(*,'(10x,"Lon of the pole of stretching "'//pstring)&
gfld%igdtmpl(24)
write(*,'(10x,"Stretching factor"'//pstring) &
gfld%igdtmpl(25)
endif
else if (gfld%igdtnum .eq. 10) then
write(*,'(5x,"Mercator Grid")')
else if (gfld%igdtnum .eq. 20 .or. gfld%igdtnum .eq. 30) then
if (gfld%igdtnum .eq. 20) then
write(*,'(5x,"Polar Stereographic Grid")')
else if (gfld%igdtnum .eq. 30) then
write(*,'(5x,"Lambert Conformal Grid")')
endif
write(*,'(10x,"Shape of the Earth"'//pstring) &
gfld%igdtmpl(1)
write(*,'(10x,"Scale factor of spher. Earth"'//pstring) &
gfld%igdtmpl(2)
write(*,'(10x,"Scaled value of spher. Earth"'//pstring) &
gfld%igdtmpl(3)
write(*,'(10x,"Scale factor of major axis"'//pstring) &
gfld%igdtmpl(4)
write(*,'(10x,"Scaled value of major axis"'//pstring) &
gfld%igdtmpl(5)
write(*,'(10x,"Scale factor of minor axis"'//pstring) &
gfld%igdtmpl(6)
write(*,'(10x,"Scaled value of minor axis"'//pstring) &
gfld%igdtmpl(7)
write(*,'(10x,"Nx"'//pstring) gfld%igdtmpl(8)
write(*,'(10x,"Ny"'//pstring) gfld%igdtmpl(9)
write(*,'(10x,"La1"'//pstring) gfld%igdtmpl(10)
write(*,'(10x,"Lo1"'//pstring) gfld%igdtmpl(11)
write(*,'(10x,"Resolution and Component",t50,":",B14.8)')&
gfld%igdtmpl(12)
write(*,'(10x,"LaD"'//pstring) gfld%igdtmpl(13)
write(*,'(10x,"LoV"'//pstring) gfld%igdtmpl(14)
write(*,'(10x,"Dx"'//pstring) gfld%igdtmpl(15)
write(*,'(10x,"Dy"'//pstring) gfld%igdtmpl(16)
write(*,'(10x,"Projection Center Flag"'//pstring) &
gfld%igdtmpl(17)
write(*,'(10x,"Scanning mode"'//pstring) &
gfld%igdtmpl(18)
if (gfld%igdtnum .eq. 30) then
write(*,'(10x,"Latin 1 "'//pstring) &
gfld%igdtmpl(19)
write(*,'(10x,"Latin 2 "'//pstring) &
gfld%igdtmpl(20)
write(*,'(10x,"Lat of southern pole of project"'//pstring)&
gfld%igdtmpl(21)
write(*,'(10x,"Lon of southern pole of project"'//pstring)&
gfld%igdtmpl(22)
endif
else if (gfld%igdtnum .eq. 40 .or. gfld%igdtnum .eq. 41) then
if (gfld%igdtnum .eq. 40) then
write(*,'(5x,"Gaussian Lat/Lon Grid")')
else if (gfld%igdtnum .eq. 41) then
write(*,'(5x,"Rotated Gaussian Lat/Lon Grid")')
else if (gfld%igdtnum .eq. 42) then
write(*,'(5x,"Stretched Gaussian Lat/Lon Grid")')
else if (gfld%igdtnum .eq. 43) then
write(*,'(5x,"Stretched and Rotated Gaussian Lat/Lon ")')
endif
else
do j = 1, gfld%igdtlen
write(*,'(5x,"Grid Definition Template entry "'//pstring) &
gfld%igdtmpl(j)
enddo
endif
! print *,'G2 SECTION 3: ',gfld%griddef,gfld%ngrdpts, &
! gfld%numoct_opt,gfld%interp_opt, &
! gfld%igdtnum
! print *,'G2 GRID TEMPLATE 3.',gfld%igdtnum,': ', &
! (gfld%igdtmpl(j),j=1,gfld%igdtlen)
if ( gfld%num_opt .eq. 0 ) then
! print *,'G2 NO Section 3 List Defining No. of Data Points.'
else
print *,'G2 Section 3 Optional List: ', &
(gfld%list_opt(j),j=1,gfld%num_opt)
endif
write(*,'(/,"GRIB2 SECTION 4 - PRODUCT DEFINITION SECTION:")')
! write(*,'(5x,"Length of Section 4"'//pstring) gfld%unknown
write(*,'(5x,"Product Definition Template Number"'//pstring)&
gfld%ipdtnum
do j = 1, gfld%ipdtlen
write(tmp4,'(i4)') j
string = '(5x,"Template Entry '//tmp4 // '"'
write(*,string//pstring) gfld%ipdtmpl(j)
enddo
! print *,'G2 PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ', &
! (gfld%ipdtmpl(j),j=1,gfld%ipdtlen)
!call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev)
!call prvtime(gfld%ipdtnum,gfld%ipdtmpl,listsec1,tabbrev)
write(*,'(5x,"Product Abbreviated Name",t50,":",a14)')&
pabbrev
! print *,'G2 TEXT: ',pabbrev,trim(labbrev)," ",trim(tabbrev)
if ( gfld%num_coord .eq. 0 ) then
! print *,'G2 NO Optional Vertical Coordinate List.'
else
print *,'G2 Section 4 Optional Coordinates: ', &
(gfld%coord_list(j),j=1,gfld%num_coord)
endif
! if ( gfld%ibmap .ne. 255 ) then
! print *,'G2 Num. of Data Points = ',gfld%ndpts, &
! ' with BIT-MAP ',gfld%ibmap
! else
! print *,'G2 Num. of Data Points = ',gfld%ndpts, &
! ' NO BIT-MAP '
! endif
write(*,'(/,"GRIB2 SECTION 5 - DATA REPRESENTATION SECTION:")')
write(*,'(5x,"Data Representation Template Number"'//pstring)&
gfld%idrtnum
do j = 1, gfld%idrtlen
write(tmp4,'(i4)') j
string = '(5x,"Template Entry '//tmp4 // '"'
write(*,string//pstring) gfld%idrtmpl(j)
enddo
! print *,'G2 DRS TEMPLATE 5.',gfld%idrtnum,': ', &
! (gfld%idrtmpl(j),j=1,gfld%idrtlen)
! if ( gfld%ipdtnum .eq. 0 ) then
! if (gfld%ipdtmpl(1) .eq. 0 ) then
! write(6,*) 'Temperature'
! else if (gfld%ipdtmpl(1) .eq. 1 ) then
! write(6,*) 'Moisture'
! else if (gfld%ipdtmpl(1) .eq. 2 ) then
! write(6,*) 'Momentum'
! else if (gfld%ipdtmpl(1) .eq. 3 ) then
! write(6,*) 'Mass'
! endif
! endif
write(*,'(/,"GRIB2 SECTION 6 - BIT-MAP SECTION:")')
write(*,'(5x,"Bit-map indicator"'//pstring) &
gfld%ibmap
fldmax=gfld%fld(1)
fldmin=gfld%fld(1)
sum=gfld%fld(1)
do j=2,gfld%ndpts
if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j)
if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j)
sum=sum+gfld%fld(j)
enddo ! gfld%ndpts
write(*,'(/,"GRIB2 SECTION 7 - DATA SECTION:")')
write(*,'(5x,"Minimum Data Value "'//rstring)&
fldmin
write(*,'(5x,"Maximum Data Value "'//rstring)&
fldmax
! print *,'G2 Data Values:'
! write(*,fmt='("G2 MIN=",f21.8," AVE=",f21.8, &
! " MAX=",f21.8)') fldmin,sum/gfld%ndpts,fldmax
if (debug_level .gt. 100 ) then
print*, 'gfld%fld = ',gfld%fld
! do j=1,gfld%ndpts
! write(*,*) j, gfld%fld(j)
! enddo
endif
endif ! Additional Print information
! ------------------------------------
if ( debug_level .le. 50) then
write(6,987) itot,gfld%discipline,gfld%ipdtmpl(1), &
gfld%ipdtmpl(2),gfld%ipdtmpl(10),gfld%ipdtmpl(12),&
gfld%ipdtmpl(15),pabbrev,hdate,fcst
987 format(2i4,i5,i4,i8,i8,i8,a10,a20,i5.2)
endif
! Deallocate arrays decoding GRIB2 record.
call gf_free(gfld)
enddo ! 1,numfields
enddo VERSION ! skgb
if (debug_level .gt. 50) &
print *, 'G2 total number of fields found = ',itot
CALL BACLOSE(junit,IOS)
ireaderr=1
else
print *,'open status failed because',ios
hdate = '9999-99-99_99:99:99'
ireaderr=2
endif ! ireaderr check
END subroutine r_grib2
!*****************************************************************************!
! Subroutine edition_num !
! !
! Purpose: !
! Read one record from the input GRIB2 file. Based on the information in !
! the GRIB2 header and the user-defined Vtable, decide whether the field in!
! the GRIB2 record is one to process or to skip. If the field is one we !
! want to keep, extract the data from the GRIB2 record, and pass the data !
! back to the calling routine. !
! !
! Argument list: !
! Input: !
! JUNIT : "Unit Number" to open and read from. Not really a Fortran !
! unit number, since we do not do Fortran I/O for the GRIB2 !
! 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. !
! GRIB2FILE: File name to open, if it is not already open. !
! !
! Output: !
! GRIB_EDITION: Set to 1 for GRIB and set to 2 for GRIB2 !
! IERR : Error flag: 0 - no error on read from GRIB2 file. !
! 1 - Hit the end of the GRIB2 file. !
! 2 - The file GRIBFLNM we tried to open does !
! not exist. !
! Author: Paula McCaslin !
! NOAA/FSL !
! Sept 2004 !
!*****************************************************************************!
SUBROUTINE edition_num(junit, gribflnm, grib_edition, ireaderr)
use grib_mod
use params
parameter(msk1=32000,msk2=4000)
character(len=1),allocatable,dimension(:) :: cgrib
integer :: listsec0(3)
integer :: listsec1(13)
character(len=*) :: gribflnm
integer :: lskip, lgrib
integer :: junit
integer :: grib_edition
integer :: i, j, ireaderr
integer :: currlen
character(len=4) :: ctemp
character(len=4),parameter :: grib='GRIB',c7777='7777'
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! SET ARGUMENTS
itot=0
icount=0
iseek=0
lskip=0
lgrib=0
currlen=0
!/* IOS Return Codes from BACIO: */
!/* 0 All was well */
!/* -1 Tried to open read only _and_ write only */
!/* -2 Tried to read and write in the same call */
!/* -3 Internal failure in name processing */
!/* -4 Failure in opening file */
!/* -5 Tried to read on a write-only file */
!/* -6 Failed in read to find the 'start' location */
!/* -7 Tried to write to a read only file */
!/* -8 Failed in write to find the 'start' location */
!/* -9 Error in close */
!/* -10 Read or wrote fewer data than requested */
!if ireaderr =1 we have hit the end of a file.
!if ireaderr =2 we have hit the end of all the files.
!if ireaderr =3 beginning characters 'GRIB' not found
! write(6,*) 'junit = ',junit,' gribflnm = ',gribflnm
! Open a byte-addressable file.
CALL BAOPENR(junit,gribflnm,IOS) ! from w3lib
! write(6,*) 'ios = ',ios
if (ios.eq.0) then
! Search opend file for the next GRIB2 messege (record).
call skgb(junit,iseek,msk1,lskip,lgrib)
! Check for EOF, or problem
if (lgrib.eq.0) then
STOP "Grib2 file or date problem, stopping in edition_num."
endif
! Check size, if needed allocate more memory.
if (lgrib.gt.currlen) then
if (allocated(cgrib)) deallocate(cgrib)
allocate(cgrib(lgrib),stat=is)
currlen=lgrib
endif
! Read a given number of bytes from unblocked file.
call baread(junit,lskip,lgrib,lengrib,cgrib)
! Check for beginning of GRIB message in the first 100 bytes
istart=0
do j=1,100
ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
if (ctemp.eq.grib ) then
istart=j
exit
endif
enddo
if (istart.eq.0) then
ireaderr=3
print*, "The beginning 4 characters >GRIB< were not found."
endif
! Unpack Section 0 - Indicator Section to extract GRIB edition field
iofst=8*(istart+5)
call gbyte(cgrib,discipline,iofst,8) ! Discipline
iofst=iofst+8
call gbyte(cgrib,grib_edition,iofst,8) ! GRIB edition number
print *, 'ungrib - grib edition num', grib_edition
CALL BACLOSE(junit,IOS)
ireaderr=1
else if (ios .eq. -4) then
print *,'edition_num: unable to open ',gribflnm
stop 'edition_num'
else
print *,'edition_num: open status failed because',ios,gribflnm
ireaderr=2
endif ! ireaderr check
END subroutine edition_num
subroutine parse_args(err, a1, h1, i1, l1, a2, h2, i2, l2, a3, h3, i3, l3, &
hlast)
integer :: err
character(len=*) , optional :: a1, a2, a3
character(len=*), optional :: h1, h2, h3
integer , optional :: i1, i2, i3
logical, optional :: l1, l2, l3
character(len=*), optional :: hlast
character(len=100) :: hold
integer :: ioff = 0
if (present(hlast)) then
ioff = -1
endif
err = 0
narg = iargc()
numarg = narg + ioff
i = 1
LOOP : do while ( i <= numarg)
ierr = 1
if (present(i1)) then
call checkiarg(i, a1, i1, ierr)
elseif (present(h1)) then
call checkharg(i, a1, h1, ierr)
elseif (present(l1)) then
call checklarg(i, a1, l1, ierr)
endif
if (ierr.eq.0) cycle LOOP
if (present(i2)) then
call checkiarg(i, a2, i2, ierr)
elseif (present(h2)) then
call checkharg(i, a2, h2, ierr)
elseif (present(l2)) then
call checklarg(i, a2, l2, ierr)
endif
if (ierr.eq.0) cycle LOOP
if (present(i3)) then
call checkiarg(i, a3, i3, ierr)
elseif (present(h3)) then
call checkharg(i, a3, h3, ierr)
elseif (present(l3)) then
call checklarg(i, a3, l3, ierr)
endif
if (ierr.eq.0) cycle LOOP
err = 1
call getarg(1, hold)
write(*, '("arg = ", A)') trim(hold)
exit LOOP
enddo LOOP
if (present(hlast)) then
if (narg.eq.0) then
err = 1
else
call getarg(narg, hlast)
endif
endif
contains
subroutine checkiarg(c, a, i, ierr)
integer :: c
character(len=*) :: a
integer :: i
character(len=100) :: hold
ierr = 1
call getarg(c, hold)
if ('-'//a.eq.trim(hold)) then
c = c + 1
call getarg(c, hold)
read(hold, *) i
c = c + 1
ierr = 0
elseif ('-'//a .eq. hold(1:len_trim(a)+1)) then
hold = hold(len_trim(a)+2: len(hold))
read(hold, *) i
c = c + 1
ierr = 0
endif
end subroutine checkiarg
subroutine checkharg(c, a, h, ierr)
integer :: c
character(len=*) :: a
character(len=*) :: h
character(len=100) :: hold
ierr = 1
call getarg(c, hold)
if ('-'//a.eq.trim(hold)) then
c = c + 1
call getarg(c, hold)
h = trim(hold)
c = c + 1
ierr = 0
elseif ('-'//a .eq. hold(1:len_trim(a)+1)) then
hold = hold(len_trim(a)+2: len(hold))
h = trim(hold)
c = c + 1
ierr = 0
endif
end subroutine checkharg
subroutine checklarg(c, a, l, ierr)
integer :: c
character(len=*) :: a
logical :: l
character(len=100) :: hold
ierr = 1
call getarg(c, hold)
if ('-'//a.eq.trim(hold)) then
l = .TRUE.
c = c + 1
ierr = 0
endif
end subroutine checklarg
end subroutine parse_args