!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RDNMCGRB ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rdnmcgrb(nx_ext,ny_ext,nz_ext, & 12,75
gribfile,grbflen, gribtime, &
gridesc, iproj_grb, gthin, &
ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit, &
pi_grb,pj_grb,ipole, di_grb,dj_grb, &
latsw,lonsw, latne,lonne, &
latrot,lonrot,angrot, &
latstr,lonstr,facstr, &
lattru1,lattru2,lontrue, &
scanmode, iscan,jscan,kscan, &
ires,iearth,icomp, &
jpenta,kpenta,mpenta,ispect,icoeff, &
xp_grb,yp_grb, xo_grb,yo_grb,zo_grb, &
rcdata,var_grb2d,var_grb3d,var_lev3d,lvldbg, iret)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! This subroutine is to read and select given variables with GRIB
! grid ID and level type from NMC GRIB files.
!
! The decoder of GRIB is from NMC, ftp:nic.fb4.noaa.gov/pub
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 01/02/1996 Initialized.
!
! MODIFICATIONS:
!
! 09/24/1998 (Dennis A. Moon, Chief Scientist, SSESCO)
! Added sorting in vertical levels to assure that the level changes
! monotonically decreasing.
!
! 1999/08/03 (Gene Bassett)
! Excluded soil levels (grid types 111 & 112) from vertical level sorting.
!
! 2000/01/05 (Eric Kemp)
! Fixed problem with century change over (ipds(21) runs from 1 to 100).
!
! 2002/09/04 (Kevin W. Thomas)
! ETA soil temp and moisture data were never sorted, as they had flags
! indicating "sigma vertical coordinates". When CONDUIT (LDM) data is
! the source, the data order isn't guaranteed, so it needs to be sorted.
!
! Duplicated data reports that sometimes occur in CONDUIT data cause fatal
! errors. Duplicates are ignored.
!
! 2004/08/23 (Kevin W. Thomas)
! Variable "grbunit" is now an allocatable character string. It is
! guaranteed to be the right size to hold a C file pointer. This will
! end the "segmentation violation" problems.
!
! 2006/10/19 (Kevin W. Thomas)
! Avoid skipping record if PDS(16)=1. This shows up in all 0 hour
! records when NCEP moved the data from one system to another.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! gribfile GRIB file name
! grbflen Length of GRIB file name
! gribtime Time of GRIB data
!
! OUTPUT:
!
! iproj_grb Map projection number of GRIB data
! trlon_grb True longitude of GRIB data (degrees E)
! latnot_grb(2) True latitude(s) of GRIB data (degrees N)
! swlat_grb Latitude of first grid point at southwest conner
! swlon_grb Longitude of first grid point at southwest conner
!
! dx_grb x-direction grid length
! dy_grb y-direction grid length
!
! iret Return flag
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'gribcst.inc' ! General GRIB definitions
INCLUDE 'mp.inc'
INTEGER :: nx_ext,ny_ext,nz_ext
INTEGER :: grbflen, grbflen_new
CHARACTER (LEN=*) :: gribfile
CHARACTER (LEN=*) :: gribtime
CHARACTER (LEN=256) :: gribfile_new
CHARACTER (LEN=1) :: opnmod
INTEGER, INTENT(IN) :: lvldbg
!
!-----------------------------------------------------------------------
!
! GRIB grid information
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=42) :: gridesc ! Grid description
INTEGER :: gdsflg ! GDS indicator
INTEGER :: iproj_grb ! Map projection indicator
INTEGER :: gthin ! Indicator of whether the grid is "thinned"
INTEGER :: ni_grb ! Number of points along x-axis
INTEGER :: nj_grb ! Number of points along y-axis
INTEGER :: np_grb ! Total number of horizontal grid points
INTEGER :: nk_grb ! Number of vertical parameters
REAL :: zk_grb(nz_ext) ! Vertical coordinate parameters
INTEGER :: npeq ! Number of lat circles from pole to equator
INTEGER :: nit(nz_ext) ! Number of x-points for thinned grid
REAL :: pi_grb ! x-coordinate of pole point
REAL :: pj_grb ! y-coordinate of pole point
INTEGER :: ipole ! Projection center flag
REAL :: di_grb ! x-direction increment or grid length
REAL :: dj_grb ! y-direction increment or grid length
REAL :: latsw ! Latitude of South West corner point
REAL :: lonsw ! Longitude of South West corner point
REAL :: latne ! Latitude of North East corner point
REAL :: lonne ! Longitude of North East corner point
REAL :: lattru1 ! Latitude (1st) at which projection is true
REAL :: lattru2 ! Latitude (2nd) at which projection is true
REAL :: lontrue ! Longitude at which projection is true
REAL :: latrot ! Latitude of southern pole of rotation
REAL :: lonrot ! Longitude of southern pole of rotation
REAL :: angrot ! Angle of rotation
REAL :: latstr ! Latitude of the pole of stretching
REAL :: lonstr ! Longitude of the pole of stretching
REAL :: facstr ! Stretching factor
INTEGER :: scanmode ! Scanning indicator
INTEGER :: iscan ! x-direction scanning indicator
INTEGER :: jscan ! y-direction scanning indicator
INTEGER :: kscan ! FORTRAN index scanning indicator
INTEGER :: ires ! Resolution direction increments indicator
INTEGER :: iearth ! Earth shape indicator: spherical or oblate?
INTEGER :: icomp ! (u,v) components decomposition indicator
INTEGER :: jpenta ! J-Pentagonal resolution parameter
INTEGER :: kpenta ! K-Pentagonal resolution parameter
INTEGER :: mpenta ! M-Pentagonal resolution parameter
INTEGER :: ispect ! Spectral representation type
INTEGER :: icoeff ! Spectral coefficient storage mode
REAL :: xp_grb ! X coordinate of sub-satellite point
REAL :: yp_grb ! Y coordinate of sub-satellite point
REAL :: xo_grb ! X coordinate of image sector origin
REAL :: yo_grb ! Y coordinate of image sector origin
REAL :: zo_grb ! Camera altitude from center of Earth
INTEGER :: iret ! Return flag
!
!-----------------------------------------------------------------------
!
! Temporary arrays to read GRIB file
!
!-----------------------------------------------------------------------
!
INTEGER :: nrecs ! number of records
INTEGER :: fsize ! Size of GRIB file
INTEGER :: stdio_size ! Size of "struct file" (C code in gribio_c.c).
INTEGER :: istatus ! Return status
CHARACTER(len=1), allocatable :: grbunit(:) ! C structure, size stdio_size
!wdt RLC
LOGICAL (KIND=1) :: ibms(nx_ext*ny_ext) ! BMS logical array for data bit map
REAL :: rcdata(nx_ext*ny_ext) ! temporary data array
REAL :: var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt) ! GRIB 2-D variables
REAL :: var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt) ! GRIB 3-D variables
INTEGER :: var_lev3d(nz_ext,n3dvs,n3dlvt)
! Levels (hybrid) for each 3-D variable
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,l,m,n, ij, ki, kj, ilev
INTEGER :: iyr,imo,iday,ihr,fhr, myr
INTEGER :: ileft, iright, iincr
INTEGER :: jleft, jright, jincr
INTEGER :: chklev, mlvtp
INTEGER :: iendn,itypec,wdlen, nwrd
INTEGER :: iret2 ! another return flag
LOGICAL :: grdflg ! flag for grid information
LOGICAL :: fexist
LOGICAL :: verbose = .FALSE.
CHARACTER(LEN=11), PARAMETER :: cmd = 'rdnmcgrib: '
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF (lvldbg > 90) verbose = .TRUE.
!
!-----------------------------------------------------------------------
!
! Find the type of endian, type of character set and length of
! machine word
!
! iendn - Integer for big-endian or little-endian
! = 0 big-endian
! = 1 little-endian
! = 2 cannot compute
! itypec - Integer for type of character set
! = 0 ascii character set
! = 1 ebcdic character set
! = 2 not ascii or ebcdic
! wdlen - Integer for words size of computer in bytes
! = 4 for 32 bit (4 bytes) computers
! = 8 for 64 bit (8 bytes) computers
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) THEN
CALL findvartype
(iendn,itypec,wdlen)
IF ( iendn == 2 ) THEN
WRITE(6,'(a)') 'Unknown endian type. Stopped in rdnmcgrb'
CALL arpsstop("endian error",1)
END IF
IF ( itypec == 2 ) THEN
WRITE(6,'(a)') 'Unknown character set. Stopped in rdnmcgrb'
CALL arpsstop("endian error",1)
END IF
!
!-----------------------------------------------------------------------
!
! initialize the var_lev3d array so we can sort the records by
! level values as we read them. We can't assume monotonic ordering
! of the levels
!
!-----------------------------------------------------------------------
!
DO k=1, nz_ext
DO n=1, n3dvs
DO m=1, n3dlvt
var_lev3d(k,n,m)= -999999
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Open the GRIB file
!
!-----------------------------------------------------------------------
!
READ (gribtime,'(i4,3i2,1x,i3)') iyr,imo,iday,ihr,fhr
myr = MOD( iyr, 100 )
IF (myr == 0) THEN
myr = 100
END IF
INQUIRE(FILE=gribfile(1:grbflen), EXIST = fexist )
IF( fexist ) THEN
gribfile_new = gribfile
grbflen_new = grbflen
GO TO 200
ENDIF
INQUIRE(FILE=trim(gribfile(1:grbflen))//'.Z', EXIST = fexist )
IF( fexist ) THEN
CALL unixcmd( 'cp '//trim(gribfile(1:grbflen))//'.Z tem_ext_file.Z' )
INQUIRE(FILE='tem_ext_file', EXIST = fexist )
IF( fexist ) call unixcmd( '/bin/rm tem_ext_file' )
CALL uncmprs( 'tem_ext_file.Z' )
gribfile_new = 'tem_ext_file'
grbflen_new = 12
GO TO 200
END IF
INQUIRE(FILE=trim(gribfile(1:grbflen))//'.gz', EXIST = fexist )
IF( fexist ) THEN
CALL unixcmd( 'cp '//trim(gribfile(1:grbflen))//'.gz tem_ext_file.gz' )
INQUIRE(FILE='tem_ext_file', EXIST = fexist )
IF( fexist ) call unixcmd( '/bin/rm tem_ext_file' )
CALL uncmprs( 'tem_ext_file.gz' )
gribfile_new = 'tem_ext_file'
grbflen_new = 12
GO TO 200
END IF
WRITE(6,'(/1x,a,/1x,a/)')
WRITE(6,'(3(1x,a))') 'GRIB file '//gribfile(1:grbflen)// &
' or its compressed version not found.', &
'Please check option extdfmt and its instruction.', &
'Program warning in RDNMCGRB.'
IRET = 1
! Note other errors go to 999 instead of 998. We can't go to 999, as
! grbunit hasn't been allocated yet.
GOTO 998
200 CONTINUE
WRITE(6,'(8A)') cmd, 'Opening GRIB ', gribfile(1:grbflen)
opnmod = 'r'
CALL gsize ( stdio_size )
ALLOCATE(grbunit(stdio_size), STAT = istatus)
CALL check_alloc_status(istatus, "grbunit")
CALL gopen ( grbunit, gribfile_new, grbflen_new, opnmod, iret )
IF ( iret /= 0 ) THEN
WRITE(6,'(a,a,/,a,i4)') &
'Error in opening file ',gribfile_new, &
'C function gopen return status: ',iret
! DEALLOCATE(grbunit)
! RETURN
IRET = 1
GOTO 999
END IF
!
!-----------------------------------------------------------------------
!
! Get the size of the GRIB file
!
!-----------------------------------------------------------------------
!
fsize = 0
CALL gseek( grbunit, fsize, 2, iret )
WRITE (6,'(2A,I10,A)') cmd, 'Size of GRIB file is', fsize,' bytes'
!
!-----------------------------------------------------------------------
!
! Scan the GRIB file to count records
!
!-----------------------------------------------------------------------
!
nrecs = 0
CALL grbscan( grbunit, fsize, nprods, rcstart, rcbytes, nrecs )
IF ( nrecs > nprods ) THEN
WRITE (6,'(a/a,i6,a,i6/a/a)') &
'The actual number of products exceeded the expected number.', &
'Number of assumed products: ', nprods, &
'Number of actual products: ', nrecs, &
'Program stopped in RDNMCGRB. Please change NPRODS in ', &
'file gribcst.inc and re-run the program.'
CALL gclose( grbunit, iret )
! STOP
CALL arpsstop("GRIB file problem",1)
ELSE IF ( nrecs <= 0 ) THEN
WRITE (6,'(a,a,a)') &
'No GRIB message found in file ',gribfile_new, &
'Program stopped in RDNMCGRB.'
CALL gclose( grbunit, iret )
! STOP
CALL arpsstop("GRIB file problem",1)
END IF
!
!-----------------------------------------------------------------------
!
! Read the GRIB file and fill in the variable arrays
!
!-----------------------------------------------------------------------
!
IF( n2dvars > n2dvs .OR. n2dlvtps > n2dlvt ) THEN
WRITE(6,'(a,i3,a,i3,a)') &
'Array var_nr2d(',n2dvs,',',n2dlvt,') declared in gribcst.inc too small.'
WRITE(6,'(a,i3,a,i3,a)') &
'It needs to be at least ',n2dvars,'x',n2dlvtps,'.'
write(6,'(a)') 'Program stopped in RDNMCGRB.'
! STOP
CALL arpsstop("GRIB file problem",1)
ENDIF
DO n=1,n2dvars
DO m=1,n2dlvtps
var_nr2d(n,m) = 0
END DO
END DO
IF( n3dvars > n3dvs .OR. n3dlvtps > n3dlvt ) THEN
write(6,'(a,i3,a,i3,a)') &
'Array var_nr3d(',n3dvs,',',n3dlvt,') declared in gribcst.inc too small.'
write(6,'(a,i3,a,i3,a)') &
'It needs to be at least ',n3dvars,'x',n3dlvtps,'.'
write(6,'(a)') 'Program stopped in RDNMCGRB.'
! STOP
CALL arpsstop("GRIB file problem",1)
ENDIF
DO m=1,n3dlvtps
DO n=1,n3dvars
var_nr3d(n,m) = 0
END DO
END DO
grdflg = .false.
DO l=1,nrecs
DO i=1,ipdsz
ipds(i) = 0
END DO
DO i=1,igdsz
igds(i) = 0
END DO
CALL gseek( grbunit, rcstart(l), 0, iret )
CALL gread( grbunit, mgrib, rcbytes(l), iret )
IF ( iendn == 1 ) THEN ! little endian machine
nwrd = (rcbytes(l)+3)/4
CALL swap4byte( mgrib, nwrd )
END IF
CALL w3fi63( mgrib, ipds, igds, ibms, rcdata, mptrs, iret )
!print *,cmd, 'PDS#',L,ipds(3),ipds(5),ipds(6),ipds(7),rcdata(1)
IF ( ipds(3) /= gridtyp ) THEN
IF (verbose) &
WRITE (6,'(/A,I5,2(A,I4))') 'Ignoring GRIB record', l, &
': Grid ID is ipds(3) =',ipds(3), ', expected gridtyp=',gridtyp
CYCLE
END IF
mlvtp = 0
chklev = 0
DO m=1,n3dlvtps
IF ( ipds(6) == levtyp3d(m) ) THEN
chklev = 1
mlvtp = m
GO TO 30
END IF
END DO
DO m=1,n2dlvtps
IF ( ipds(6) == levtyp2d(m) ) THEN
IF (levtyp2d(m) == 105) THEN ! T_2m, QV_2m, U_10m, V_10m
SELECT CASE (ipds(5))
CASE (11, 51, 52) ! Must at 2m
IF (ipds(7) /= 2) EXIT
CASE (33, 34) ! Must at 10m
IF (ipds(7) /= 10) EXIT
END SELECT
END IF
chklev = 2
mlvtp = m
EXIT
END IF
END DO
30 CONTINUE
IF ( chklev == 0 ) THEN
IF (verbose) &
WRITE (6,'(/A,I5,A,I4)') 'Ignoring GRIB record', l, &
': Type of level not expected, ipds(6) =',ipds(6)
CYCLE
END IF
IF ( ipds(16) /= 0 .AND. ipds(16) /= 10 .AND. &
ipds(16) /= 1 .AND. ipds(16) /= 4 ) THEN
IF (verbose) &
WRITE(6,'(/A,I5,A,I5)') 'Ignoring GRIB record', l, &
': average, or time difference - ',ipds(16)
CYCLE
ELSE IF ( ipds(8) /= myr .OR. ipds(9) /= imo .OR. &
ipds(10) /= iday .OR. ipds(11) /= ihr .OR. &
ipds(14) /= fhr ) THEN
IF (ipds(16) == 4 .AND. ipds(15) == fhr) THEN
! good here for accumulated rain
ELSE
IF (verbose) &
WRITE (6,'(/A,I5,A,5I2.2,"f",I3.3,3A)') &
'Ignoring GRIB record', l, ': time of record (', &
ipds(21)-1,(ipds(j),j=8,11),ipds(14), &
') does not match expected (', gribtime, ")"
CYCLE
END IF
END IF
gdsflg = IAND(ipds(4),128)/128
IF ( .NOT. grdflg ) THEN
CALL grbgrid(gridtyp, gdsflg, igds, &
gridesc, iproj_grb, gthin, &
ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit, &
pi_grb,pj_grb,ipole, di_grb,dj_grb, &
latsw,lonsw, latne,lonne, &
latrot,lonrot,angrot, &
latstr,lonstr,facstr, &
lattru1,lattru2,lontrue, &
scanmode, iscan,jscan,kscan, &
ires,iearth,icomp, &
jpenta,kpenta,mpenta,ispect,icoeff, &
xp_grb,yp_grb, xo_grb,yo_grb,zo_grb)
WRITE (6,'(2A)') cmd, gridesc
grdflg = .true.
END IF
IF ( igds(1) /= mproj_grb ) THEN
IF (verbose) &
WRITE(6,'(/A,I5,2(A,I4))') 'Ignoring GRIB record', l, &
': Map projection is igds(1) =', igds(1), &
', expected mproj_grb =', mproj_grb
CYCLE
END IF
IF ( igds(2) /= nx_ext .OR. igds(3) /= ny_ext) THEN
WRITE (6,'(a/a,i4,2x,a,i4/a,i4,2x,a,i4/a)') &
'Horizontal dimension was not consistent,', &
' igds(2) = ', igds(2), ' igds(3) = ', igds(3), &
' nx_ext = ', nx_ext , ' ny_ext = ', ny_ext, &
'Program stopped in RDNMCGRB.'
! STOP
CALL arpsstop("GRIB data problem",1)
END IF
IF ( iscan == 0 ) THEN
ileft = 1
iright = nx_ext
iincr = 1
ELSE
ileft = nx_ext
iright = 1
iincr = -1
END IF
IF ( jscan /= 0 ) THEN
jleft = 1
jright = ny_ext
jincr = 1
ELSE
jleft = ny_ext
jright = 1
jincr = -1
END IF
IF ( kscan == 0 ) THEN
ki = 1
kj = nx_ext
ELSE IF ( kscan == 1 ) THEN
ki = ny_ext
kj = 1
ELSE
IF (verbose) &
WRITE(6,'(/A,I5,2(A,I4))') 'Ignoring GRIB record', l, &
'Unknown scanning mode, kscan = ', kscan
CYCLE
END IF
IF (verbose) &
WRITE (6,'(/A,I5,3(A,I5),$)') 'GRIB record', l, &
': param =', ipds(5), ', level type =', ipds(6), ', level =', ipds(7)
IF ( chklev == 1 ) THEN ! 3-D variable
DO n=1,n3dvars
! DO m=1,n3dlvtps ! we already knew which level type it is,
! then why do this loop?
IF ( ipds(5) == var_id3d(n,mlvtp) ) THEN
IF (verbose) WRITE (6,'(A)') ' accept 3D.'
var_nr3d(n,mlvtp) = var_nr3d(n,mlvtp) + 1
IF (levtyp3d(mlvtp) == 111 .OR. levtyp3d(mlvtp) == 112) THEN
!
!-----------------------------------------------------------------------
!
! Insert the level and data values into the var_lev3d and var_grb3d
! in the order that they occur in the data file (e.g. below surface
! coordinates).
!
! This data must be sorted. When input is CONDUIT (LDM), data isn't
! guaranteed to be in the proper order in the file!
!
!-----------------------------------------------------------------------
!
! ilev = var_nr3d(n,mlvtp)
ilev= 0
DO k=1,nz_ext
IF ( ipds(7) == var_lev3d(k,n,mlvtp) ) THEN
WRITE(6,*) 'Duplicate report, non-pressure coordinates'
var_nr3d(n,mlvtp) = var_nr3d(n,mlvtp) - 1
GOTO 555
ENDIF
IF ( ipds(7) < var_lev3d(k,n,mlvtp) .OR. &
var_lev3d(k,n,mlvtp) < 0.0) THEN
ilev = k
EXIT
END IF
END DO
IF ( ilev == 0 ) THEN
PRINT *, 'couldnt locate level ',ipds(7)
PRINT *, 'for var ID#',var_id3d(n,mlvtp)
PRINT *, ' n = ', n, ' m = ', mlvtp
PRINT *, ' levtyp3d = ',levtyp3d(mlvtp)
! STOP
CALL arpsstop("GRIB data problem",1)
END IF
ELSE
!
!-----------------------------------------------------------------------
!
! Insert the level and data values into the var_lev3d and var_grb3d
! arrays in order of decreasing level value,
!
! For example if the level type is pressure levels the levels
! will be sorted in order of decreasing pressure
!
! ilev is the level at which to insert the current data
!
!-----------------------------------------------------------------------
!
ilev= 0
DO k=1,nz_ext
IF ( ipds(7) == var_lev3d(k,n,mlvtp) ) then
WRITE(6,*) 'Duplicate report, pressure coordinates'
var_nr3d(n,mlvtp) = var_nr3d(n,mlvtp) - 1
GOTO 555
ENDIF
IF ( ipds(7) > var_lev3d(k,n,mlvtp) ) THEN
ilev = k
EXIT
END IF
END DO
IF ( ilev == 0 ) THEN
PRINT *, 'couldnt locate level ',ipds(7)
PRINT *, 'for var ID#',var_id3d(n,mlvtp)
PRINT *, ' n = ', n, ' m = ', mlvtp
! STOP
CALL arpsstop("GRIB data problem",1)
END IF
ENDIF
! slide all the following levels down one to make room
DO k=nz_ext,ilev+1, -1
var_lev3d(k,n,mlvtp)= var_lev3d(k-1,n,mlvtp)
DO j=jleft,jright,jincr
DO i=ileft,iright,iincr
var_grb3d(i,j,k,n,mlvtp)= var_grb3d(i,j,k-1,n,mlvtp)
END DO
END DO
END DO
!
! now insert the data at the appropriate level
!
var_lev3d(ilev,n,mlvtp) = ipds(7)
DO j=jleft,jright,jincr
DO i=ileft,iright,iincr
ij = ki*iincr*(i-ileft) + kj*jincr*(j-jleft) + 1
var_grb3d(i,j,ilev,n,mlvtp) = rcdata(ij)
END DO
END DO
END IF
555 CONTINUE
! END DO
END DO
ELSE IF ( chklev == 2 ) THEN ! 2-D variable
IF ( ipds(16) == 0 .OR. ipds(16) == 10 .OR. &
ipds(16) == 1 .OR. ipds(16) == 4 ) THEN
! Don't use average,
! or time difference fields.
DO n=1,n2dvars
! DO m=1,n2dlvtps
IF ( ipds(5) == var_id2d(n,mlvtp) ) THEN
IF (verbose) WRITE (6,'(a)') ' accept 2D.'
var_nr2d(n,mlvtp) = var_nr2d(n,mlvtp) + 1
DO i=ileft,iright,iincr
DO j=jleft,jright,jincr
ij = ki*iincr*(i-ileft) + kj*jincr*(j-jleft) + 1
var_grb2d(i,j,n,mlvtp) = rcdata(ij)
END DO
END DO
END IF
! END DO
END DO
END IF
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Set good status
!
!-----------------------------------------------------------------------
!
iret = 0
999 CONTINUE
! Don't trash out "iret".
CALL gclose( grbunit, iret2 )
DEALLOCATE( grbunit )
998 CONTINUE
ENDIF ! IF (myproc == 0)
CALL mpupdatei(iproj_grb,1)
CALL mpupdatei(gthin,1)
CALL mpupdatei(ni_grb,1)
CALL mpupdatei(nj_grb,1)
CALL mpupdatei(np_grb,1)
CALL mpupdatei(nk_grb,1)
CALL mpupdater(zk_grb,nz_ext)
CALL mpupdatei(npeq,1)
CALL mpupdater(nit,nz_ext)
CALL mpupdater(pi_grb,1)
CALL mpupdater(pj_grb,1)
CALL mpupdatei(ipole,1)
CALL mpupdater(di_grb,1)
CALL mpupdater(dj_grb,1)
CALL mpupdater(latsw,1)
CALL mpupdater(lonsw,1)
CALL mpupdater(latne,1)
CALL mpupdater(lonne,1)
CALL mpupdater(latrot,1)
CALL mpupdater(lonrot,1)
CALL mpupdater(angrot,1)
CALL mpupdater(latstr,1)
CALL mpupdater(lonstr,1)
CALL mpupdater(facstr,1)
CALL mpupdater(lattru1,1)
CALL mpupdater(lattru2,1)
CALL mpupdater(lontrue,1)
CALL mpupdatei(scanmode,1)
CALL mpupdatei(iscan,1)
CALL mpupdatei(jscan,1)
CALL mpupdatei(kscan,1)
CALL mpupdatei(ires,1)
CALL mpupdatei(iearth,1)
CALL mpupdatei(icomp,1)
CALL mpupdatei(jpenta,1)
CALL mpupdatei(kpenta,1)
CALL mpupdatei(mpenta,1)
CALL mpupdatei(ispect,1)
CALL mpupdatei(icoeff,1)
CALL mpupdater(xp_grb,1)
CALL mpupdater(yp_grb,1)
CALL mpupdater(xo_grb,1)
CALL mpupdater(yo_grb,1)
CALL mpupdater(zo_grb,1)
CALL mpupdater(rcdata,nx_ext*ny_ext)
CALL mpupdater(var_grb2d,nx_ext*ny_ext*n2dvs*n2dlvt)
CALL mpupdater(var_grb3d,nx_ext*ny_ext*nz_ext*n3dvs*n3dlvt)
CALL mpupdatei(var_lev3d,nz_ext*n3dvs*n3dlvt)
CALL mpupdatei(iret,1)
CALL mpupdatei(n2dvars,1)
CALL mpupdatei(n2dlvtps,1)
CALL mpupdatei(var_nr2d,n2dvs*n2dlvt)
CALL mpupdatei(n3dvars,1)
CALL mpupdatei(var_nr3d,n3dvs*n3dlvt)
RETURN
END SUBROUTINE rdnmcgrb
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GRBSCAN ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE grbscan(grbunit,fsize,nprods, & 3
rcstart,rcbytes,nrecs)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! This subroutine is to scan the NMC GRIB files.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 01/03/1996 Initialized.
!
! MODIFICATIONS:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! grbunit Unit of GRIB file opened to scan for GRIB records
! nprods Maxmum record number
! fsize Size of the GRIB file
!
! OUTPUT:
!
! rcstart Starting position (bytes) for each record
! rcbytes Record length in bytes
! nrecs Record number actually scanned
! iret Return flag
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: grbunit
INTEGER :: fsize
INTEGER :: nprods
INTEGER :: rcbytes(nprods) ! record length in bytes
INTEGER :: rcstart(nprods) ! record starting byte in a GRIB file
INTEGER :: nrecs
INTEGER :: iret
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: nbytes
CHARACTER (LEN=1) :: buffer(8)
!
!-----------------------------------------------------------------------
!
! External function
!
!-----------------------------------------------------------------------
!
INTEGER :: char2i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
nbytes = 0
nrecs = 0
100 IF ( nbytes < fsize ) THEN
CALL gseek( grbunit,nbytes,0,iret )
IF ( iret /= 0 ) GO TO 200
CALL gread( grbunit,buffer,8,iret )
IF ( iret /= 8 ) GO TO 200
IF ( char2i(buffer(1)) /= 71 .OR. & ! test ASCII 'G'
char2i(buffer(2)) /= 82 .OR. & ! test ASCII 'R'
char2i(buffer(3)) /= 73 .OR. & ! test ASCII 'I'
char2i(buffer(4)) /= 66 .OR. & ! test ASCII 'B'
char2i(buffer(8)) /= 1 ) THEN ! test Edition #, 1
nbytes = nbytes + 1
GO TO 100
END IF
nrecs = nrecs + 1
IF ( nrecs <= nprods ) THEN
rcstart(nrecs) = nbytes
rcbytes(nrecs) = char2i(buffer(5))*65536 &
+ char2i(buffer(6))*256 &
+ char2i(buffer(7))
nbytes = nbytes + rcbytes(nrecs)
END IF
GO TO 100
END IF
200 RETURN
END SUBROUTINE grbscan
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE rdgrbdims ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rdgrbdims(gribfile,grbflen,ni_grb,nj_grb,iret) 1,7
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! This subroutine is to read dimension size from NMC GRIB #1 file.
!
! The decoder of GRIB is from NMC, ftp:nic.fb4.noaa.gov/pub
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yunheng Wang
! 08/12/2005 Initialized.
!
! MODIFICATIONS:
! 08/17/2007 Yunheng Wang
! Interface changed.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! gribfile GRIB file name
! grbflen Length of GRIB file name
!
! OUTPUT:
!
! ni_grb Grid size in west-east direction
! nj_grb Grid size in south-north direction
!
! iret Return flag
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
CHARACTER(LEN=*), INTENT(IN) :: gribfile
INTEGER, INTENT(IN) :: grbflen
INTEGER, INTENT(OUT) :: ni_grb ! Number of points along x-axis
INTEGER, INTENT(OUT) :: nj_grb ! Number of points along y-axis
INTEGER, INTENT(OUT) :: iret ! Return flag
!
!-----------------------------------------------------------------------
!
! Temporary arrays to read GRIB file
!
!-----------------------------------------------------------------------
!
INTEGER :: nrecs ! number of records
INTEGER :: iproj ! Map projection indicator
INTEGER :: fsize ! Size of GRIB file
INTEGER :: stdio_size ! Size of "struct file" (C code in gribio_c.c).
CHARACTER(LEN=1), ALLOCATABLE :: grbunit(:) ! C structure, size stdio_size
LOGICAL (KIND=1), ALLOCATABLE :: ibms(:) ! BMS logical array for data bit map
REAL, ALLOCATABLE :: rcdata(:) ! temporary data array
!
!-----------------------------------------------------------------------
!
! Misc. Local variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=1) :: opnmod
INTEGER :: iendn,itypec,wdlen, nwrd
INTEGER :: l, i
LOGICAL, PARAMETER :: verbose = .FALSE.
CHARACTER(LEN=11), PARAMETER :: cmd = 'rdgrbdims: '
INCLUDE 'gribcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ALLOCATE(ibms(nbufsz/4), STAT = iret)
CALL check_alloc_status(iret, 'rdgrbdims:ibms')
ALLOCATE(rcdata(nbufsz/4), STAT = iret)
CALL check_alloc_status(iret, 'rdgrbdims:rcdata')
!
!-----------------------------------------------------------------------
!
! Find the type of endian, type of character set and length of
! machine word
!
! iendn - Integer for big-endian or little-endian
! = 0 big-endian
! = 1 little-endian
! = 2 cannot compute
! itypec - Integer for type of character set
! = 0 ascii character set
! = 1 ebcdic character set
! = 2 not ascii or ebcdic
! wdlen - Integer for words size of computer in bytes
! = 4 for 32 bit (4 bytes) computers
! = 8 for 64 bit (8 bytes) computers
!
!-----------------------------------------------------------------------
!
CALL findvartype(iendn,itypec,wdlen)
IF ( iendn == 2 ) THEN
WRITE(6,'(a)') 'Unknown endian type. Stopped in rdnmcgrb'
iret = -1
GOTO 999
END IF
IF ( itypec == 2 ) THEN
WRITE(6,'(a)') 'Unknown character set. Stopped in rdnmcgrb'
iret = -1
GOTO 999
END IF
!
!-----------------------------------------------------------------------
!
! Open the GRIB file
!
!-----------------------------------------------------------------------
!
WRITE(6,'(8A)') cmd, 'Opening GRIB ', TRIM(gribfile)
opnmod = 'r'
CALL gsize ( stdio_size )
ALLOCATE(grbunit(stdio_size), STAT = iret)
CALL check_alloc_status(iret, "rdgrbdim:grbunit")
CALL gopen ( grbunit, gribfile, grbflen, opnmod, iret )
IF ( iret /= 0 ) THEN
WRITE(6,'(a,a,/,a,i4)') 'Error in opening file ',gribfile, &
'C function gopen return status: ',iret
DEALLOCATE(grbunit)
GOTO 999
END IF
!
!-----------------------------------------------------------------------
!
! Get the size of the GRIB file
!
!-----------------------------------------------------------------------
!
fsize = 0
CALL gseek( grbunit, fsize, 2, iret )
WRITE (6,'(2A,I10,A)') cmd, 'Size of GRIB file is', fsize,' bytes'
!
!-----------------------------------------------------------------------
!
! Scan the GRIB file to count records
!
!-----------------------------------------------------------------------
!
nrecs = 0
CALL grbscan( grbunit, fsize, nprods, rcstart, rcbytes, nrecs )
IF ( nrecs > nprods ) THEN
WRITE (6,'(a/a,i6,a,i6/a/a)') &
'The actual number of products exceeded the expected number.', &
'Number of assumed products: ', nprods, &
'Number of actual products: ', nrecs, &
'Program stopped in RDNMCGRB. Please change NPRODS in ', &
'file gribcst.inc and re-run the program.'
CALL gclose( grbunit, iret )
iret = -3
GOTO 999
ELSE IF ( nrecs <= 0 ) THEN
WRITE (6,'(a,a,a)') &
'No GRIB message found in file ',gribfile, &
'Program stopped in RDNMCGRB.'
CALL gclose( grbunit, iret )
iret = -3
GOTO 999
END IF
!
!-----------------------------------------------------------------------
!
! Read the first record in GRIB file
!
!-----------------------------------------------------------------------
!
DO l = 1,nrecs
DO i=1,ipdsz
ipds(i) = 0
END DO
DO i=1,igdsz
igds(i) = 0
END DO
CALL gseek( grbunit, rcstart(l), 0, iret )
CALL gread( grbunit, mgrib, rcbytes(l), iret )
IF ( iendn == 1 ) THEN ! little endian machine
nwrd = (rcbytes(l)+3)/4
CALL swap4byte( mgrib, nwrd )
END IF
CALL w3fi63( mgrib, ipds, igds, ibms, rcdata, mptrs, iret )
iproj = IABS(igds(1))
IF (iproj == 0 .OR. iproj == 3) THEN
ni_grb = igds(2)
nj_grb = igds(3)
ELSE
WRITE(6,'()') 'Unknow map projection: ',iproj
iret = -5
EXIT
END IF
iret = 0 ! set good status
EXIT
END DO
!
!-----------------------------------------------------------------------
!
! close grib file
!
!-----------------------------------------------------------------------
!
CALL gclose( grbunit, iret )
DEALLOCATE( grbunit )
999 CONTINUE
DEALLOCATE(ibms,rcdata)
RETURN
END SUBROUTINE rdgrbdims