!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GRIBENC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE gribenc(itype,wdlen,grbpkbit,npts,fvar,ivar, & 3,5
ipds,igds,ibdshd,pds,gds, &
nbufsz,bds,ibds,msglen,mgrib, &
tem1,item1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Make a GRIB message for variable fvar or ivar for given PDS and
! GDS.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 11/05/1995
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! itype Data type: 0 for floating array and 1 for integer array
! npts Number of grid points
! fvar Floating array to be written into GRIB file
! ivar Integer array to be written into GRIB file
! ipds Integer array of PDS
! igds Integer array of GDS
! pds PDS (GRIB Section 1).
! gds GDS (GRIB Section 3).
!
! nbusz Buffer size of a GRIB message array
!
! OUTPUT:
!
! bds BDS except the first 11 octets header
! mgrib Buffer carrying the GRIB message
! msglen Length of the GRIB message
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: itype ! Data type: 0 - floating, 1 - integer
INTEGER :: wdlen ! Length of machine word in bytes
INTEGER :: grbpkbit ! Bit length in packing GRIB data
INTEGER :: npts ! Number of grid points
REAL :: fvar(npts) ! Floating array
INTEGER :: ivar(npts) ! Integer array
INTEGER :: ipds(25)
INTEGER :: igds(25)
INTEGER :: ibdshd(4)
CHARACTER (LEN=1) :: pds(28) ! PDS
CHARACTER (LEN=1) :: gds(42) ! GDS
CHARACTER (LEN=1) :: bdshd(11) ! BDS header
INTEGER :: nbufsz ! Size of GRIB message array
INTEGER :: ibds(nbufsz/4) ! Identical to BDS
CHARACTER (LEN=1) :: bds(nbufsz) ! BDS
INTEGER :: msglen ! Length of GRIB message
CHARACTER (LEN=1) :: mgrib(nbufsz) ! GRIB message array
REAL :: tem1(npts) ! working array
INTEGER :: item1(npts) ! working array
!
!-----------------------------------------------------------------------
!
! Local GRIB parameters
!
!-----------------------------------------------------------------------
!
INTEGER :: gribid(4) ! GRIB indicator 'GRIB'
DATA gribid/ 71,82,73,66/ ! ASCII code for 'G', 'R', 'I', 'B'
INTEGER :: gribend(4)
DATA gribend/ 55,55,55,55/ ! ASCII code for '7', '7', '7', '7'
INTEGER :: bdslen
INTEGER :: pdslen
INTEGER :: gdslen
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i
INTEGER :: npos,npts1
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Making PDS - GRIB Section 1
!
!-----------------------------------------------------------------------
!
CALL mkpds
(ipds,pds)
!
!-----------------------------------------------------------------------
!
! Making GDS - GRIB Section 2
!
!-----------------------------------------------------------------------
!
CALL mkgds(igds,gds,npts1)
IF ( npts /= npts1 ) THEN
WRITE (6,'(a/a,i8,a,i8/a)') &
'The point number was not correct', &
'npts = ',npts1, 'while nx*ny = ',npts, &
'Model stopped in GRIBENC'
CALL arpsstop(' ',1)
END IF
pdslen = ICHAR(pds(1))*65536 + ICHAR(pds(2))*256 + ICHAR(pds(3))
gdslen = ICHAR(gds(1))*65536 + ICHAR(gds(2))*256 + ICHAR(gds(3))
!
!-----------------------------------------------------------------------
!
! Making BDS - GRIB Section 4
!
!-----------------------------------------------------------------------
!
CALL mkbds(itype,wdlen,grbpkbit,npts,fvar,ivar, &
pds,ibdshd,bdshd, &
nbufsz,ibds,bdslen)
msglen = 8 + pdslen + gdslen + bdslen + 11 + 4 ! no bit map section
! write (6,'(a,i3,a,i3,a,i8,a,i8)')
! : ' pdslen = ', pdslen, ' gdslen = ', gdslen,
! : ' bdslen = ', bdslen+11, ' msglen = ', msglen
!
!-----------------------------------------------------------------------
!
! Output sections to the buffer.
!
! Move SECTION 0 (8 BYTES) into mgrib
!
!-----------------------------------------------------------------------
!
npos = 0
DO i = 1,4
mgrib(i) = CHAR(gribid(i))
END DO
mgrib(5) = CHAR(MOD(msglen/65536,256)) ! Total length of
mgrib(6) = CHAR(MOD(msglen/256, 256)) ! the GRIB message
mgrib(7) = CHAR(MOD(msglen, 256)) ! in three octets
mgrib(8) = CHAR(01) ! Edition 1
npos = npos + 8
!
!-----------------------------------------------------------------------
!
! Move SECTION 1 - 'PDS' into mgrib
!
!-----------------------------------------------------------------------
!
IF ( pdslen > 0) THEN
DO i=1,pdslen
mgrib(npos+i) = pds(i)
END DO
ELSE
WRITE (6,'(a)') 'PDS length <= 0, pdslen = ', pdslen
END IF
npos = npos + pdslen
!
!-----------------------------------------------------------------------
!
! Move SECTION 2 - 'GDS' into mgrib
!
!-----------------------------------------------------------------------
!
IF ( gdslen > 0) THEN
DO i=1,gdslen
mgrib(npos+i) = gds(i)
END DO
ELSE
WRITE (6,'(a)') 'GDS length <= 0, gdslen = ', gdslen
END IF
npos = npos + gdslen
!
!-----------------------------------------------------------------------
!
! No SECTION 3 (BMS) to move
!
!-----------------------------------------------------------------------
!
! IF ( bmslen.gt.0) THEN
! DO 230 i=1,bmslen
! mgrib(npos+i) = bms(i)
!230 CONTINUE
! END IF
!
! npos = npos + bmslen
!
!-----------------------------------------------------------------------
!
! Move SECTION 4 - 'BDS' into mgrib
!
!-----------------------------------------------------------------------
!
DO i=1,11
mgrib(npos+i) = bdshd(i)
END DO
npos = npos + 11
IF ( bdslen > 0) THEN
DO i=1,bdslen
mgrib(npos+i) = bds(i)
END DO
END IF
npos = npos + bdslen
!
!-----------------------------------------------------------------------
!
! End with SECTION 5 ('7777' where '7's ASCII code is 55)
!
!-----------------------------------------------------------------------
!
DO i=1,4
mgrib(npos+i) = CHAR(gribend(i))
END DO
npos = npos + 4
IF ( npos /= msglen ) THEN
WRITE (6,'(a/a,i8,a,i8/a)') &
'The GRIB message length was incorrect.', &
'npos = ',npos, ' msglen = ',msglen, &
'Program stopped in GRIBENC'
CALL arpsstop(' ',1)
END IF
RETURN
END SUBROUTINE gribenc
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MKPDS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE mkpds(ipds,pds) 1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Make PDS from ipds
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 11/05/1995
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! ipds Integer array of PDS
!
! OUTPUT:
!
! pds PDS (GRIB Section 1).
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: ipds(25)
CHARACTER (LEN=1) :: pds(28) ! PDS
INTEGER :: k,levflg,level,dscl
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
pds(1) = CHAR(MOD(ipds(1)/65536,256))
pds(2) = CHAR(MOD(ipds(1)/256,256))
pds(3) = CHAR(MOD(ipds(1),256))
pds(4) = CHAR(ipds(2))
pds(5) = CHAR(ipds(3))
pds(6) = CHAR(ipds(4))
pds(7) = CHAR(ipds(5))
pds(8) = CHAR(ior(ishft(ipds(6),7),ishft(ipds(7),6)))
pds(9) = CHAR(ipds(8))
pds(10) = CHAR(ipds(9))
levflg = ipds(9)
!
!-----------------------------------------------------------------------
!
! check if level is in two octets or one
!
!-----------------------------------------------------------------------
!
IF ( (levflg >= 1.AND.levflg <= 100).OR.levflg == 102.OR. &
levflg == 103.OR.levflg == 105.OR.levflg == 107.OR. &
levflg == 109.OR.levflg == 111.OR.levflg == 113.OR. &
levflg == 115.OR.levflg == 125.OR.levflg == 160.OR. &
levflg == 200.OR.levflg == 201 ) THEN
level = ipds(11) ! one level stored in pds(11) & pds(12)
IF ( level < 0 ) THEN
level = - level ! positive integer for the abs(level)
level = ior(level,32768) ! the first bit of 16 bits for sign
END IF
pds(11) = CHAR(MOD(level/256,256))
pds(12) = CHAR(MOD(level,256))
ELSE
pds(11) = CHAR(ipds(10)) ! level 1 stored in pds(11)
pds(12) = CHAR(ipds(11)) ! level 2 stored in pds(12)
END IF
pds(13) = CHAR(ipds(12))
pds(14) = CHAR(ipds(13))
pds(15) = CHAR(ipds(14))
pds(16) = CHAR(ipds(15))
pds(17) = CHAR(ipds(16))
pds(18) = CHAR(ipds(17))
!
!-----------------------------------------------------------------------
!
! Check if time P1 is in two octets or one
!
!-----------------------------------------------------------------------
!
IF (ipds(20) == 10) THEN
pds(19) = CHAR(MOD(ipds(18)/256,256))
pds(20) = CHAR(MOD(ipds(18),256))
ELSE
pds(19) = CHAR(ipds(18))
pds(20) = CHAR(ipds(19))
END IF
pds(21) = CHAR(ipds(20))
pds(22) = CHAR(MOD(ipds(21)/256,256))
pds(23) = CHAR(MOD(ipds(21),256))
pds(24) = CHAR(ipds(22))
pds(25) = CHAR(ipds(23))
pds(26) = CHAR(ipds(24))
dscl = ipds(25)
IF (dscl < 0) THEN
dscl = -dscl ! If D-scale less than 0,
dscl = ior(dscl,32768) ! set the highest bit of two bytes to 1
END IF
pds(27) = CHAR(MOD(dscl/256,256))
pds(28) = CHAR(MOD(dscl, 256))
!
!-----------------------------------------------------------------------
!
! Set the rest bytes of PDS to zero
!
!-----------------------------------------------------------------------
!
IF (ipds(1) > 28) THEN
k = ipds(1)
DO i = 29,k
pds(i) = CHAR(0)
END DO
ELSE IF ( ipds(1) < 28 ) THEN
WRITE (6,'(a,i3/a)') &
'The PDS length must be no less than 28 octets, ipds(1) = ', &
ipds(1), 'Program stopped in MKPDS'
END IF
RETURN
END SUBROUTINE mkpds
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MKGDS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE mkgds(igds,gds,npts) 1,1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Make GDS from igds
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 11/05/1995
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! igds Integer array of GDS
!
! OUTPUT:
!
! gds GDS (GRIB Section 2)
! npts Number of grid points
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: igds(25)
CHARACTER (LEN=1) :: gds(42) ! GDS
INTEGER :: npts
!
!-----------------------------------------------------------------------
!
! Local GRIB parameters
!
!-----------------------------------------------------------------------
!
INTEGER :: gdslen
INTEGER :: lat,lon
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! ARPS only support 3 types of projection
!
!-----------------------------------------------------------------------
!
IF ( igds(3) == 5 ) THEN
gdslen = 32
ELSE IF (igds(3) == 1 .OR. igds(3) == 3 ) THEN
gdslen = 42
ELSE ! grid type not valid
WRITE (6,'(a)') &
'ARPS does not support this type of projection grid: ', &
igds(3), 'Program stopped in MKGDS'
CALL arpsstop(' ',1)
END IF
!
!-----------------------------------------------------------------------
!
! Number of grid points
!
!-----------------------------------------------------------------------
!
npts = igds(4)*igds(5) ! npts = nx * ny
!
!-----------------------------------------------------------------------
!
! Put GDS length in bytes 1,2,3
!
!-----------------------------------------------------------------------
!
gds(1) = CHAR(MOD(gdslen/65536,256))
gds(2) = CHAR(MOD(gdslen/ 256,256))
gds(3) = CHAR(MOD(gdslen ,256))
gds(4) = CHAR(igds(1)) ! nz use one byte
gds(5) = CHAR(igds(2)) ! = 255, N/A for PV & PL
gds(6) = CHAR(igds(3)) ! data representation
gds(17) = CHAR(igds(8)) ! resolution and u&v component flag
!
!-----------------------------------------------------------------------
!
! Mercator projection
!
!-----------------------------------------------------------------------
!
IF ( igds(3) == 1 ) THEN
gds(7) = CHAR(MOD(igds(4)/256,256)) ! nx uses two bytes
gds(8) = CHAR(MOD(igds(4), 256))
gds(9) = CHAR(MOD(igds(5)/256,256)) ! ny uses two bytes
gds(10) = CHAR(MOD(igds(5), 256))
lat = igds(6) ! swlat
IF ( lat < 0 ) THEN
lat = -lat
lat = ior(lat,8388608) ! minus sign put in the highest bit
END IF
gds(11) = CHAR(MOD(lat/65536,256))
gds(12) = CHAR(MOD(lat/256, 256))
gds(13) = CHAR(MOD(lat, 256))
lon = igds(7) ! swlon
IF ( lon < 0 ) THEN
lon = -lon
lon = ior(lon,8388608) ! minus sign put in the highest bit
END IF
gds(14) = CHAR(MOD(lon/65536,256))
gds(15) = CHAR(MOD(lon/256, 256))
gds(16) = CHAR(MOD(lon, 256))
lat = igds(9) ! latitude increment
IF ( lat < 0 ) THEN
lat = -lat
lat = ior(lat,8388608) ! minus sign put in the highest bit
END IF
gds(18) = CHAR(MOD(lat/65536,256))
gds(19) = CHAR(MOD(lat/ 256,256))
gds(20) = CHAR(MOD(lat ,256))
lon = igds(10) ! longitude increment
IF ( lon < 0 ) THEN
lon = -lon
lon = ior(lon,8388608) ! minus sign put in the highest bit
END IF
gds(21) = CHAR(MOD(lon/65536,256))
gds(22) = CHAR(MOD(lon/256, 256))
gds(23) = CHAR(MOD(lon, 256))
gds(24) = CHAR(MOD(igds(13)/65536,256)) ! true latitude
gds(25) = CHAR(MOD(igds(13)/256, 256))
gds(26) = CHAR(MOD(igds(13), 256))
gds(27) = CHAR(00)
gds(28) = CHAR(igds(14))
gds(29) = CHAR(MOD(igds(12)/65536,256))
gds(30) = CHAR(MOD(igds(12)/256, 256))
gds(31) = CHAR(MOD(igds(12), 256))
gds(32) = CHAR(MOD(igds(11)/65536,256))
gds(33) = CHAR(MOD(igds(11)/256, 256))
gds(34) = CHAR(MOD(igds(11), 256))
gds(35) = CHAR(00)
gds(36) = CHAR(00)
gds(37) = CHAR(00)
gds(38) = CHAR(00)
gds(39) = CHAR(00)
gds(40) = CHAR(00)
gds(41) = CHAR(00)
gds(42) = CHAR(00)
!
!-----------------------------------------------------------------------
!
! Lambert Conformal projection
!
!-----------------------------------------------------------------------
!
ELSE IF ( igds(3) == 3 ) THEN
gds( 7) = CHAR(MOD(igds(4)/256,256)) ! nx
gds( 8) = CHAR(MOD(igds(4), 256))
gds( 9) = CHAR(MOD(igds(5)/256,256)) ! ny
gds(10) = CHAR(MOD(igds(5), 256))
lat = igds(6)
IF ( lat < 0 ) THEN
lat = -lat
lat = ior(lat,8388608)
END IF
gds(11) = CHAR(MOD(lat/65536,256))
gds(12) = CHAR(MOD(lat/256, 256))
gds(13) = CHAR(MOD(lat, 256))
lon = igds(7)
IF ( lon < 0 ) THEN
lon = -lon
lon = ior(lon,8388608)
END IF
gds(14) = CHAR(MOD(lon/65536,256))
gds(15) = CHAR(MOD(lon/256, 256))
gds(16) = CHAR(MOD(lon, 256))
lon = igds(9)
IF ( lon < 0 ) THEN
lon = -lon
lon = ior(lon,8388608)
END IF
gds(18) = CHAR(MOD(lon/65536,256))
gds(19) = CHAR(MOD(lon/256, 256))
gds(20) = CHAR(MOD(lon, 256))
gds(21) = CHAR(MOD(igds(10)/65536,256))
gds(22) = CHAR(MOD(igds(10)/ 256,256))
gds(23) = CHAR(MOD(igds(10), 256))
gds(24) = CHAR(MOD(igds(11)/65536,256))
gds(25) = CHAR(MOD(igds(11)/256, 256))
gds(26) = CHAR(MOD(igds(11), 256))
gds(27) = CHAR(igds(12))
gds(28) = CHAR(igds(13))
gds(29) = CHAR(MOD(igds(15)/65536,256))
gds(30) = CHAR(MOD(igds(15)/256, 256))
gds(31) = CHAR(MOD(igds(15), 256))
gds(32) = CHAR(MOD(igds(16)/65536,256))
gds(33) = CHAR(MOD(igds(16)/256, 256))
gds(34) = CHAR(MOD(igds(16), 256))
lat = igds(17)
IF ( lat < 0 ) THEN
lat = -lat
lat = ior(lat,8388608)
END IF
gds(35) = CHAR(MOD(lat/65536,256))
gds(36) = CHAR(MOD(lat/256, 256))
gds(37) = CHAR(MOD(lat, 256))
lon = igds(18)
IF ( lon < 0 ) THEN
lon = -lon
lon = ior(lon,8388608)
END IF
gds(38) = CHAR(MOD(lon/65536,256))
gds(39) = CHAR(MOD(lon/256, 256))
gds(40) = CHAR(MOD(lon, 256))
gds(41) = CHAR(00)
gds(42) = CHAR(00)
!
!-----------------------------------------------------------------------
!
! Polar stereographic projection
!
!-----------------------------------------------------------------------
!
ELSE IF ( igds(3) == 5 ) THEN
gds( 7) = CHAR(MOD(igds(4)/256,256))
gds( 8) = CHAR(MOD(igds(4), 256))
gds( 9) = CHAR(MOD(igds(5)/256,256))
gds(10) = CHAR(MOD(igds(5), 256))
lat = igds(6)
IF ( lat < 0 ) THEN
lat = -lat
lat = ior(lat,8388608)
END IF
gds(11) = CHAR(MOD(lat/65536,256))
gds(12) = CHAR(MOD(lat/256, 256))
gds(13) = CHAR(MOD(lat, 256))
lon = igds(7)
IF ( lon < 0 ) THEN
lon = -lon
lon = ior(lon,8388608)
END IF
gds(14) = CHAR(MOD(lon/65536,256))
gds(15) = CHAR(MOD(lon/256, 256))
gds(16) = CHAR(MOD(lon, 256))
lon = igds(9)
IF ( lon < 0 ) THEN
lon = -lon
lon = ior(lon,8388608)
END IF
gds(18) = CHAR(MOD(lon/65536,256))
gds(19) = CHAR(MOD(lon/256, 256))
gds(20) = CHAR(MOD(lon, 256))
gds(21) = CHAR(MOD(igds(10)/65536,256))
gds(22) = CHAR(MOD(igds(10)/256, 256))
gds(23) = CHAR(MOD(igds(10), 256))
gds(24) = CHAR(MOD(igds(11)/65536,256))
gds(25) = CHAR(MOD(igds(11)/256, 256))
gds(26) = CHAR(MOD(igds(11), 256))
gds(27) = CHAR(igds(12))
gds(28) = CHAR(igds(13))
gds(29) = CHAR(00)
gds(30) = CHAR(00)
gds(31) = CHAR(00)
gds(32) = CHAR(00)
END IF
RETURN
END SUBROUTINE mkgds
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MKBDS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE mkbds(itype,wdlen,grbpkbit,npts,fvar,ivar, & 1,3
pds,ibdshd,bdshd, &
nbufsz,ibds,bdslen)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Pack the data and make the GRIB BDS Section.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 11/05/1995
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! itype Data type: 0 for floating array and 1 for integer array
! wdlen Length of a computer word in bits
! grbpkbit DATA packing length in bits
!
! npts Number of grid points
! fvar Floating array to be written into GRIB file
! ivar Integer array to be written into GRIB file
!
! ibdshd Integer array for BDS header
!
! pds PDS (GRIB Section 1).
!
! nbfusz Buffer size of a GRIB message array
!
! OUTPUT:
!
! ibds Identical to BDS
! bdshd BDS first 11 octets header
! bdslen Length of bds
!
! TEMPORARY:
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: itype ! Data type: 0 - floating, 1 - integer
INTEGER :: grbpkbit ! DATA packing length in bits
INTEGER :: wdlen ! Length of a computer word in bits
INTEGER :: npts ! Array size
REAL :: fvar(npts) ! Floating array
INTEGER :: ivar(npts) ! Integer array
CHARACTER (LEN=1) :: pds(28) ! PDS
INTEGER :: ibdshd(4) ! Integer array for BDS header
CHARACTER (LEN=1) :: bdshd(11) ! BDS header
INTEGER :: nbufsz ! Size of GRIB message array
INTEGER :: ibds(nbufsz/wdlen) ! Identical to BDS
INTEGER :: bdslen ! Length of BDS
!
!-----------------------------------------------------------------------
!
! Local GRIB parameters
!
!-----------------------------------------------------------------------
!
INTEGER :: dscl, bscl
REAL :: scale, fmax,fmin
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: bwidth, nbits, iskip1, iskip2
INTEGER :: i, j
REAL :: tema
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
iskip1 = 0
iskip2 = 0
ibdshd(3) = itype ! set the floating or integer data flag
!
!-----------------------------------------------------------------------
!
! Perform decimal scale to either floating or integet data array
!
!-----------------------------------------------------------------------
!
dscl = ICHAR(pds(27))*256 + ICHAR(pds(28))
IF ( IAND(dscl,32768) /= 0 ) THEN
dscl = - IAND(dscl,32761)
END IF
scale = 10.0**dscl
IF ( itype == 0 ) THEN
DO i = 1,npts
fvar(i) = fvar(i)*scale
END DO
ELSE IF ( itype == 1 ) THEN
DO i = 1,npts
fvar(i) = FLOAT(ivar(i))*scale
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Find the minimum and maximum value of the array and subtract the
! minimum
!
!-----------------------------------------------------------------------
!
CALL a3dmax0(fvar,1,npts,1,npts,1,1,1,1,1,1,1,1, fmax,fmin)
!wdt 2001/05/24 Eliminate constant field hack. RLC
!!tmp arpscntl - try having no contsant fields to fix ncl:
! IF (fmax.eq.fmin) THEN
! fmax = fmin + 1.0e-10
! fvar(1) = fmax
! ENDIF
DO i=1,npts
fvar(i) = fvar(i) - fmin
END DO
IF ( (fmax-fmin) == 0.0 ) THEN
bwidth = 0
bscl = 0
scale = 1.0
ELSE IF ( itype /= 0 ) THEN
bwidth = 8
bscl = 0
scale = 1.0
ELSE
bwidth = grbpkbit ! use the value specified by user
tema = (fmax-fmin)/(2.**(bwidth+1)-1.)
IF ( tema /= 0.0 ) THEN
tema = LOG(tema)/LOG(2.0) + 2
END IF
bscl = MIN( INT(tema), INT(tema+SIGN(1.0,tema)) )
scale = 2.0**bscl
END IF
!
!-----------------------------------------------------------------------
!
! Do binary scale to the data
!
!-----------------------------------------------------------------------
!
DO i=1,npts
ivar(i) = nint( fvar(i)/scale )
END DO
!
!-----------------------------------------------------------------------
!
! Calculate the length of BDS which must be an even number of bytes
!
!-----------------------------------------------------------------------
!
nbits = npts*bwidth + 11*8
j = MOD( nbits, 16 )
IF ( j /= 0 ) THEN
nbits = nbits + 16 -j
END IF
bdslen = nbits/8
!
!-----------------------------------------------------------------------
!
! Make the header of BDS
!
!-----------------------------------------------------------------------
!
bdshd(1) = CHAR(MOD(bdslen/65536,256))
bdshd(2) = CHAR(MOD(bdslen/256, 256))
bdshd(3) = CHAR(MOD(bdslen, 256))
bdshd(4) = CHAR( ibdshd(1)*128 + ibdshd(2)*64 &
+ ibdshd(3)*32 + ibdshd(4)*16 )
IF ( bscl < 0 ) THEN
bscl = 32768 - bscl
END IF
bdshd(5) = CHAR(MOD(bscl/256,256))
bdshd(6) = CHAR(MOD(bscl, 256))
!
!-----------------------------------------------------------------------
!
! Convert the reference floating value to IBM GRIB representation
!
!-----------------------------------------------------------------------
!
CALL flt2ibm( fmin, bdshd(7), wdlen*8 )
! write (6,'(a,i6,a,e16.10,a,i3,a,i4,a,i8,a,i8)')
! :' E = ',bscl, ' R = ',fmin, ' L = ',bwidth,
! :' A = ',iexp, ' B = ',imant, ' M = ',nint((fmax-fmin)/scale)
bdshd(11) = CHAR(MOD(bwidth,256))
IF ( bwidth /= 0 ) THEN
CALL grbsbytes( ibds,ivar,iskip1,bwidth,iskip2,npts )
ELSE
ibds(1) = 0
END IF
bdslen = bdslen - 11
RETURN
END SUBROUTINE mkbds
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FLT2IBM ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE flt2ibm(rfloat,chribm,kbits) 18
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Convert floating point number from machine
! representation to GRIB representation.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: John Hennessy, ECMWF
! 06/18/1991
!
! MODIFICATIONS:
!
! 11/14/1995 (Yuhe Liu)
! Changed name, converted to ARPS standard format, and added some
! document
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! rfloat Floating point number to be converted.
!
! kbits Number of bits in computer word.
!
! OUTPUT:
!
! chribm 32 bits IBM floating format in 4-byte character string
!
!-----------------------------------------------------------------------
!
! METHOD:
!
! Floating point number represented as 8 bit signed
! exponent and 24 bit mantissa in integer values.
!
! REFERENCE:
!
! WMO Manual on Codes re GRIB representation.
!
! Common block variables.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL :: rfloat
CHARACTER (LEN=1) :: chribm(4)
INTEGER :: kbits
INTEGER :: iexp
INTEGER :: ISIGN
INTEGER :: kexp
INTEGER :: kmant
REAL :: zeps
REAL :: zref
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF ( kbits == 32 ) THEN
zeps = 1.0E-8
ELSE
zeps = 1.0E-12
END IF
zref = rfloat
!
!-----------------------------------------------------------------------
!
! Sign of value.
!
!-----------------------------------------------------------------------
!
ISIGN = 0
IF (zref < 0.) THEN
ISIGN = 128
zref = - zref
END IF
IF ( zref < 1.0E-20 ) zref = 0.0
!
!-----------------------------------------------------------------------
!
! Convert value of 0.0.
!
!-----------------------------------------------------------------------
!
IF (zref == 0.0) THEN
kexp = 0
kmant = 0
iexp = 0
ELSE
!
!-----------------------------------------------------------------------
!
! Convert other values.
!
! Exponent first
!
!-----------------------------------------------------------------------
!
iexp = INT(ALOG(zref)*(1.0/ALOG(16.0))+64.0+1.0+zeps)
IF ( iexp < 0 ) iexp = 0
IF ( iexp > 127 ) iexp = 127
!
!-----------------------------------------------------------------------
!
! Mantissa.
!
!-----------------------------------------------------------------------
!
kmant = nint (zref/16.0**(iexp-70))
!
!-----------------------------------------------------------------------
!
! Check that mantissa value does not exceed 24 bits.
! 16777215 = 2**24 - 1
!
!-----------------------------------------------------------------------
!
IF (kmant > 16777215) THEN
iexp = iexp + 1
kmant = nint (zref/16.0**(iexp-70))
!
!-----------------------------------------------------------------------
!
! Check for mantissa overflow. If so, set mantissa to zero
!
!-----------------------------------------------------------------------
!
IF (kmant > 16777215) THEN
WRITE (6,'(a,e20.12)') 'Bad mantissa value: overflow', &
rfloat
WRITE (6,'(a,i2,a,i10,a,i10)') &
'isign = ', ISIGN, ' iexp = ', iexp, ' kmant = ', kmant
kmant = 0
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Add sign bit to exponent.
!
!-----------------------------------------------------------------------
!
END IF
kexp = iexp + ISIGN
chribm(1) = CHAR(MOD(kexp,256))
chribm(2) = CHAR(MOD(kmant/65536,256))
chribm(3) = CHAR(MOD(kmant/256, 256))
chribm(4) = CHAR(MOD(kmant, 256))
RETURN
END SUBROUTINE flt2ibm
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE WRTHISHD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE wrthishd(nx,ny,nz,nchanl,grdbas, fmtver, x,y,z, & 1,17
wdlen,nbufsz,mgrib,nbytes)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Write the header of ARPS history dump by using GRIB representation
! to floating point numbers.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 11/27/1995
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! nchanl I/O unit of history dump file
!
! grdbas Flag for grid and base state dump
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
!
! fmtver Format and version identical
!
! wdlen Length in bytes of a machine word
! nbufsz Buffer size of GRIB message
!
! OUTPUT:
!
! mgrib String containing the header of ARPS history dump
! nbytes Length in bytes of the header
!
! WORK ARRAY:
!
!-----------------------------------------------------------------------
!
! The following parameters are passed into this subroutine through
! a common block in globcst.inc, and they determine which
! variables are output.
!
! grdout =0 or 1. If grdout=0, grid variables are not dumped.
! basout =0 or 1. If basout=0, base state variables are not dumped.
! varout =0 or 1. If varout=0, model perturbation variables are not dumped.
! mstout =0 or 1. If mstout=0, water variables are not dumped.
! iceout =0 or 1. If iceout=0, qi, qs and qh are not dumped.
! rainout=0 or 1. If rainout=0, rain variables are not dumped.
! prcout =0 or 1. If prcout=0, precipitation rates are not dumped.
! tkeout =0 or 1. If tkeout=0, tke is not dumped.
! trbout =0 or 1. If trbout=0, turbulence parameter km is not dumped.
! sfcout =0 or 1. If sfcout=0, surface variables are not dumped.
! landout=0 or 1. If landout=0, surface propertty arrays are not dumped.
! radout =0 or 1. If radout=0, radiation arrays are not dumped
! flxout =0 or 1. If flxout=0, surface fluxes are not dumped.
!
! These following parameters are also passed in through common
! blocks in globcst.inc.
!
! runname,curtim,umove,vmove,xgrdorg,ygrdorg
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: nchanl ! I/O unit
INTEGER :: grdbas ! Flag for grid and base state dump
REAL :: x (nx) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL :: y (ny) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL :: z (nz) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
CHARACTER (LEN=40) :: fmtver
INTEGER :: wdlen ! Length in bytes of machine word
INTEGER :: nbufsz
INTEGER :: nbytes ! Length in bytes of the header
CHARACTER (LEN=1) :: mgrib(nbufsz) ! Header of ARPS GRIB file
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
INTEGER :: idummy
INTEGER :: r0exp,r0mant
CHARACTER (LEN=1) :: chrtmp(4)
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
idummy = 0
r0exp = 0 ! for rdummy = 0.0
r0mant = 0 ! for rdummy = 0.0
nbytes = 3
DO i=1,40
mgrib(nbytes+i) = fmtver(i:i)
END DO
nbytes = nbytes + 40
DO i=1,80
mgrib(nbytes+i) = runname(i:i)
END DO
nbytes = nbytes + 80
mgrib(nbytes+1) = CHAR(nocmnt)
nbytes = nbytes + 1
IF( nocmnt > 0 ) THEN
DO j=1,nocmnt
DO i=1,80
mgrib(nbytes+i) = cmnt(j)(i:i)
END DO
nbytes = nbytes + 80
END DO
END IF
CALL flt2ibm(curtim,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
mgrib(nbytes+1) = CHAR(MOD(nx/256,256))
mgrib(nbytes+2) = CHAR(MOD(nx, 256))
nbytes = nbytes + 2
mgrib(nbytes+1) = CHAR(MOD(ny/256,256))
mgrib(nbytes+2) = CHAR(MOD(ny, 256))
nbytes = nbytes + 2
mgrib(nbytes+1) = CHAR(MOD(nz/256,256))
mgrib(nbytes+2) = CHAR(MOD(nz, 256))
nbytes = nbytes + 2
IF( grdbas == 1 ) THEN
mgrib(nbytes+1) = CHAR(1)
mgrib(nbytes+2) = CHAR(1)
mgrib(nbytes+3) = CHAR(0)
mgrib(nbytes+4) = CHAR(mstout)
mgrib(nbytes+5) = CHAR(0)
mgrib(nbytes+6) = CHAR(0)
mgrib(nbytes+7) = CHAR(0)
mgrib(nbytes+8) = CHAR(0)
mgrib(nbytes+9) = CHAR(landout)
mgrib(nbytes+10) = CHAR(totout)
mgrib(nbytes+11) = CHAR(0)
mgrib(nbytes+12) = CHAR(idummy)
mgrib(nbytes+13) = CHAR(idummy)
mgrib(nbytes+14) = CHAR(mapproj)
mgrib(nbytes+15) = CHAR(month)
mgrib(nbytes+16) = CHAR(day)
mgrib(nbytes+17) = CHAR(MOD(year/256,256))
mgrib(nbytes+18) = CHAR(MOD(year, 256))
mgrib(nbytes+19) = CHAR(hour)
mgrib(nbytes+20) = CHAR(minute)
mgrib(nbytes+21) = CHAR(second)
nbytes = nbytes + 21
ELSE
mgrib(nbytes+1) = CHAR(grdout)
mgrib(nbytes+2) = CHAR(basout)
mgrib(nbytes+3) = CHAR(varout)
mgrib(nbytes+4) = CHAR(mstout)
mgrib(nbytes+5) = CHAR(iceout)
mgrib(nbytes+6) = CHAR(trbout)
mgrib(nbytes+7) = CHAR(sfcout)
mgrib(nbytes+8) = CHAR(rainout)
mgrib(nbytes+9) = CHAR(landout)
mgrib(nbytes+10) = CHAR(totout)
mgrib(nbytes+11) = CHAR(tkeout)
mgrib(nbytes+12) = CHAR(idummy)
mgrib(nbytes+13) = CHAR(idummy)
mgrib(nbytes+14) = CHAR(mapproj)
mgrib(nbytes+15) = CHAR(month)
mgrib(nbytes+16) = CHAR(day)
mgrib(nbytes+17) = CHAR(MOD(year/256,256))
mgrib(nbytes+18) = CHAR(MOD(year, 256))
mgrib(nbytes+19) = CHAR(hour)
mgrib(nbytes+20) = CHAR(minute)
mgrib(nbytes+21) = CHAR(second)
nbytes = nbytes + 21
END IF
CALL flt2ibm(umove,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
CALL flt2ibm(vmove,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
CALL flt2ibm(xgrdorg,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
CALL flt2ibm(ygrdorg,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
CALL flt2ibm(trulat1,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
CALL flt2ibm(trulat2,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
CALL flt2ibm(trulon,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
CALL flt2ibm(sclfct,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
mgrib(nbytes+1) = CHAR(r0exp)
mgrib(nbytes+2) = CHAR(r0mant)
mgrib(nbytes+3) = CHAR(r0mant)
mgrib(nbytes+4) = CHAR(r0mant)
nbytes = nbytes + 4
mgrib(nbytes+1) = CHAR(r0exp)
mgrib(nbytes+2) = CHAR(r0mant)
mgrib(nbytes+3) = CHAR(r0mant)
mgrib(nbytes+4) = CHAR(r0mant)
nbytes = nbytes + 4
mgrib(nbytes+1) = CHAR(r0exp)
mgrib(nbytes+2) = CHAR(r0mant)
mgrib(nbytes+3) = CHAR(r0mant)
mgrib(nbytes+4) = CHAR(r0mant)
nbytes = nbytes + 4
mgrib(nbytes+1) = CHAR(r0exp)
mgrib(nbytes+2) = CHAR(r0mant)
mgrib(nbytes+3) = CHAR(r0mant)
mgrib(nbytes+4) = CHAR(r0mant)
nbytes = nbytes + 4
mgrib(nbytes+1) = CHAR(r0exp)
mgrib(nbytes+2) = CHAR(r0mant)
mgrib(nbytes+3) = CHAR(r0mant)
mgrib(nbytes+4) = CHAR(r0mant)
nbytes = nbytes + 4
mgrib(nbytes+1) = CHAR(r0exp)
mgrib(nbytes+2) = CHAR(r0mant)
mgrib(nbytes+3) = CHAR(r0mant)
mgrib(nbytes+4) = CHAR(r0mant)
nbytes = nbytes + 4
mgrib(nbytes+1) = CHAR(r0exp)
mgrib(nbytes+2) = CHAR(r0mant)
mgrib(nbytes+3) = CHAR(r0mant)
mgrib(nbytes+4) = CHAR(r0mant)
nbytes = nbytes + 4
CALL flt2ibm(tstop,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
CALL flt2ibm(thisdmp,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
CALL flt2ibm(latitud,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
CALL flt2ibm(ctrlat,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
CALL flt2ibm(ctrlon,chrtmp,wdlen*8)
mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
IF ( totout /= 0 ) THEN
mgrib(nbytes+1) = CHAR(nstyp)
nbytes = nbytes + 1
!
!-----------------------------------------------------------------------
!
! Reserve 20 integers for future compatibility.
!
!-----------------------------------------------------------------------
!
mgrib(nbytes+1) = CHAR(prcout)
nbytes = nbytes + 1
mgrib(nbytes+1) = CHAR(radout)
nbytes = nbytes + 1
mgrib(nbytes+1) = CHAR(flxout)
nbytes = nbytes + 1
mgrib(nbytes+1) = CHAR(0) ! snowcvr not output anymore
nbytes = nbytes + 1
mgrib(nbytes+1) = CHAR(snowout)
nbytes = nbytes + 1
DO i=6,20
mgrib(nbytes+1) = CHAR(idummy)
nbytes = nbytes + 1
END DO
DO i=0,16,4
mgrib(nbytes+1) = CHAR(r0exp)
mgrib(nbytes+2) = CHAR(r0mant)
mgrib(nbytes+3) = CHAR(r0mant)
mgrib(nbytes+4) = CHAR(r0mant)
nbytes = nbytes + 4
END DO
END IF
!
!-----------------------------------------------------------------------
!
! If grdout=1 or grdbas=1, write out grid variables
!
!-----------------------------------------------------------------------
!
IF(grdout == 1 .OR. grdbas == 1 ) THEN
DO i=1,nx
CALL flt2ibm( mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
END DO
DO j=1,ny
CALL flt2ibm( mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
END DO
DO k=1,nz
CALL flt2ibm( mgrib(nbytes+1) = chrtmp(1)
mgrib(nbytes+2) = chrtmp(2)
mgrib(nbytes+3) = chrtmp(3)
mgrib(nbytes+4) = chrtmp(4)
nbytes = nbytes + 4
END DO
END IF
mgrib(1) = CHAR(MOD(nbytes/65536,256))
mgrib(2) = CHAR(MOD(nbytes/256, 256))
mgrib(3) = CHAR(MOD(nbytes, 256))
WRITE (nchanl) (mgrib(i),i=1,nbytes)
RETURN
END SUBROUTINE wrthishd