!
!   Author: Suranjana Saha
c--------------------------------------------------------------------------
c    compile using : make NAME=bufrslupa
c--------------------------------------------------------------------------
c     the mnemonic tables are in /nwprod/prepobs/ucl/prep.bufrtable
c---------------------------------------------------------------------------
       CHARACTER*1   kdbug,kdbuga,kall,kprs,kmas,kwnd
       CHARACTER*2   ktol
       CHARACTER*6   catname
       CHARACTER*255 infile,outfile
       CHARACTER*5   kst,klev
       CHARACTER*8   klonw,klone,klatn,klats,kdhr
       CHARACTER*2   kchr
       real(8)       bmiss

       dimension plev(28)

       bmiss=10e10; call setbmiss(bmiss) ! this sets bufrlib missing value to 10e10

       call getenv("idbug",kdbug)
       read(kdbug,'(i1)') idbug
       write(*,*) "idbug= ",idbug
c
       call getenv("idbuga",kdbuga)
       read(kdbuga,'(i1)') idbuga
       write(*,*) "idbuga= ",idbuga
c
       call getenv("iall",kall)
       read(kall,'(i1)') iall
       write(*,*) "iall= ",iall
C
       call getenv("iwnd",kwnd)
       read(kwnd,'(i1)') iwnd
       write(*,*) "iwnd= ",iwnd
C
       call getenv("imas",kmas)
       read(kmas,'(i1)') imas
       write(*,*) "imas= ",imas
C
       call getenv("fchr",kchr)
       read(kchr,'(i2)') ichr
       fchr=float(ichr)
       write(*,*) "fchr= ",fchr
C
       call getenv("fdhr",kdhr)
       read(kdhr,'(f8.1)') fdhr
       write(*,*) "fdhr= ",fdhr
c
       call getenv("nst",kst)
       read(kst,'(i5)') nst
       write(*,*) "nst= ",nst
c
       call getenv("infile",infile)
       write(*,*) "infile= ",infile
c
       call getenv("outfile",outfile)
       write(*,*) "outfile= ",outfile
c
       call getenv("catname",catname)
       write(*,*) "catname= ",catname
c
       call getenv("flonw",klonw)
       read(klonw,'(f8.1)') flonw
       write(*,*) "flonw= ",flonw
c
       call getenv("flone",klone)
       read(klone,'(f8.1)') flone
       write(*,*) "flone= ",flone
c
       call getenv("flatn",klatn)
       read(klatn,'(f8.1)') flatn
       write(*,*) "flatn= ",flatn
c
       call getenv("flats",klats)
       read(klats,'(f8.1)') flats
       write(*,*) "flats= ",flats
c
       call getenv("ilev",klev)
       read(klev,'(i4)') ilev
       write(*,*) "ilev= ",ilev
c
       call getenv("iprs",kprs)
       read(kprs,'(i1)') iprs
       write(*,*) "iprs= ",iprs
c
       call getenv("itol",ktol)
       read(ktol,'(i2)') itol
       write(*,*) "itol= ",itol
c
c...  define mandatory levels....
       if(iprs.eq.1) then
         nlevs=21
         plev(1)= 1000
         plev(2)= 925
         plev(3)= 850
         plev(4)= 700
         plev(5)= 500
         plev(6)= 400
         plev(7)= 300
         plev(8)= 250
         plev(9)=200
         plev(10)=150
         plev(11)=100
         plev(12)=70
         plev(13)=50
         plev(14)=30
         plev(15)=20
         plev(16)=10
         plev(17)=7
         plev(18)=5
         plev(19)=3
         plev(20)=2
         plev(21)=1
       endif
c
c...  define sigma levels....
       if(iprs.eq.0) then
         nlevs=28
         plev(1)= 995
         plev(2)= 982
         plev(3)= 964
         plev(4)= 943
         plev(5)= 916
         plev(6)= 884
         plev(7)= 846
         plev(8)= 801
         plev(9)= 751
         plev(10)=694
         plev(11)=633
         plev(12)=568
         plev(13)=502
         plev(14)=436
         plev(15)=372
         plev(16)=312
         plev(17)=258
         plev(18)=210
         plev(19)=168
         plev(20)=133
         plev(21)=103
         plev(22)= 78
         plev(23)= 58
         plev(24)= 42
         plev(25)= 29
         plev(26)= 18
         plev(27)= 10
         plev(28)=  3
       endif
c
         if(iprs.eq.2) nlevs=ilev
c
       if(flonw.lt.0.) flonw=360.+flonw
       if(flone.lt.0.) flone=360.+flone
       fhrs=fchr-fdhr
       fhre=fchr+fdhr
       print *,fchr,fdhr,fhrs,fhre
c
       call stats(infile,outfile,flonw,flone,flatn,flats,
     * iprs,nlevs,plev,idbug,idbuga,itol,iall,nst,fhrs,fhre,
     * catname,imas,iwnd)
c
       stop
       end
       subroutine stats(infile,outfile,flonw,flone,flatn,flats,
     * iprs,nlevs,plev,idbug,idbuga,itol,iall,nst,fhrs,fhre,
     * catname,imas,iwnd)
c
       CHARACTER*80 HSTR,PSTR,ZSTR,TSTR,QSTR,USTR,VSTR,PSTR1
       CHARACTER*8  SUBSET
C
       character*8   tstn(nst),stnidc,id
       CHARACTER*255 infile,outfile
c
       CHARACTER*10 indate
       CHARACTER*6  catname
       CHARACTER*8  STNID,idj,idi
       CHARACTER*1  labn
C
       real(8)    HDR(14)
       EQUIVALENCE  (HDR(1),STNID)
       real(8)    POB(4,255),ZOB(4,255),TOB(4,255),QOB(4,255)
       real(8)    UOB(4,255),VOB(4,255),PSOB(4)
       DIMENSION    PPR(4,nlevs),ZPR(4,nlevs),TPR(4,nlevs),QPR(4,nlevs)
       DIMENSION    UPR(4,nlevs),VPR(4,nlevs),PSPR(4)
       DIMENSION    tlon(nst),tlat(nst),thr(nst),typ(nst),elvt(nst)
       DIMENSION    typint(nst)
C
       DIMENSION PO(nlevs,nst)
       DIMENSION POM(nlevs,nst),POW(nlevs,nst)
       DIMENSION PSO(nst),PSA(nst),PSF(nst),PSQ(nst)
       DIMENSION ZO(nlevs,nst),ZA(nlevs,nst),ZF(nlevs,nst),ZQ(nlevs,nst)
       DIMENSION TO(nlevs,nst),TA(nlevs,nst),TF(nlevs,nst),TQ(nlevs,nst)
       DIMENSION QO(nlevs,nst),QA(nlevs,nst),QF(nlevs,nst),QQ(nlevs,nst)
       DIMENSION UO(nlevs,nst),UA(nlevs,nst),UF(nlevs,nst),WQ(nlevs,nst)
       DIMENSION VO(nlevs,nst),VA(nlevs,nst),VF(nlevs,nst)
c
       DIMENSION POG(nlevs)
       DIMENSION ZOG(nlevs),ZAG(nlevs),ZFG(nlevs),ZQG(nlevs)
       DIMENSION TOG(nlevs),TAG(nlevs),TFG(nlevs),TQG(nlevs)
       DIMENSION QOG(nlevs),QAG(nlevs),QFG(nlevs),QQG(nlevs)
       DIMENSION UOG(nlevs),UAG(nlevs),UFG(nlevs),WQG(nlevs)
       DIMENSION VOG(nlevs),VAG(nlevs),VFG(nlevs)
c
       REAL    plev(nlevs)
       integer kfillm(nlevs,nst),kfillw(nlevs,nst)
       integer mmfill(nst),mwfill(nst)
       integer kmmax(nst),kwmax(nst)
c
       DATA HSTR
     ./'SID XOB YOB DHR ELV TYP T29 ITP SQN RQM DUP PRG SRC RUD'/
C
       DATA STNID/'        '/
c
c...   get the ob, guess, analysis and quality mark ....
c
       DATA PSTR1 /'POB PAN PFC PQM CAT=0'/
       DATA PSTR /'POB PAN PFC PQM'/
       DATA ZSTR /'ZOB ZAN ZFC ZQM'/
       DATA TSTR /'TOB TAN TFC TQM'/
       DATA QSTR /'QOB QAN QFC QQM'/
c      DATA USTR /'DDO UAN UFC WQM'/
c      DATA VSTR /'FFO VAN VFC WQM'/
       DATA USTR /'UOB UAN UFC WQM'/
       DATA VSTR /'VOB VAN VFC WQM'/
C
       DATA BMISS /10E10/
       DATA rmiss /99999./
       DATA LUNIN /11/
C----------------------------------------------------
       indate='          '
c
       print *,' infile ',infile
       print *,' outfile ',outfile
       print *,' catname ',catname
       print *,' flonw ',flonw
       print *,' flone ',flone
       print *,' flatn ',flatn
       print *,' flats ',flats
       iodd=0
       if(flonw.gt.flone) then
       iodd=1
       flone=flone+360.
       print *,' iodd ',iodd
       print *,' flone ',flone
       endif
c
       open(11,file=infile,form='unformatted')
c
       open(51,file=outfile,form='unformatted')
c
C  INITIALIZE
C  ----------
       DO IST=1,nst
       PSO(IST)   = bmiss
       PSA(IST)   = bmiss
       PSF(IST)   = bmiss
       PSQ(IST)   = bmiss
       ENDDO
C
       DO IST=1,nst
       DO L=1,nlevs
       PO(L,IST)   = bmiss
       POM(L,IST)   = bmiss
       POW(L,IST)   = bmiss
       ZO(L,IST)   = bmiss
       TO(L,IST)   = bmiss
       QO(L,IST)   = bmiss
       UO(L,IST)   = bmiss
       VO(L,IST)   = bmiss
       ZA(L,IST)   = bmiss
       TA(L,IST)   = bmiss
       QA(L,IST)   = bmiss
       UA(L,IST)   = bmiss
       VA(L,IST)   = bmiss
       ZF(L,IST)   = bmiss
       TF(L,IST)   = bmiss
       QF(L,IST)   = bmiss
       UF(L,IST)   = bmiss
       VF(L,IST)   = bmiss
       ZQ(L,IST)   = bmiss
       TQ(L,IST)   = bmiss
       QQ(L,IST)   = bmiss
       WQ(L,IST)   = bmiss
       KFILLM(L,IST)   = 0
       KFILLW(L,IST)   = 0
       ENDDO
       MMFILL(IST)   = 0
       MWFILL(IST)   = 0
       ENDDO
C
       nvar=0
       if((imas.eq.1).and.(iwnd.eq.1)) nvar=1
       print *,'nvar ',nvar
c
       CALL DATELEN(10)
       CALL OPENBF(LUNIN,'IN ',LUNIN)
c
       IST=0
!      DO WHILE(IREADFT(LUNIN,SUBSET,IDATE).EQ.0)
       DO WHILE(IREADMG(LUNIN,SUBSET,IDATE).EQ.0)
       write(indate,'(i10)') idate
c
c... look for only upper-air data... if you want anything else, the
c    following headers can be used :
c
c    ADPSFC SFCSHP  AIRCFT AIRCAR SATWND UPABOG 
c    SFCBOG SATEMP SATBOG SYNDAT PROFLR SPSSMI
c
       IF(SUBSET.EQ.catname) THEN
C
         DO WHILE(IREADSB(LUNIN).EQ.0)
c
           CALL UFBINT(LUNIN,HDR,14,  1,NLEV,HSTR)
c
           hr=hdr(4)
           flon=hdr(2)
           flat=hdr(3)
           typx=hdr(6)
           typintx=hdr(8)
           elvx=hdr(5)
           kkx = idint(hdr(6)*.01)
c
       if((nvar.eq.0).and.(imas.eq.1).and.(kkx.ne.1)) go to 443
       if((nvar.eq.0).and.(iwnd.eq.1).and.(kkx.ne.2)) go to 443
c
c... check for window hour...
       hrx=hr*10.
       fhrsx=fhrs*10.
       fhrex=fhre*10.
           if((hrx.ge.fhrsx).and.(hrx.le.fhrex)) then
c
c... look for all stations inside the lat/lon box specified...
c
       flonx=flon
       if((iodd.eq.1).and.(flonx.lt.flonw)) flonx=flonx+360.
c
       if((flonx.ge.flonw).and.(flonx.le.flone).and.
     *    (flat.le.flatn).and.(flat.ge.flats)) then
c
c...  first time....
       if(ist.eq.0) then
       IST=IST+1
       tstn(IST)='        '
       tstn(IST)=stnid
       tlon(IST)=flon
       tlat(IST)=flat
       thr(IST)=hr
       elvt(IST)=elvx
       typint(IST)=typintx
       if(kkx.eq.1) then
       typ(IST)=typx-100
       else
       typ(IST)=typx-200
       endif
       go to 10
       endif
c
c..   check for unique station....
       if(nvar.eq.1) then
       do n=1,ist
       dlon=abs(flon-tlon(n))
       dlat=abs(flat-tlat(n))
       dthr=abs(hr-thr(n))
       if((dlat.le.0.01).and.(dlon.le.0.01).and.(dthr.le.0.001)) then
       mmf=mmfill(n)
       mwf=mwfill(n)
c
c...  check if both wind and mass report have been read..
       if((mmf.eq.1).and.(mwf.eq.1)) go to 443
       if((mmf.eq.0).and.(mwf.eq.0)) go to 12
       if((mmf.eq.1).or.(mwf.eq.1))  go to 10
c
       endif
       enddo
       endif
c
c...  if not, then add a new station....
 12    continue
       IST=IST+1
       if(ist.gt.nst) then
       print *,'number of stations exceeded'
       print *,'increase number to more than ',nst
       go to 1555
       endif
       tstn(IST)='        '
       tstn(IST)=stnid
       tlon(IST)=flon
       tlat(IST)=flat
       thr(IST)=hr
       elvt(IST)=elvx
       typint(IST)=typintx
       if(kkx.eq.1) then
       typ(IST)=typx-100
       else
       typ(IST)=typx-200
       endif
c
 10    continue
c
       do nn=1,4
       psob(nn)=bmiss
       enddo
c
       if((kkx.eq.1).and.(imas.eq.1)) then
       do ll=1,255
       do nn=1,4
       pob(nn,ll)=bmiss
       zob(nn,ll)=bmiss
       tob(nn,ll)=bmiss
       qob(nn,ll)=bmiss
       enddo
       enddo
       do ll=1,nlevs
       do nn=1,4
       ppr(nn,ll)=bmiss
       zpr(nn,ll)=bmiss
       tpr(nn,ll)=bmiss
       qpr(nn,ll)=bmiss
       enddo
       enddo
       endif
c
       if((kkx.eq.2).and.(iwnd.eq.1)) then
       do ll=1,255
       do nn=1,4
       uob(nn,ll)=bmiss
       vob(nn,ll)=bmiss
       enddo
       enddo
       do ll=1,nlevs
       do nn=1,4
       ppr(nn,ll)=bmiss
       zpr(nn,ll)=bmiss
       tpr(nn,ll)=bmiss
       qpr(nn,ll)=bmiss
       upr(nn,ll)=bmiss
       vpr(nn,ll)=bmiss
       enddo
       enddo
       endif
c
c... get surface pressure ob
       if(kkx.eq.1) then
       CALL UFBINT(LUNIN,PSOB, 4,1,NLEV,PSTR1)
       if(nlev.eq.1) then
       do i=1,4
       if(psob(i).ne.bmiss) pspr(i)=psob(i)
       if(pspr(i).gt.rmiss) pspr(i)=bmiss
       enddo
       pso(IST) =  pspr(1)
       psa(IST) =  pspr(2)
       psf(IST) =  pspr(3)
       psq(IST) =  pspr(4)
       endif
       write(6,610) ist,indate,stnid,flon,flat,typx,typintx,elvx,hr,
     * pso(ist),psa(ist),psf(ist),psq(ist)
       endif
c
c... get pressure at all levels..
       CALL UFBINT(LUNIN,POB, 4,255,NLEV,PSTR)
C..
       lx=0
       DO l=1,nlev
c
       pdata=POB(1,l)
c
       if(iprs.le.1) then
       if(kkx.eq.1) k=GETLEV(pdata,plev,itol,nlevs,kfillm,ist,nst)
       if(kkx.eq.2) k=GETLEV(pdata,plev,itol,nlevs,kfillw,ist,nst)
       else
       k=0
       lx=lx+1
       if(lx.gt.nlevs) then
       print *,'number of levels more than ilev ..cap at level ',pdata
       go to 443
       endif
       if((kkx.eq.1).and.(kfillm(lx,ist).eq.0)) k=lx
       if((kkx.eq.2).and.(kfillw(lx,ist).eq.0)) k=lx
       endif
c
c..  check proper level...
c
       if(k.gt.0) then
c
       if(iprs.le.1) then
       PO(K,IST) =  plev(k)
       else
       if(nvar.eq.0) then
       po(k,ist)=pob(1,l)
       else
       if(kkx.eq.1) POM(K,IST) =  pob(1,l)
       if(kkx.eq.2) POW(K,IST) =  pob(1,l)
       endif
       endif
c
c...  check for mass report....
       if((kkx.eq.1).and.(imas.eq.1)) then
       kmmax(ist)=k
       kfillm(k,ist)=1
       mmfill(ist)=1
c
       CALL UFBINT(LUNIN,ZOB, 4,255,NLEV,ZSTR)
       CALL UFBINT(LUNIN,TOB, 4,255,NLEV,TSTR)
       CALL UFBINT(LUNIN,QOB, 4,255,NLEV,QSTR)
c
       do i=1,4
       if(zob(i,l).ne.bmiss) zpr(i,k)=zob(i,l)
       if(zpr(i,k).gt.rmiss) zpr(i,k)=bmiss
       if(tob(i,l).ne.bmiss) tpr(i,k)=tob(i,l)
       if(tpr(i,k).gt.rmiss) tpr(i,k)=bmiss
       if(qob(i,l).ne.bmiss) qpr(i,k)=qob(i,l)
       if(qpr(i,k).gt.rmiss) qpr(i,k)=bmiss
       enddo
c
       ZO(K,IST) =  ZPR(1,K)
       ZQ(K,IST) =  ZPR(4,K)
       TQ(K,IST) =  TPR(4,K)
       QQ(K,IST) =  QPR(4,K)
c
       TO(K,IST) =  TPR(1,K)
       QO(K,IST) =  QPR(1,K)
       if(TO(K,IST).ne.bmiss) TO(K,IST) =  TO(K,IST)+273.159
       IF(QO(K,IST).ne.bmiss) QO(K,IST) =  QO(K,IST)*1.E-3
       TV =  TO(K,IST)
       Q  =  QO(K,IST)
       if(Q.ne.bmiss) Q=Q*1.E-3
       if((Q.ne.bmiss).and.(TV.ne.bmiss)) TO(K,IST)=TV/(1.+.608*Q)
c
       if(iall.eq.1) then
       ZA(K,IST) =  ZPR(2,K)
       ZF(K,IST) =  ZPR(3,K)
c
       TA(K,IST) =  TPR(2,K)
       QA(K,IST) =  QPR(2,K)
       if(TA(K,IST).ne.bmiss) TA(K,IST) =  TA(K,IST)+273.159
       IF(QA(K,IST).ne.bmiss) QA(K,IST) =  QA(K,IST)*1.E-3
       TV = TA(K,IST)
       Q = QA(K,IST)
       if(Q.ne.bmiss) Q=Q*1.E-3
       if((Q.ne.bmiss).and.(TV.ne.bmiss)) TA(K,IST)=TV/(1.+.608*Q)
c
       TF(K,IST) =  TPR(3,K)
       QF(K,IST) =  QPR(3,K)
       if(TF(K,IST).ne.bmiss) TF(K,IST) =  TF(K,IST)+273.159
       IF(QF(K,IST).ne.bmiss) QF(K,IST) =  QF(K,IST)*1.E-3
       TV =  TF(K,IST)
       Q  =  QF(K,IST)
       if(Q.ne.bmiss) Q=Q*1.E-3
       if((Q.ne.bmiss).and.(TV.ne.bmiss)) TF(K,IST)=TV/(1.+.608*Q)
       endif
c
       if(idbug.eq.1) then
       pval=po(k,ist)
       typx=typ(ist)+100
       if((iprs.gt.1).and.(nvar.eq.1)) pval=pom(k,ist)
       if(iall.eq.1) then
c
       write(6,611) ist,indate,stnid,flon,flat,typx,typintx,elvx,hr,
     * k,pval,zo(k,ist),za(k,ist),zf(k,ist),zq(k,ist)
       write(6,612) ist,indate,stnid,flon,flat,typx,typintx,elvx,hr,
     * k,pval,to(k,ist),ta(k,ist),tf(k,ist),tq(k,ist)
       write(6,613) ist,indate,stnid,flon,flat,typx,typintx,elvx,hr,
     * k,pval,qo(k,ist),qa(k,ist),qf(k,ist),qq(k,ist)
c
       else
c
       write(6,711) ist,indate,stnid,flon,flat,typx,typintx,elvx,hr,
     * k,pval,zo(k,ist),zq(k,ist)
       write(6,712) ist,indate,stnid,flon,flat,typx,typintx,elvx,hr,
     * k,pval,to(k,ist),tq(k,ist)
       write(6,713) ist,indate,stnid,flon,flat,typx,typintx,elvx,hr,
     * k,pval,qo(k,ist),qq(k,ist)
c
       endif
       endif
c
       endif
c
c... now for wind report...
       if((kkx.eq.2).and.(iwnd.eq.1)) then
       kwmax(ist)=k
       kfillw(k,ist)=1
       mwfill(ist)=1
c
       CALL UFBINT(LUNIN,UOB, 4,255,NLEV,USTR)
       CALL UFBINT(LUNIN,VOB, 4,255,NLEV,VSTR)
c
       do i=1,4
       if(uob(i,l).ne.bmiss) upr(i,k)=uob(i,l)
       if(upr(i,k).gt.rmiss) upr(i,k)=bmiss
       if(vob(i,l).ne.bmiss) vpr(i,k)=vob(i,l)
       if(vpr(i,k).gt.rmiss) vpr(i,k)=bmiss
       enddo
c
       UO(K,IST) =  UPR(1,K)
       VO(K,IST) =  VPR(1,K)
       WQ(K,IST) =  UPR(4,K)
c
       if(iall.eq.1) then
       UA(K,IST) =  UPR(2,K)
       UF(K,IST) =  UPR(3,K)
       VA(K,IST) =  VPR(2,K)
       VF(K,IST) =  VPR(3,K)
       endif
c
       if(idbug.eq.1) then
       pval=po(k,ist)
       typx=typ(ist)+200
       if((iprs.gt.1).and.(nvar.eq.1)) pval=pow(k,ist)
       if(iall.eq.1) then
c
       write(6,614) ist,indate,stnid,flon,flat,typx,typintx,elvx,hr,
     * k,pval,uo(k,ist),ua(k,ist),uf(k,ist),wq(k,ist)
       write(6,615) ist,indate,stnid,flon,flat,typx,typintx,elvx,hr,
     * k,pval,vo(k,ist),va(k,ist),vf(k,ist)
c
       else
c
       write(6,714) ist,indate,stnid,flon,flat,typx,typintx,elvx,hr,
     * k,pval,uo(k,ist),wq(k,ist)
       write(6,715) ist,indate,stnid,flon,flat,typx,typintx,elvx,hr,
     * k,pval,vo(k,ist)
c
       endif
       endif
c
       endif
c
c...  end proper level section
       endif
c
c...  end quality mark section
c      endif
c
c... end level-loop
		   ENDDO
c
c... end lat/lon box section
       endif
c
c... end window-hour section
       endif
c
 443   continue
c
       ENDDO
c
c...  end report section
       endif
C
       ENDDO
c
 1555  continue
       nstns=IST
       print *,indate,' number of stations ',nstns
c
c.... get unique station ids...
c
       do j=1,nstns
        n=0
        idj=tstn(j)
        do i=j+1,nstns
        idi=tstn(i)
        if(idj.eq.idi) then
        n=n+1
        print *,' unique station ',idj,idi,j,i,n
        write(labn,'(i1)') n
        idi(6:6)=labn
        tstn(i)=idi
        endif
        enddo
       enddo
c
       nstf=0
       do ist=1,nstns
c
c...   find a common set of levels
c
       if((iprs.gt.1).and.(nvar.eq.1)) then
c
       im=1
       iw=1
       nl=0
       do k=1,nlevs
c
       zog(k)=bmiss
       tog(k)=bmiss
       qog(k)=bmiss
       zag(k)=bmiss
       tag(k)=bmiss
       qag(k)=bmiss
       zfg(k)=bmiss
       tfg(k)=bmiss
       qfg(k)=bmiss
       zqg(k)=bmiss
       tqg(k)=bmiss
       qqg(k)=bmiss
       uog(k)=bmiss
       vog(k)=bmiss
       uag(k)=bmiss
       vag(k)=bmiss
       ufg(k)=bmiss
       vfg(k)=bmiss
       wqg(k)=bmiss
c
       pm=pom(im,ist)
       pw=pow(iw,ist)
       if(pm.eq.bmiss.and.pw.eq.bmiss) go to 11
       if(pm.ne.bmiss.and.pw.eq.bmiss) then
       pog(k)=pm
       zog(k)=zo(im,ist)
       tog(k)=to(im,ist)
       qog(k)=qo(im,ist)
       zag(k)=za(im,ist)
       tag(k)=ta(im,ist)
       qag(k)=qa(im,ist)
       zfg(k)=zf(im,ist)
       tfg(k)=tf(im,ist)
       qfg(k)=qf(im,ist)
       zqg(k)=zq(im,ist)
       tqg(k)=tq(im,ist)
       qqg(k)=qq(im,ist)
       im=im+1
       nl=nl+1
       go to 11
       endif
       if(pm.eq.bmiss.and.pw.ne.bmiss) then
       pog(k)=pw
       uog(k)=uo(iw,ist)
       vog(k)=vo(iw,ist)
       uag(k)=ua(iw,ist)
       vag(k)=va(iw,ist)
       ufg(k)=uf(iw,ist)
       vfg(k)=vf(iw,ist)
       wqg(k)=wq(iw,ist)
       iw=iw+1
       nl=nl+1
       go to 11
       endif
c
       nl=nl+1
       if(pm.gt.pw) then
       pog(k)=pm
       zog(k)=zo(im,ist)
       tog(k)=to(im,ist)
       qog(k)=qo(im,ist)
       zag(k)=za(im,ist)
       tag(k)=ta(im,ist)
       qag(k)=qa(im,ist)
       zfg(k)=zf(im,ist)
       tfg(k)=tf(im,ist)
       qfg(k)=qf(im,ist)
       zqg(k)=zq(im,ist)
       tqg(k)=tq(im,ist)
       qqg(k)=qq(im,ist)
       im=im+1
       endif
       if(pm.eq.pw) then
       pog(k)=pm
       zog(k)=zo(im,ist)
       tog(k)=to(im,ist)
       qog(k)=qo(im,ist)
       zag(k)=za(im,ist)
       tag(k)=ta(im,ist)
       qag(k)=qa(im,ist)
       zfg(k)=zf(im,ist)
       tfg(k)=tf(im,ist)
       qfg(k)=qf(im,ist)
       zqg(k)=zq(im,ist)
       tqg(k)=tq(im,ist)
       qqg(k)=qq(im,ist)
       uog(k)=uo(iw,ist)
       vog(k)=vo(iw,ist)
       uag(k)=ua(iw,ist)
       vag(k)=va(iw,ist)
       ufg(k)=uf(iw,ist)
       vfg(k)=vf(iw,ist)
       wqg(k)=wq(iw,ist)
       im=im+1
       iw=iw+1
       endif
       if(pm.lt.pw) then
       pog(k)=pw
       uog(k)=uo(iw,ist)
       vog(k)=vo(iw,ist)
       uag(k)=ua(iw,ist)
       vag(k)=va(iw,ist)
       ufg(k)=uf(iw,ist)
       vfg(k)=vf(iw,ist)
       wqg(k)=wq(iw,ist)
       iw=iw+1
       endif
 11    continue
       enddo
c
       else
c
c... check how many levels are there...
       nl=0
       do l=1,nlevs
       pdata=po(l,ist)
       if(ist.eq.1) print *,l,pdata
       if(pdata.ne.bmiss) then
       nl=nl+1
       pog(nl)=po(l,ist)
       if(imas.eq.1) then
       zog(nl)=zo(l,ist)
       tog(nl)=to(l,ist)
       qog(nl)=qo(l,ist)
       zag(nl)=za(l,ist)
       tag(nl)=ta(l,ist)
       qag(nl)=qa(l,ist)
       zfg(nl)=zf(l,ist)
       tfg(nl)=tf(l,ist)
       qfg(nl)=qf(l,ist)
       zqg(nl)=zq(l,ist)
       tqg(nl)=tq(l,ist)
       qqg(nl)=qq(l,ist)
       endif
       if(iwnd.eq.1) then
       uog(nl)=uo(l,ist)
       vog(nl)=vo(l,ist)
       uag(nl)=ua(l,ist)
       vag(nl)=va(l,ist)
       ufg(nl)=uf(l,ist)
       vfg(nl)=vf(l,ist)
       wqg(nl)=wq(l,ist)
       endif
       endif
       enddo
c
       endif
c
       print *,indate,' number of levels for station ',ist,' ',nl
       if(nl.eq.0) go to 1234
       nstf=nstf+1
c
c....  write grads data out...
c
       id(1:8)='        '
       id=tstn(ist)
       rlon=tlon(ist)
       rlat=tlat(ist)
       hr=thr(ist)
       typx=typ(ist)
       typintx=typint(ist)
       elvx=elvt(ist)
       psog=pso(ist)
       psag=psa(ist)
       psfg=psf(ist)
       psqg=psq(ist)
       iflag=0
       id(7:7)=' '
       id(8:8)='\0'
c
       rt=0.0
       write(51) id,rlat,rlon,rt,nl,iflag
c
       if(iall.eq.1) then
c
       if((imas.eq.1).and.(iwnd.eq.1)) then
       write(51) (pog(l),
     *            typx,typintx,elvx,psog,psag,psfg,psqg,pog(l),
     *            zog(l),zag(l),zfg(l),zqg(l),
     *            tog(l),tag(l),tfg(l),tqg(l),
     *            qog(l),qag(l),qfg(l),qqg(l),
     *            uog(l),uag(l),ufg(l),wqg(l),
     *            vog(l),vag(l),vfg(l),l=1,nl)

       if(idbuga.eq.1) then
       print *,id,rlat,rlon,rt,nl,iflag
       do l=1,nl
       write(6,411) indate,hr,psog,psag,psfg,psqg,pog(l),
     *              zog(l),zag(l),zfg(l),zqg(l),
     *              tog(l),tag(l),tfg(l),tqg(l),
     *              qog(l),qag(l),qfg(l),qqg(l),
     *              uog(l),uag(l),ufg(l),wqg(l),
     *              vog(l),vag(l),vfg(l)
       enddo
       endif
       endif
       if((imas.eq.1).and.(iwnd.eq.0)) then
       write(51) (pog(l),
     *            typx,typintx,elvx,psog,psag,psfg,psqg,pog(l),
     *            zog(l),zag(l),zfg(l),zqg(l),
     *            tog(l),tag(l),tfg(l),tqg(l),
     *            qog(l),qag(l),qfg(l),qqg(l),l=1,nl)
       if(idbuga.eq.1) then
       print *,id,rlat,rlon,rt,nl,iflag
       do l=1,nl
       write(6,412) indate,hr,psog,psag,psfg,psqg,pog(l),
     *              zog(l),zag(l),zfg(l),zqg(l),
     *              tog(l),tag(l),tfg(l),tqg(l),
     *              qog(l),qag(l),qfg(l),qqg(l)
       enddo
       endif
       endif
       if((imas.eq.0).and.(iwnd.eq.1)) then
       write(51) (pog(l),
     *            typx,typintx,elvx,psog,psag,psfg,psqg,pog(l),
     *            uog(l),uag(l),ufg(l),wqg(l),
     *            vog(l),vag(l),vfg(l),l=1,nl)
       if(idbuga.eq.1) then
       print *,id,rlat,rlon,rt,nl,iflag
       do l=1,nl
       write(6,413) indate,hr,psog,psag,psfg,psqg,pog(l),
     *              uog(l),uag(l),ufg(l),wqg(l),
     *              vog(l),vag(l),vfg(l)
       enddo
       endif
       endif
c
       else
c
       if((imas.eq.1).and.(iwnd.eq.1)) then
       write(51) (pog(l),
     *            typx,typintx,elvx,psog,psqg,pog(l),
     *            zog(l),zqg(l),
     *            tog(l),tqg(l),
     *            qog(l),qqg(l),
     *            uog(l),wqg(l),
     *            vog(l),l=1,nl)
       if(idbuga.eq.1) then
       print *,id,rlat,rlon,rt,nl,iflag
       do l=1,nl
       write(6,511) indate,hr,psog,psqg,pog(l),
     *              zog(l),zqg(l),
     *              tog(l),tqg(l),
     *              qog(l),qqg(l),
     *              uog(l),wqg(l),
     *              vog(l)
       enddo
       endif
       endif
       if((imas.eq.1).and.(iwnd.eq.0)) then
       write(51) (pog(l),
     *            typx,typintx,elvx,psog,psqg,pog(l),
     *            zog(l),zqg(l),
     *            tog(l),tqg(l),
     *            qog(l),qqg(l),l=1,nl)
       if(idbuga.eq.1) then
       print *,id,rlat,rlon,rt,nl,iflag
       do l=1,nl
       write(6,512) indate,hr,psog,psqg,pog(l),
     *              zog(l),zqg(l),
     *              tog(l),tqg(l),
     *              qog(l),qqg(l)
       enddo
       endif
       endif
       if((imas.eq.0).and.(iwnd.eq.1)) then
       write(51) (pog(l),
     *            typx,typintx,elvx,psog,psqg,pog(l),
     *            uog(l),wqg(l),
     *            vog(l),l=1,nl)
       if(idbuga.eq.1) then
       print *,id,rlat,rlon,rt,nl,iflag
       do l=1,nl
       write(6,513) indate,hr,psog,psqg,pog(l),
     *              uog(l),wqg(l),
     *              vog(l)
       enddo
       endif
       endif
c
       endif
c
 1234  continue
c... end station-loop
        enddo
c
       print *,indate,' total number of stations for grads ',nstf
c
c..  write footer to end this time level
       id='       '
       rlat=0.
       rlon=0.
       rt=0.
       nlev=0
       iflag=0
       write(51) id,rlat,rlon,rt,nlev,iflag
c
       CALL CLOSBF(LUNIN)
       close(51)
c
 411   format(a10,1x,f5.2,1x,4f6.1,1x,f6.1,' mb ',1x,
     *        4(3e10.2,1x,f3.0,/),3e10.2)
 412   format(a10,1x,f5.2,1x,4f6.1,1x,f6.1,' mb ',1x,
     *        3(3e10.2,1x,f3.0,/))
 413   format(a10,1x,f5.2,1x,4f6.1,1x,f6.1,' mb ',1x,
     *        3e10.2,1x,f3.0,1x,3e10.2)
c
 511   format(a10,1x,f5.2,1x,2f6.1,1x,f6.1,' mb ',1x,
     *        2(e10.2,1x,f3.0),1x,e10.2)
 512   format(a10,1x,f5.2,1x,2f6.1,1x,f6.1,' mb ',1x,
     *        3(e10.2,1x,f3.0))
 513   format(a10,1x,f5.2,1x,2f6.1,1x,f6.1,' mb ',1x,
     *        e10.2,1x,f3.0,1x,e10.2)
c
 610   format(i6,1x,a10,1x,a8,1x,f6.2,1x,f6.2,1x,f4.0,1x,f4.0,1x,f6.1,
     * 1x,f5.2,1x,'ps',1x,3f10.2,1x,'qf ',f4.0)
 611   format(i6,1x,a10,1x,a8,1x,f6.2,1x,f6.2,1x,f4.0,1x,f4.0,1x,f6.1,
     * 1x,f5.2,1x,'z',1x,i2,1x,f6.1,1x,3f10.2,1x,'qf ',f4.0)
 612   format(i6,1x,a10,1x,a8,1x,f6.2,1x,f6.2,1x,f4.0,1x,f4.0,1x,f6.1,
     * 1x,f5.2,1x,'t',1x,i2,1x,f6.1,1x,3f10.2,1x,'qf ',f4.0)
 613   format(i6,1x,a10,1x,a8,1x,f6.2,1x,f6.2,1x,f4.0,1x,f4.0,1x,f6.1,
     * 1x,f5.2,1x,'q',1x,i2,1x,f6.1,1x,3f10.2,1x,'qf ',f4.0)
 614   format(i6,1x,a10,1x,a8,1x,f6.2,1x,f6.2,1x,f4.0,1x,f4.0,1x,f6.1,
     * 1x,f5.2,1x,'u',1x,i2,1x,f6.1,1x,3f10.2,1x,'qf ',f4.0)
 615   format(i6,1x,a10,1x,a8,1x,f6.2,1x,f6.2,1x,f4.0,1x,f4.0,1x,f6.1,
     * 1x,f5.2,1x,'v',1x,i2,1x,f6.1,1x,3f10.2)
c
 711   format(i6,1x,a10,1x,a8,1x,f6.2,1x,f6.2,1x,f4.0,1x,f4.0,1x,f6.1,
     * 1x,f5.2,1x,'z',1x,i2,1x,f6.1,1x,f10.2,1x,'qf ',f4.0)
 712   format(i6,1x,a10,1x,a8,1x,f6.2,1x,f6.2,1x,f4.0,1x,f4.0,1x,f6.1,
     * 1x,f5.2,1x,'t',1x,i2,1x,f6.1,1x,f10.2,1x,'qf ',f4.0)
 713   format(i6,1x,a10,1x,a8,1x,f6.2,1x,f6.2,1x,f4.0,1x,f4.0,1x,f6.1,
     * 1x,f5.2,1x,'q',1x,i2,1x,f6.1,1x,f10.2,1x,'qf ',f4.0)
 714   format(i6,1x,a10,1x,a8,1x,f6.2,1x,f6.2,1x,f4.0,1x,f4.0,1x,f6.1,
     * 1x,f5.2,1x,'u',1x,i2,1x,f6.1,1x,f10.2,1x,'qf ',f4.0)
 715   format(i6,1x,a10,1x,a8,1x,f6.2,1x,f6.2,1x,f4.0,1x,f4.0,1x,f6.1,
     * 1x,f5.2,1x,'v',1x,i2,1x,f6.1,1x,f10.2)
c
       return
       end
C-----------------------------------------------------------------------
       INTEGER FUNCTION nfill(C)
       CHARACTER*(*) C
       NFILL=LEN(C)
       DO 1 J=1,NFILL
       IF(C(J:J).EQ.' ') THEN
       NFILL=J-1
       RETURN
       ENDIF
 1     CONTINUE
       RETURN
       END
      FUNCTION GETLEV(P,plev,itol,nlevs,kfill,ist,nst)
c
      dimension plev(nlevs)
      integer kfill(nlevs,nst)
c
       IP = NINT(P)
c
c..  first check if exact match exists...
       do nl=1,nlevs
       getlev=0
       lsig=plev(nl)
c
       if((ip.eq.lsig).and.(kfill(nl,ist).eq.0)) then
       getlev=nl
c      print *,'from getlev exact ',nl,p,ip,kk,getlev
       return
       endif
c
c      print *,'from getlev exact ',nl,p,ip,kk,getlev
       enddo
c
       if(itol.eq.0) return
c
       do nl=1,nlevs
       getlev=0
       lsig=plev(nl)
c
       do kk=lsig-itol,lsig+itol
       if((ip.eq.kk).and.(kfill(nl,ist).eq.0)) then
       getlev=nl
c      print *,'from getlev tol ',nl,p,ip,kk,getlev
       return
       endif
       enddo
c
c      print *,'from getlev tol ',nl,p,ip,kk,getlev
       enddo
c
      RETURN
      END