c
c	added variance as a type
c
c=================================================================
c
c	change center/model id (in common block)
c
	subroutine setctr(ictr, isctr, imodel)
c
	integer ictr, isctr, imodel
	integer center, subcen, model
	common /gribc/center, subcen, model
	center = ictr
	subcen = isctr
	model = imodel
	return
	end
c=================================================================
c
	block data gribitc
	integer center, subcen, model
	common /gribc/center, subcen, model
c	nmc
	data center/7/
c	reanalysis
	data subcen/1/
c	who knows
	data model/80/
	end
c
c=================================================================
	subroutine ezgbw1(data,nx,ny,iparm,islola,level,
     2	iyear,imonth,iday,ihour,iunit)
c
c
c	ez grib write 1			W. Ebisuzaki
c
c	sets many defaults - so you don't have to 
c		lola - grid,  gaussian (future)
c	writes only 1 field that is for an instantaneous time
c
c
c	real data(nx,ny)	- data to written as grib record
c	integer iparm		- grib parameter number ex. 33 -> U wind
c	logical islola		- .true. -> lola  .false. -> gaussian grid
c	real level		- 0.0 -> 1.0 (sigma), 1.+ ->2000 (prs)
c				  -1 -> surface
c
c
	real data(nx,ny), level
	integer iparm,iyear,imonth,iday,ihour,iunit
	logical islola, ldummy(1)

	call ezgbw(data,nx,ny,ldummy,.false.,iparm,islola,level,
     2	iyear,imonth,iday,ihour,iunit,0,0,'hour','inst')

	return
	end

c=================================================================
	subroutine ezgbw2(data,nx,ny,iparm,islola,level,
     2	iyear,imonth,iday,ihour,iunit,itime1,itime2,timeun,type)
c
c	ez grib write 2			W. Ebisuzaki
c
c	sets many defaults - so you don't have to 
c		lola - grid,  gaussian (future)
c	writes only 1 field that is for a time average/accumulation/inst
c
c
	real data(nx,ny), level
	integer iparm,iyear,imonth,iday,ihour,iunit
	logical islola, ldummy(1)
	integer itime1,itime2
	character*(*) timeun,type

	call ezgbw(data,nx,ny,ldummy,.false., iparm,islola,level,
     2	iyear,imonth,iday,ihour,iunit,itime1,itime2,timeun,type)
	return
	end
c=================================================================
	subroutine ezgbw(data,nx,ny,bit,hasbit,iparm,islola,level,
     2	iyear,imonth,iday,ihour,iunit,itime1,itime2,timeun,type)
c
c
c	ez grib write 				W. Ebisuzaki
c
c	sets many defaults - so you don't have to 
c		lola - grid,  gaussian (future)
c
c	real data(nx,ny)	- data to written as grib record
c	logical bit(*), hasbit  - bitmap of defined values
c	integer iparm		- grib parameter number ex. 33 -> U wind
c	logical islola		- .true. -> lola  .false. -> gaussian grid
c	real level		- 0.0 -> 1.0 (sigma), 1.+ ->2000 (prs)
c				  -1 -> surface
c	integer itime1, itime2	- time average from time1 to time2
c	character*(*) timeun	- 'hour', 'day', 'month'
c	character*(*) type	- 'ave', 'acc', 'inst', 'climo', 'var'
c
c note: type='climo'
c	itime1 = length of mean
c	itime2 = number of years of data
c
	parameter (ilpds=28)
	parameter (mxbit=16)
	parameter (iptv=2)
	parameter (mxnxny=384*190)
c	parameter (ngrib = 100 + ilpds + mxnxny*(mxbit+1)/8+10)
	parameter (ngrib= 1000 + ilpds + mxnxny*(mxbit+1)/8)

	real data(nx,ny), level
	character*1 grib(ngrib)
	integer iyear, imonth, iday, ihour
	character*(*) timeun, type
	logical hasbit, bit(*)

	integer levelt, level1, level2, itime1, itime2
	integer iftu, time1, time2, timerg, nave, nmiss
	integer dscale(255)

	integer gridty, lengrb
	logical init, islola

        integer center, subcen, model
        common /gribc/center, subcen, model

	data init/.false./

	save init, dscale

	if (nx*ny.gt.mxnxny) then
	    write(*,*) '**************'
	    write(*,*) '*** ezgbwr ***'
	    write(*,*) 'array too small: edit, recompile'
	    write(*,*) '*** Fatal Error ***'
	    call exit(99)
	    stop
	endif

	if (.not.init) then
	    call idsdef(iptv,dscale)
	    init = .true.
	endif

	nave = 0
	nmiss = 0

c	gridtyp 0 -> lola   4 -> gaussian grid
	gridty = 0
	if (.not.islola) gridty = 4

c	set colatitude
	if (.not.islola) then
	    if (ny.eq.94) then
c		T62 Gaussian lat
	        colat1 = (90.0-88.542)*3.141592654/180.0
	    elseif (ny.eq.190) then
c		T126 Gaussian lat
	        colat1 = (90.0-89.227)*3.141592654/180.0
	    else
		write(*,*) '*** Missing Gaussian colat in gribit ***'
		write(*,*) '*** please add ***'
		write(*,*) '*** 0 used ***'
	        colat1 = 0.0
	    endif
	endif

c	set level
	level1 = 0
	if (level.ge.2000.0) then
c	    meters above ground - 2000
	    level2 = level - 2000.0
	    levelt = 105
	elseif (level.gt.1.0) then
c	    must be pressure level
	    level2 = level
	    levelt = 100
	elseif (level.ge.0.0 .and. level.le.1.0) then
c	    must be sigma level
	    level2 = 10000.0 * level
	    levelt = 107
	elseif (level.eq.-1.0) then
c	    must be the surface
	    level2 = 0
	    levelt = 1
	elseif (level.gt.-256.0.and.level.lt.-1.0) then
	    rlevel = -level
	    levelt = int(rlevel)

            if (levelt.le.100.or.levelt.eq.102.or.levelt.eq.103 .or.
     1            levelt.eq.105.or.levelt.eq.107.or.levelt.eq.109 .or.
     1            levelt.eq.111.or.levelt.eq.113.or.levelt.eq.125 .or.
     1            levelt.eq.160.or.levelt.eq.200.or.levelt.eq.201) then
                level1 = 0
	        level2 = int(1000*(rlevel - int(rlevel)) + 0.5)
	    else
	        level2 = int(1000*(rlevel - int(rlevel)) + 0.5)
	        level1 = level2 / 256
	        level2 = mod(level2,256)
            endif

	else if (level.eq.-257.0) then
c	    sigma 0..1
	    levelt = 108
	    level2 = 100
	else
c	    ????
	    write(*,*) 'ezgbw2: unknown level?', level
	    call exit(99)
	    stop
	endif

c	for an average/accum.
	if (timeun.eq.'hour'.or.timeun.eq.'HOUR') then
	    iftu = 1
	elseif (timeun.eq.'day'.or.timeun.eq.'DAY') then
	    iftu = 2
	elseif (timeun.eq.'month'.or.timeun.eq.'MONTH') then
	    iftu = 3
	elseif (timeun.eq.'year'.or.timeun.eq.'YEAR') then
	    iftu = 4
	    write(*,*) 'ezwgb: check year parameter!!'
	    stop
	else
	    write(*,*) 'ezwgb timeun:',timeun
	    call exit(99)
	endif

	if (type.eq.'inst'.or.type.eq.'INST'.or.itime1.eq.itime2) then
	    timerg = 10
	    time1 = 0
	    time2 = min0(itime1, itime2)
	else if (type.eq.'acc'.or.type.eq.'ACC') then
	    timerg = 4
	    time1 = min0(itime1, itime2)
	    time2 = max0(itime1, itime2)
	else if (type.eq.'ave'.or.type.eq.'AVE') then
	    timerg = 3
	    time1 = min0(itime1, itime2)
	    time2 = max0(itime1, itime2)
	else if (type.eq.'climo'.or.type.eq.'CLIMO') then
c	    time1 = length of average
c	    time2 = number of years of data
	    timerg = 51
	    time1 = 0
	    time2 = itime1
	    nave = itime2
	else if (type.eq.'var'.or.type.eq.'VAR') then
	    timerg = 118
	    time1 = min0(itime1, itime2)
	    time2 = max0(itime1, itime2)
	else
	     write(*,*) 'ezwgb2 type:',type
	     call exit(99)
	endif

	ids = dscale(iparm)
	if (type.eq.'var'.or.type.eq.'VAR') ids = 2*ids

c	bitmap?
	ibms = 0
	if (hasbit) ibms = 1

	call gribit(data,bit,gridty,nx,ny,mxbit,colat1,ilpds,
     1		iptv,center,subcen,model,ibms,iparm,
     2		levelt, level1, level2,
     3		iyear, imonth, iday, ihour, 
     4		iftu, time1, time2, timerg,
     5		nave, nmiss, ids, grib, lengrb, ierr)

	if (ierr.eq.0) call binwrt(iunit, grib, lengrb)

	return
	end
c
c=====================================================================
c
	subroutine ezgbwk(data,lbm,kpds,kgds,iunit)
c
c	write grib record, kpds, kgds must be defined by caller
c
	parameter (mxbit=16)
	parameter (ilpds = 28)
	parameter (mxnxny=384*190)
c	parameter (ngrib = 100 + ilpds + mxnxny*(mxbit+1)/8+10)
	parameter (ngrib= 1000 + ilpds + mxnxny*(mxbit+1)/8)

	real data(*)
	logical lbm(*)
	integer kpds(*), kgds(*)
	integer center, subcen, model
	character*1 grib(ngrib)

c	data(nx,ny) -- data to write

c	lbm(nx,ny) -- optional bitmap

c	lola/gaussian
	idrt = kgds(1)

c	dim of data(nx,ny)
	nx = kgds(2)
	ny = kgds(3)

c	maximum number of bits

c	colatitude
c	colat1 = (90e3 - (kpds(6)))*acos(-1.0)/(180e3)
	colat1 = (90e3 - (kgds(4)))*acos(-1.0)/(180e3)

c	length of pds (usually 28)
c	ilpds

c	parameter table (usually 1)
	iptv = kpds(19)
c	forecast center
	center = kpds(1)
c	subcenter
	subcen = kpds(23)
c	generating model
	model = kpds(2)
c	has bitmap?
	ibms = mod(kpds(4)/64,2)
	write(*,*) 'ezgbwk1 ibms=',ibms

c	input parameter
	ipu = kpds(5)

c	level indicator
	itl=kpds(6)
c	old code
c        if(itl.LE.100.OR.itl.EQ.102.OR.itl.EQ.103.OR.
c     &    itl.EQ.105.OR.itl.EQ.107.OR.itl.EQ.111.OR.
c     &    itl.EQ.160)  then
c	mod 2/28/95
	if (itl.le.100 .or. itl.eq.102 .or. itl.eq.103 .or.
     1		itl.eq.105 .or. itl.eq.107 .or. itl.eq.109 .or.
     1		itl.eq.111 .or. itl.eq.113 .or. itl.eq.125 .or.
     1		itl.eq.160 .or. itl.eq.200 .or. itl.eq.201) then

            il1=0
            il2=kpds(7)
        else
            il1=mod(kpds(7)/256,256)
            il2=mod(kpds(7),256)
        endif

c	date
c	iyr = kpds(8)
	iyr = kpds(8) + (kpds(21)-1)*100

	imo = kpds(9)
	idy = kpds(10)
	ihr = kpds(11)
c	forecast unit, and times, and units
	iftu = kpds(13)
	ip1 = kpds(14)
	ip2 = kpds(15)
	itr = kpds(16)

c	number in/missing from average
	ina = kpds(17)
	inm = kpds(20)
c	integer decimal scaling
	ids = kpds(22)

c	write(*,*) 'idrt=',idrt,' nx/ny=',nx,ny,' mxbit=',mxbit
c	write(*,*) 'colat1=',colat1,' ilpds=',ilpds,' iptv=',iptv
c	write(*,*) 'icen=',icen,' igen=',igen,' ibms=',ibms,' param=',ipu
c	write(*,*) 'times:',itl,il1,il2
c	write(*,*) 'date:',iyr,imo,idy,ihr
c	write(*,*) 'forcast times:',iftu,ip1,ip2,itr
c	write(*,*) 'no. for ave=',ina,' missing=',inm,' dec scale=',ids
	call gribit(data,lbm,idrt,nx,ny,mxbit,colat1,ilpds,iptv,
     a      center,subcen,model,
     1	    ibms,ipu,itl,il1,il2,iyr,imo,idy,ihr,iftu,ip1,ip2,itr,
     2	    ina,inm,ids,grib,lengrb,ierr)

	if (ierr.eq.0) call binwrt(iunit, grib, lengrb)

	return
	end
c=====================================================================
	subroutine binwrt(iunit,string,n)
c
c	binary write
c
	character*1 string(n)
	write(*,*) '>> binwrt'
	write(iunit) string
	return
	end
c=====================================================================
CFPP$ NOCONCUR R
      SUBROUTINE GRIBIT(F,LBM,IDRT,IM,JM,MXBIT,COLAT1,
     &                  ILPDS,IPTV,ICEN,ISUBCN,IGEN,IBMS,IPU,ITL,
     &                  IL1,IL2,
     &                  IYR,IMO,IDY,IHR,IFTU,IP1,IP2,ITR,INA,INM,IDS,
     &                  GRIB,LGRIB,IERR)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    GRIBIT      CREATE GRIB MESSAGE
C   PRGMMR: IREDELL          ORG: W/NMC23    DATE: 92-10-31
C
C ABSTRACT: CREATE A GRIB MESSAGE FROM A FULL FIELD.
C   AT PRESENT, ONLY GLOBAL LATLON GRIDS AND GAUSSIAN GRIDS ARE ALLOWED.
C
C PROGRAM HISTORY LOG:
C   92-10-31  IREDELL
C
C USAGE:    CALL GRIBIT(F,LBM,IDRT,IM,JM,MXBIT,COLAT1,
C    &                  ILPDS,IPTV,ICEN,IGEN,IBMS,IPU,ITL,IL1,IL2,
C    &                  IYR,IMO,IDY,IHR,IFTU,IP1,IP2,ITR,INA,INM,IDS,
C    &                  GRIB,LGRIB,IERR)
C   INPUT ARGUMENT LIST:
C     F        - REAL (IM*JM) FIELD DATA TO PACK INTO GRIB MESSAGE
C     LBM      - LOGICAL (IM*JM) BITMAP TO USE IF IBMS=1
C     IDRT     - INTEGER DATA REPRESENTATION TYPE
C                (0 FOR LATLON OR 4 FOR GAUSSIAN)
C     IM       - INTEGER LONGITUDINAL DIMENSION
C     JM       - INTEGER LATITUDINAL DIMENSION
C     MXBIT    - INTEGER MAXIMUM NUMBER OF BITS TO USE (0 FOR NO LIMIT)
C     COLAT1   - REAL FIRST COLATITUDE OF GAUSSIAN GRID IF IDRT=4
C     ILPDS    - INTEGER LENGTH OF THE PDS (USUALLY 28)
C     IPTV     - INTEGER PARAMETER TABLE VERSION (USUALLY 1)
C     ICEN     - INTEGER FORECAST CENTER (USUALLY 7)
C     ISUBCN   - INTEGER SUBCENTER
C     IGEN     - INTEGER MODEL GENERATING CODE
C     IBMS     - INTEGER BITMAP FLAG (0 FOR NO BITMAP)
C     IPU      - INTEGER PARAMETER AND UNIT INDICATOR
C     ITL      - INTEGER TYPE OF LEVEL INDICATOR
C     IL1      - INTEGER FIRST LEVEL VALUE (0 FOR SINGLE LEVEL)
C     IL2      - INTEGER SECOND LEVEL VALUE
C    &                  IYR,IMO,IDY,IHR,IFTU,IP1,IP2,ITR,INA,INM,IDS,
C    &                  GRIB,LGRIB,IERR)
C     IYR      - INTEGER YEAR (YYYY)
C     IMO      - INTEGER MONTH
C     IDY      - INTEGER DAY
C     IHR      - INTEGER HOUR
C     IFTU     - INTEGER FORECAST TIME UNIT (1 FOR HOUR)
C     IP1      - INTEGER FIRST TIME PERIOD
C     IP2      - INTEGER SECOND TIME PERIOD (0 FOR SINGLE PERIOD)
C     ITR      - INTEGER TIME RANGE INDICATOR (10 FOR SINGLE PERIOD)
C     INA      - INTEGER NUMBER INCLUDED IN AVERAGE
C     INM      - INTEGER NUMBER MISSING FROM AVERAGE
C     IDS      - INTEGER DECIMAL SCALING
C
C   OUTPUT ARGUMENT LIST:
C     GRIB     - CHARACTER (LGRIB) GRIB MESSAGE
C     LGRIB    - INTEGER LENGTH OF GRIB MESSAGE
C                (NO MORE THAN 100+ILPDS+IM*JM*(MXBIT+1)/8)
C     IERR     - INTEGER ERROR CODE (0 FOR SUCCESS)
C
C SUBPROGRAMS CALLED:
C   GTBITS     - COMPUTE NUMBER OF BITS AND ROUND DATA APPROPRIATELY
C   W3FI72     - ENGRIB DATA INTO A GRIB1 MESSAGE
C
C ATTRIBUTES:
C   LANGUAGE: CRAY FORTRAN
C
C$$$
      REAL F(IM*JM)
      LOGICAL LBM(IM*JM)
      CHARACTER GRIB(*)
      INTEGER IBM(IM*JM*IBMS+1-IBMS),IPDS(25),IGDS(18),IBDS(9)
      REAL FR(IM*JM)
      CHARACTER PDS(ILPDS)
c	write(*,*) '>>gribit'
      NF=IM*JM
      LONI=NINT(360.E3/IM)
      IF(IDRT.EQ.0) THEN
        LAT1=NINT(90.E3)
        LATI=NINT(180.E3/(JM-1))
        IF(IM.EQ.144.AND.JM.EQ.73) THEN
          IGRID=2
        ELSEIF(IM.EQ.360.AND.JM.EQ.181) THEN
          IGRID=3
        ELSE
          IGRID=255
        ENDIF
      ELSEIF(IDRT.EQ.4) THEN
        LAT1=NINT(90.E3-180.E3/ACOS(-1.)*COLAT1)
        LATI=JM/2
        IGRID=255
      ELSE
        IERR=40
        RETURN
      ENDIF
      IPDS(01)=ILPDS    ! LENGTH OF PDS
      IPDS(02)=IPTV     ! PARAMETER TABLE VERSION ID
      IPDS(03)=ICEN     ! CENTER ID
      IPDS(04)=IGEN     ! GENERATING MODEL ID
      IPDS(05)=IGRID    ! GRID ID
      IPDS(06)=1        ! GDS FLAG
      IPDS(07)=IBMS     ! BMS FLAG
      IPDS(08)=IPU      ! PARAMETER UNIT ID
      IPDS(09)=ITL      ! TYPE OF LEVEL ID
      IPDS(10)=IL1      ! LEVEL 1 OR 0
      IPDS(11)=IL2      ! LEVEL 2
c     IPDS(12)=IYR      ! YEAR
      IPDS(12)=mod((IYR-1),100) + 1
      IPDS(13)=IMO      ! MONTH
      IPDS(14)=IDY      ! DAY
      IPDS(15)=IHR      ! HOUR
      IPDS(16)=0        ! MINUTE
      IPDS(17)=IFTU     ! FORECAST TIME UNIT ID
      IPDS(18)=IP1      ! TIME PERIOD 1
      IPDS(19)=IP2      ! TIME PERIOD 2 OR 0
      IPDS(20)=ITR      ! TIME RANGE INDICATOR
      IPDS(21)=INA      ! NUMBER IN AVERAGE
      IPDS(22)=INM      ! NUMBER MISSING
c     IPDS(23)=20       ! CENTURY
      IPDS(23)=((IYR-IPDS(12))/100) + 1
      IPDS(24)=ISUBCN   ! SUBCENTER
      IPDS(25)=IDS      ! DECIMAL SCALING
      IGDS(01)=0        ! NUMBER OF VERTICAL COORDS
      IGDS(02)=255      ! VERTICAL COORD FLAG
      IGDS(03)=IDRT     ! DATA REPRESENTATION TYPE
      IGDS(04)=IM       ! EAST-WEST POINTS
      IGDS(05)=JM       ! NORTH-SOUTH POINTS
      IGDS(06)=LAT1     ! LATITUDE OF ORIGIN
      IGDS(07)=0        ! LONGITUDE OF ORIGIN
      IGDS(08)=128      ! RESOLUTION FLAG
      IGDS(09)=-LAT1    ! LATITUDE OF END
      IGDS(10)=-LONI    ! LONGITUDE OF END
      IGDS(11)=LATI     ! LAT INCREMENT OR GAUSSIAN LATS
      IGDS(12)=LONI     ! LONGITUDE INCREMENT
      IGDS(13)=0        ! SCANNING MODE FLAGS
c     IGDS(14:18)=0     ! NOT USED
	call ifill(igds(14),5,0)
c     IBDS(1:9)=0       ! BDS FLAGS
	call ifill(ibds,9,0)

      NBM=NF
c	write(*,*) 'gribit:ibms=',ibms
      IF(IBMS.NE.0) THEN
        NBM=0
        DO I=1,NF
          IF(LBM(I)) THEN
            IBM(I)=1
            NBM=NBM+1
          ELSE
            IBM(I)=0
          ENDIF
        ENDDO
        IF(NBM.EQ.NF) IPDS(7)=0
      ENDIF
c	write(*,*) 'gribit2 nbm=',nbm
      IF(NBM.EQ.0) THEN
        DO I=1,NF
          FR(I)=0.
        ENDDO
        NBIT=0
      ELSE
c	write(*,*)'gribit3 ipds(7)=',ipds(7),' ids=', ids
c	write(*,*) 'f(1), fr(1)=',f(1),fr(1),' nbit=',nbit
c	write(*,*) 'nf=#=',nf
	fmin = f(1)
	fmax = fmin
	do ii = 1, nf
	    fmin = amin1(fmin,f(ii))
	    fmax = amax1(fmax,f(ii))
	enddo
c	write(*,*) 'fmin/fmax=',fmin,fmax

        CALL GTBITS(IPDS(7),IDS,NF,IBM,F,FR,FMIN,FMAX,NBIT)
c	write(*,*) 'fmin/max=',fmin,fmax,' nbit=',nbit
        IF(MXBIT.GT.0) NBIT=MIN(NBIT,MXBIT)
      ENDIF
c	write(*,*)'gribit4'
      CALL W3FI72(0,FR,0,NBIT,0,IPDS,PDS,
     &            1,255,IGDS,0,0,IBM,NF,IBDS,
     &            NFO,GRIB,LGRIB,IERR)
c	write(*,*)'gribit5'
      RETURN
      END

CFPP$ NOCONCUR R
      SUBROUTINE GTBITS(IBM,IDS,LEN,MG,G,GROUND,GMIN,GMAX,NBIT)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    GTBITS      COMPUTE NUMBER OF BITS AND ROUND FIELD.
C   PRGMMR: IREDELL          ORG: W/NMC23    DATE: 92-10-31
C
C ABSTRACT: THE NUMBER OF BITS REQUIRED TO PACK A GIVEN FIELD
C   AT A PARTICULAR DECIMAL SCALING IS COMPUTED USING THE FIELD RANGE.
C   THE FIELD IS ROUNDED OFF TO THE DECIMAL SCALING FOR PACKING.
C   THE MINIMUM AND MAXIMUM ROUNDED FIELD VALUES ARE ALSO RETURNED.
C   GRIB BITMAP MASKING FOR VALID DATA IS OPTIONALLY USED.
C
C PROGRAM HISTORY LOG:
C   92-10-31  IREDELL
C
C USAGE:    CALL GTBITS(IBM,IDS,LEN,MG,G,GMIN,GMAX,NBIT)
C   INPUT ARGUMENT LIST:
C     IBM      - INTEGER BITMAP FLAG (=0 FOR NO BITMAP)
C     IDS      - INTEGER DECIMAL SCALING
C                (E.G. IDS=3 TO ROUND FIELD TO NEAREST MILLI-VALUE)
C     LEN      - INTEGER LENGTH OF THE FIELD AND BITMAP
C     MG       - INTEGER (LEN) BITMAP IF IBM=1 (0 TO SKIP, 1 TO KEEP)
C     G        - REAL (LEN) FIELD
C
C   OUTPUT ARGUMENT LIST:
C     GROUND   - REAL (LEN) FIELD ROUNDED TO DECIMAL SCALING
C     GMIN     - REAL MINIMUM VALID ROUNDED FIELD VALUE
C     GMAX     - REAL MAXIMUM VALID ROUNDED FIELD VALUE
C     NBIT     - INTEGER NUMBER OF BITS TO PACK
C
C SUBPROGRAMS CALLED:
C   ISRCHNE  - FIND FIRST VALUE IN AN ARRAY NOT EQUAL TO TARGET VALUE
C
C ATTRIBUTES:
C   LANGUAGE: CRAY FORTRAN
C
C$$$
      DIMENSION MG(LEN),G(LEN),GROUND(LEN)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  ROUND FIELD AND DETERMINE EXTREMES WHERE BITMAP IS ON
      DS=10.**IDS
      IF(IBM.EQ.0) THEN
        GROUND(1)=NINT(G(1)*DS)/DS
        GMAX=GROUND(1)
        GMIN=GROUND(1)
        DO I=2,LEN
          GROUND(I)=NINT(G(I)*DS)/DS
          GMAX=MAX(GMAX,GROUND(I))
          GMIN=MIN(GMIN,GROUND(I))
        ENDDO
      ELSE
        I1=ISRCHNE(LEN,MG,1,0)
        IF(I1.GT.0.AND.I1.LE.LEN) THEN
          GROUND(I1)=NINT(G(I1)*DS)/DS
          GMAX=GROUND(I1)
          GMIN=GROUND(I1)
          DO I=I1+1,LEN
            IF(MG(I).NE.0) THEN
              GROUND(I)=NINT(G(I)*DS)/DS
              GMAX=MAX(GMAX,GROUND(I))
              GMIN=MIN(GMIN,GROUND(I))
            ENDIF
          ENDDO
        ELSE
          GMAX=0.
          GMIN=0.
        ENDIF
      ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  COMPUTE NUMBER OF BITS
      NBIT=LOG((GMAX-GMIN)*DS+0.9)/LOG(2.)+1.
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      END
