	program gribrh
c
c	a program to go from rel humdity/temp to specific hum
c
c	this program is built along the lines of
c	grib_covar.f
c
c	file 50: seq grib file-x -- temp
c	file 51: index grib file-temp
c	file 52: seq grib file- rh
c	file 53: index grib file-rh
c	file 54: list of times [yyyymmddhh]
c	file 55: level
c
c	file 1: output file seq grib - specific humidity
c
c	v1.1 - dynamic sizing of precision
c
	parameter (maxv=100000)

	logical lbmsx(maxv), lbmsy(maxv)
	real datax(maxv), datay(maxv)
	integer icount

	integer kpdsx(100),kgdsx(100),reclx,offsetx
	integer kpdsy(100),kgdsy(100),recly,offsety
	integer kpdsx6, kpdsx7

	integer idx, idy, newid, time, time2, year, mon, hour, day
	real level
	logical hasbtx, hasbty

	call baopen(50,'fort.50',iret)
	if (iret.ne.0) then
	    write(*,*) 'error opening fort.50'
	    stop 7
	endif
	call baopen(51,'fort.51',iret)
	if (iret.ne.0) then
	    write(*,*) 'error opening fort.51'
	    stop 7
	endif
	call baopen(52,'fort.52',iret)
	if (iret.ne.0) then
	    write(*,*) 'error opening fort.52'
	    stop 7
	endif
	call baopen(53,'fort.53',iret)
	if (iret.ne.0) then
	    write(*,*) 'error opening fort.53'
	    stop 7
	endif

	read(55,*) level

c	find kpds(5..7) parameters

	call atok5('TMP',idx)
	call atok5('RH',idy)
	call atok5('SPFH',newid)

	if (idx.eq.-1.or.idy.eq.-1.or.newid.eq.-1) then
	    write(*,*) 'illegal variable name'
	    call exit(8)
	endif

	call setlev(kpdsx,kgdsx,level)
	kpdsx6=kpdsx(6)
	kpdsx7=kpdsx(7)

c	*** read index file ***
	inumx = 0
	inumy = 0
	icount = 0

250	read(54,*,end=500) time
	    if (time.le.0) goto 500

	    time2 = time
	    hour = mod(time,100)
	    time = time/100
	    day = mod(time,100)
	    time = time/100
	    mon = mod(time,100)
	    time = time/100
	    year = mod(time,10000)

            do 260 i = 1, 100
                kpdsx(i) = -1
                kpdsy(i) = -1
260         continue
            kpdsx(5) = idx
            kpdsx(6) = kpdsx6
            kpdsx(7) = kpdsx7
            kpdsx(8) = mod(year-1,100)+1
            kpdsx(21) = (year-kpdsx(8))/100 + 1
            kpdsx(9) = mon
            kpdsx(10) = day
            kpdsx(11) = hour
 
            kpdsy(5) = idy
            kpdsy(6) = kpdsx6
            kpdsy(7) = kpdsx7
            kpdsy(8) = mod(year-1,100)+1
            kpdsy(21) = (year-kpdsy(8))/100 + 1
            kpdsy(9) = mon
            kpdsy(10) = day
            kpdsy(11) = hour

            call pdscan(51,kpdsx,inumx,offsetx,reclx)
            if (reclx.eq.0) then
                write(*,*) 'missing x, time ', time2
		call exit(8)
                goto 250
            endif

            call pdscab(53,kpdsy,inumy,offsety,recly)
            if (recly.eq.0) then
                write(*,*) 'missing y, time ', time2
		call exit(8)
                goto 250
            endif

            call rdgb(50,reclx,offsetx,kpdsx,kgdsx,ndatax,lbmsx,datax)
            call rdgb(52,recly,offsety,kpdsy,kgdsy,ndatay,lbmsy,datay)

	    call getbit(kpdsx, kgdsx, hasbtx)
	    call getbit(kpdsy, kgdsy, hasbty)

	    nx = kgdsx(2)
	    ny = kgdsx(3)
	    nxny = nx * ny
	    if (nxny.ne.kgdsy(2)*kgdsy(3)) then
		write(*,*) 'unequal dimensions!'
		call exit(8)
	    endif

	    if (hasbtx.or.hasbty) then
		write(*,*) 'expecting no bitmap - but has bitmap'
		write(*,*) ' modify code'
		call exit(8)
	    endif

	    do 450 i = 1, nxny
		call trh2sh(datax(i),datay(i),level,shum)
	        datax(i) = shum
450	    continue

c	    old way -- fixed scaling
c           call ezgbw1(datax,nx,ny,newid,.true.,level,
c    1              year,mon,day,hour,1)

c	    new way -- find best decimal scaling factor

	    shmin = datax(1)
	    shmax = datax(1)
	    do i = 2, nxny
		shmin = amin1(shmin,datax(i))
		shmax = amax1(shmax,datax(i))
	    enddo
	    shdiff = shmax - shmin
	    iscale = 0
	    if (shdiff .eq. 0.0) then
		write(*,*) 'constant spec. humidity field'
	    else
460		if (shdiff.lt.999.9) then
		    iscale = iscale + 1
		    shdiff = shdiff * 10.0
		    goto 460
		endif

470		if (shdiff.gt.9999.9) then
		    iscale = iscale - 1
		    shdiff = shdiff * 0.1
		    goto 470
		endif
	    endif

		
	    kpdsx(5) = newid
	    kpdsx(22) = iscale
            call ezgbwk(datax,lbmsx,kpdsx,kgdsx,1)

	    icount = icount + 1
	goto 250

500	continue

	write(*,*) 'no. spec hum fields calculated = ',icount
	if (icount.eq.0) call exit(8)
	call exit(0)
	stop
	end
