	program gribslp
c
c	a program to find sea-level pressure
c
c	uses same method as GDAS
c
c	this program is built along the lines of
c	grib_covar.f
c
c	file 50: seq grib file-x -- Z1000
c	file 51: index grib file-x
c	file 52: seq grib file- Z500
c	file 53: index grib file-y
c	file 54: list of times [yyyymmddhh]
c
c	file 1: output file seq grib - slp
c
	parameter (maxv=100000)

	logical lbmsx(maxv), lbmsy(maxv)
	real z1000(maxv), z500(maxv), slp(maxv)
	integer icount

	integer kpdsx(100),kgdsx(100),reclx,offsetx
	integer kpdsy(100),kgdsy(100),recly,offsety
	integer kpdsx6, kpdsx7
	integer kpdsy6, kpdsy7
	integer kpdsz6, kpdsz7

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

	write(*,*) 'grib_slp'

	call baopen(50,'fort.50',iret)
	if (iret.ne.0) then
	    write(*,*) 'error: could not open unit 50'
	    call exit(8)
	endif

	call baopen(51,'fort.51',iret)
	if (iret.ne.0) then
	    write(*,*) 'error: could not open unit 51'
	    call exit(8)
	endif

	call baopen(52,'fort.52',iret)
	if (iret.ne.0) then
	    write(*,*) 'error: could not open unit 52'
	    call exit(8)
	endif

	call baopen(53,'fort.53',iret)
	if (iret.ne.0) then
	    write(*,*) 'error: could not open unit 52'
	    call exit(8)
	endif

c	find kpds(5..7) parameters

	call atok5('HGT',idx)
	call atok5('HGT',idy)
	call atok5('PRES',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,1000.0)
	kpdsx6=kpdsx(6)
	kpdsx7=kpdsx(7)
	call setlev(kpdsy,kgdsy,500.0)
	kpdsy6=kpdsy(6)
	kpdsy7=kpdsy(7)
c	mean sea level
	kpdsz6 = 102
	kpdsz7 = 0
	levelz = -102.0
	icount = 0

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

250	read(54,*,end=500) time
	write(*,*) 'time=',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) = kpdsy6
            kpdsy(7) = kpdsy7
            kpdsy(8) = mod(year-1,100)+1
            kpdsy(21) = (year-kpdsy(8))/100 + 1
            kpdsy(9) = mon
            kpdsy(10) = day
            kpdsy(11) = hour

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

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

	write(*,*) '>>rdgb1'
            call rdgb(50,reclx,offsetx,kpdsx,kgdsx,n,lbmsx,z1000)
	write(*,*) '>>rdgb2'
            call rdgb(52,recly,offsety,kpdsy,kgdsy,n,lbmsy,z500)

	    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
c		*** gdas - in mb ***
c		slp(i)=1000.0*exp(z1000(i)/
c     2			(1.5422885*(z500(i)-z1000(i))))
c		*** grib - in Pascals ***
		slp(i)=100000.0*exp(z1000(i)/
     2			(1.5422885*(z500(i)-z1000(i))))
450	    continue

c	    set center, subcenter, model
	    call setctr(kpdsx(1), kpdsx(23), kpdsx(2))

	write(*,*) '>>ezgbw1'
            call ezgbw1(slp,nx,ny,newid,.true.,levelz,
     1              year,mon,day,hour,1)


	    icount = icount + 1
	goto 250

500	continue

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