         subroutine get_ob_lists(iter,miter,cgrid,igust,ivis,ipblh,idist)
!
!  to compile:
!    xlf90 get_ob_lists.f -L/nwprod/lib -lw3_8
!
!  abstract: 
!    read in from original unformatted gsi diagnostic files 
!    and write out formatted streamlined files. Valid for
!    gsi release 200609.
  

         implicit none

         integer(4),intent(in):: iter,miter
         integer(4),intent(in):: igust,ivis,ipblh,idist
         character(60),intent(in):: cgrid

         integer(4),parameter::nflds=6+4 !add gust,vis,pblh,dist

         integer(4),parameter::lun_t=10
         integer(4),parameter::lun_q=11
         integer(4),parameter::lun_ps=12
         integer(4),parameter::lun_u=13
         integer(4),parameter::lun_v=14
         integer(4),parameter::lun_spd=15
         integer(4),parameter::lun_gust=16
         integer(4),parameter::lun_vis=17
         integer(4),parameter::lun_pblh=18
         integer(4),parameter::lun_dist=19
         integer(4),parameter::nrejectmax=20000
         integer(4),parameter::nmplower=-10000
         integer(4),parameter::nmpupper=+100

         character(8),allocatable,dimension(:):: cdiagbuf
         character(8),allocatable,dimension(:):: cprvstg
         character(8),allocatable,dimension(:):: csprvstg
         character(3) otype
         character(8) cstation,cstationmax
         character(8) cprovider,csubprovider
         character(96) cheader
         character(3) clun3
         character(2) ctvts
         character(80) filename
         character(15) clistorig
         character(18) caux18
         character(15) cblank15
         character(1)  cblank1
         character(1)  chvar1
         character(10) chvar10
         character(90) cfmt

         integer(4) idate,nchar,nreal,i,ii,mype,lun,lun2,n,m
         integer(4) itype
         integer(4) ntrjs,nqrjs,nprjs,nwrjs

         integer(4) ntot(nflds),nmp(nmplower:nmpupper,nflds)
         integer(4) kuse,ifld
         integer(4) kk

         real(4),allocatable,dimension(:,:)::rdiagbuf
         real(4) rlat,rlon,oberr,oberr2,ob,ob_model,ddiff, & 
                 dudiff,dvdiff,rmuse,xx,yy
         real(4) uob,uob_model,vob,vob_model,rfactor,qtflg
         real(8) rlat8,rlon8,xx8,yy8
         real(8) da8,alat18,elon18,elonv8,alatan8
         real(4) shgt0  !station height
         real(4) hgt0   !observation elevation
         real(4) hgt    !model terrain at ob location
         real(4) slm    !dominant surface type
         real(4) dtime
         real(4) rtvts
         real(4) percnumber
         character(90) t_rjlist(nrejectmax)
         character(90) q_rjlist(nrejectmax)
         character(90) p_rjlist(nrejectmax)
         character(90) w_rjlist(nrejectmax)
         logical muse
         logical tlistexist,qlistexist,plistexist,wlistexist
         logical diagexist
         logical fexist, r15887_diagfile_fmt

         namelist/diagfile_fmt/r15887_diagfile_fmt

         data r15887_diagfile_fmt/.true./

         inquire(file='diagfile_fmt_input',exist=fexist)
         if (fexist) then
            open (55,file='diagfile_fmt_input',form='formatted')
            read(55,*) r15887_diagfile_fmt
            close(55)
         endif

         print*,'in get_ob_lists:: cgrid,iter,miter,r15887_diagfile_fmt=',trim(cgrid),iter,miter,r15887_diagfile_fmt

         if (iter .eq. 1) clun3='ges'

         if (iter .gt. 1) then
            if((iter-1).eq.miter) then
               clun3='anl'
              else
               clun3='   '
               write(clun3(1:2),'(i2.2)')iter
            endif
         endif
         print*,'in get_ob_lists:: clun3=',trim(clun3)


         diagexist=.false.
         inquire(file='diag_conv_'//trim(clun3)//'.dat',exist=diagexist)
         if (diagexist) open (7,file='diag_conv_'//trim(clun3)//'.dat',form='unformatted')
         if (.not.diagexist) return

         if (iter .eq. 1) clun3='01'   !keep old filenaming

         open(lun_t,file='t_obs.listing_iter_'//trim(clun3),form='formatted') !output file

         open(lun_q,file='q_obs.listing_iter_'//trim(clun3),form='formatted') !output file

         open(lun_ps,file='ps_obs.listing_iter_'//trim(clun3),form='formatted') !output file

         open(lun_u,file='u_obs.listing_iter_'//trim(clun3),form='formatted') !output file
 
         open(lun_v,file='v_obs.listing_iter_'//trim(clun3),form='formatted') !output file

         open(lun_spd,file='spd_obs.listing_iter_'//trim(clun3),form='formatted') !output file

         if (igust > 0) open(lun_gust,file='gust_obs.listing_iter_'//trim(clun3),form='formatted') !output file

         if (ivis  > 0) open(lun_vis,file='vis_obs.listing_iter_'//trim(clun3),form='formatted')   !output file

         if (ipblh > 0) open(lun_pblh,file='pblh_obs.listing_iter_'//trim(clun3),form='formatted') !output file

         if (idist > 0) open(lun_dist,file='dist_obs.listing_iter_'//trim(clun3),form='formatted') !output file

         cheader='stnname   obtype   lat(dg)  lon(dg E)   dtime      oberr         ob         guess         rmuse'
         write(lun_t,'(a)')   'RTMA TEMPERATURE OBS' 
         write(lun_t,'(a)')   'UNITS of  oberr, ob, and guess are K. Ob is used only if rmuse=+1.0 or +2.0' 
         write(lun_t,'(a)')   'rmuse=-100. ==> user chose to monitor this ob and see how well it agrees with the guess' 
         write(lun_t,'(a)')   'rmuse=-150. ==> ob is being monitored. Internally selected based on MADIS QC flag values'
         write(lun_t,'(a)')   'rmuse=-5000. ==> ob was in the reject list ; rmuse=-5100. ==> in diurnal list' 
         write(lun_t,'(a)')   'dtime is the hour relative to the valid analysis time. For example, dtime=-0.1' 
         write(lun_t,'(a)')   '       means 0.1h (i.e. 6 minutes) before the valid analysis time'
         write(lun_t,'(a)')   'rejectlist: list of sub-standard obs where ob was found. It can be (i) static (sta),'
         write(lun_t,'(a)')   '            (ii) from weather forecast office (wfo), (iii) global based on MADIS QC stats (glb), or' 
         write(lun_t,'(a)')   '            (iv) dynamic (dyn). Note that the ob can be in more than one reject list.'
         write(lun_t,'(a)')   'tv=virtual T, and ts=sensible T. Note that the GSI can assimilate obs as either tv or ts' 
         write(lun_t,'(a)')   '                              '
         write(lun_t,'(a)')   cheader//' rejectlist  tv-or-ts'

         write(lun_q,'(a)')   'RTMA SPECIFIC HUMIDITY OBS' 
         write(lun_q,'(a)')   'UNITS of  oberr, ob, and guess are g/Kg. Ob is used only if rmuse=+1.0 or +2.0' 
         write(lun_q,'(a)')   'rmuse=-100. ==> user chose to monitor this ob and see how well it agrees with the guess' 
         write(lun_q,'(a)')   'rmuse=-150. ==> ob is being monitored. Internally selected based on MADIS QC flag values'
         write(lun_q,'(a)')   'rmuse=-5000. ==> ob was in the reject list ; rmuse=-5100. ==> in diurnal list'
         write(lun_q,'(a)')   'dtime is the hour relative to the valid analysis time. For example, dtime=-0.1' 
         write(lun_q,'(a)')   '       means 0.1h (i.e. 6 minutes) before the valid analysis time'
         write(lun_q,'(a)')   'rejectlist: lists of sub-standard obs where ob was found. It can be (i) static (sta),'
         write(lun_q,'(a)')   '            (ii) from weather forecast office (wfo), (iii) global based on MADIS QC stats (glb), or' 
         write(lun_q,'(a)')   '            (iv) dynamic (dyn). Note that the ob can be in more than one reject list'
         write(lun_q,'(a)')   '                              '
         write(lun_q,'(a)')   cheader//' rejectlist'

         write(lun_ps,'(a)')   'RTMA SURFACE PRESSURE OBS' 
         write(lun_ps,'(a)')   'UNITS of  oberr, ob, and guess are Pa. Ob is used only if rmuse=+1.0 or +2.0' 
         write(lun_ps,'(a)')   'rmuse=-100. ==> user chose to monitor this ob and see how well it agrees with the guess' 
         write(lun_ps,'(a)')   'rmuse=-150. ==> ob is being monitored. Internally selected based on MADIS QC flag values'
         write(lun_ps,'(a)')   'rmuse=-5000. ==> ob was in the reject list'
         write(lun_ps,'(a)')   'dtime is the hour relative to the valid analysis time. For example, dtime=-0.1' 
         write(lun_ps,'(a)')   '       means 0.1h (i.e. 6 minutes) before the valid analysis time'
         write(lun_ps,'(a)')   'rejectlist: list of sub-standard obs where ob was found. It can be (i) static (sta),'
         write(lun_ps,'(a)')   '            (ii) from weather forecast office (wfo), (iii) global based on MADIS QC stats (glb), or' 
         write(lun_ps,'(a)')   '            (iv) dynamic (dyn). Note that the ob can be in more than one reject list'
         write(lun_ps,'(a)')   '                              '
         write(lun_ps,'(a)')   cheader//' rejectlist'

         write(lun_u,'(a)')   'RTMA U-WIND OBS' 
         write(lun_u,'(a)')   'UNITS of  oberr, ob, and guess are m/s. Ob is used only if rmuse=+1.0 or +2.0' 
         write(lun_u,'(a)')   'rmuse=-100. ==> user chose to monitor this ob and see how well it agrees with the guess' 
         write(lun_u,'(a)')   'rmuse=-150. ==> ob is being monitored. Internally selected based on MADIS QC flag values'
         write(lun_u,'(a)')   'rmuse=-5000.==> this non-mesonet wind was was in the reject list' 
         write(lun_u,'(a)')   'rmuse=-6000.==> this mesonet wind did not belong to any of the GSD uselists, and neither was it'
         write(lun_u,'(a)')   '                in the reject list'
         write(lun_u,'(a)')   'rmuse=-6100.==> this mesonet wind was in the reject list and in at least one of the GSD uselists'
         write(lun_u,'(a)')   'rmuse=-6200.==> this mesonet wind was in the reject list and in neither one of the GSD uselists'
         write(lun_u,'(a)')   'dtime is the hour relative to the valid analysis time. For example, dtime=-0.1' 
         write(lun_u,'(a)')   '       means 0.1h (i.e. 6 minutes) before the valid analysis time'
         write(lun_u,'(a)')   'rejectlist: list of sub-standard obs where ob was found. It can be (i) static (sta),'
         write(lun_u,'(a)')   '            (ii) from weather forecast office (wfo), (iii) global based on MADIS QC stats (glb), or' 
         write(lun_u,'(a)')   '            (iv) dynamic (dyn). Note that the ob can be in more than one reject list'
         write(lun_u,'(a)')   '                              '
         write(lun_u,'(a)')   cheader//' rejectlist'

         write(lun_v,'(a)')   'RTMA V-WIND OBS' 
         write(lun_v,'(a)')   'UNITS of  oberr, ob, and guess are m/s. Ob is used only if rmuse=+1.0 or +2.0' 
         write(lun_v,'(a)')   'rmuse=-100. ==> user chose to monitor this ob and see how well it agrees with the guess' 
         write(lun_v,'(a)')   'rmuse=-150. ==> ob is being monitored. Internally selected based on MADIS QC flag values'
         write(lun_v,'(a)')   'rmuse=-5000.==> this non-mesonet wind was was in the reject list' 
         write(lun_v,'(a)')   'rmuse=-6000.==> this mesonet wind did not belong to any of the GSD uselists, and neither was it'
         write(lun_v,'(a)')   '                in the reject list'
         write(lun_v,'(a)')   'rmuse=-6100.==> this mesonet wind was in the reject list and on at least one of the GSD uselists'
         write(lun_v,'(a)')   'rmuse=-6200.==> this mesonet wind was in the reject list and on neither one of the GSD uselists'
         write(lun_v,'(a)')   'dtime is the hour relative to the valid analysis time. For example, dtime=-0.1' 
         write(lun_v,'(a)')   '       means 0.1h (i.e. 6 minutes) before the valid analysis time'
         write(lun_v,'(a)')   'rejectlist: list of sub-standard obs where ob was found. It can be (i) static (sta),'
         write(lun_v,'(a)')   '            (ii) from weather forecast office (wfo), (iii) global based on MADIS QC stats (glb), or' 
         write(lun_v,'(a)')   '            (iv) dynamic (dyn). Note that the ob can be in more than one reject list' 
         write(lun_v,'(a)')   '                              '
         write(lun_v,'(a)')   cheader//' rejectlist'

         write(lun_spd,'(a)')   'RTMA WIND SPEED OBS' 
         write(lun_spd,'(a)')   'UNITS of  oberr, ob, and guess are m/s. Ob is used only if rmuse=+1.0 or +2.0' 
         write(lun_spd,'(a)')   'rmuse=-100. ==> user chose to monitor this ob and see how well it agrees with the guess' 
         write(lun_spd,'(a)')   'rmuse=-150. ==> ob is being monitored. Internally selected based on MADIS QC flag values'
         write(lun_spd,'(a)')   'rmuse=-500. ==> ob was in the reject list'
         write(lun_spd,'(a)')   'dtime is the hour relative to the valid analysis time. For example, dtime=-0.1' 
         write(lun_spd,'(a)')   '       means 0.1h (i.e. 6 minutes) before the valid analysis time'
         write(lun_spd,'(a)')   'rejectlist: list of sub-standard obs where ob was found. It can be (i) static (sta),'
         write(lun_spd,'(a)')   '            (ii) from weather forecast office (wfo), (iii) global based on MADIS QC stats (glb), or' 
         write(lun_spd,'(a)')   '            (iv) dynamic (dyn). Note that the ob can be in more than one reject list'
         write(lun_spd,'(a)')   '                              '
         write(lun_spd,'(a)')   cheader//' rejectlist'

         if (igust > 0) then
             write(lun_gust,'(a)')   'RTMA GUST-WIND OBS' 
             write(lun_gust,'(a)')   'UNITS of  oberr, ob, and guess are m/s. Ob is used only if rmuse=+1.0 or +2.0' 
             write(lun_gust,'(a)')   'rmuse=-100. ==> user chose to monitor this ob and see how well it agrees with the guess' 
             write(lun_gust,'(a)')   'rmuse=-150. ==> ob is being monitored. Internally selected based on MADIS QC flag values'
             write(lun_gust,'(a)')   'rmuse=-5000.==> this non-mesonet wind was was in the reject list' 
             write(lun_gust,'(a)')   'rmuse=-6000.==> this mesonet wind did not belong to any of the GSD uselists, and neither was it'
             write(lun_gust,'(a)')   '                in the reject list'
             write(lun_gust,'(a)')   'rmuse=-6100.==> this mesonet wind was in the reject list and on at least one of the GSD uselists'
             write(lun_gust,'(a)')   'rmuse=-6200.==> this mesonet wind was in the reject list and on neither one of the GSD uselists'
             write(lun_gust,'(a)')   'dtime is the hour relative to the valid analysis time. For example, dtime=-0.1' 
             write(lun_gust,'(a)')   '       means 0.1h (i.e. 6 minutes) before the valid analysis time'
             write(lun_gust,'(a)')   'rejectlist: list of sub-standard obs where ob was found. It can be (i) static (sta),'
             write(lun_gust,'(a)')   '            (ii) from weather forecast office (wfo), (iii) global based on MADIS QC stats (glb), or' 
             write(lun_gust,'(a)')   '            (iv) dynamic (dyn). Note that the ob can be in more than one reject list' 
             write(lun_gust,'(a)')   '                              '
             write(lun_gust,'(a)')   cheader//' rejectlist'
         endif

         if (ivis > 0) then
             write(lun_vis,'(a)')   'RTMA VISIBILITY OBS' 
             write(lun_vis,'(a)')   'UNITS of  oberr, ob, and guess is m. Ob is used only if rmuse=+1.0 or +2.0' 
             write(lun_vis,'(a)')   'rmuse=-100. ==> user chose to monitor this ob and see how well it agrees with the guess' 
             write(lun_vis,'(a)')   'rmuse=-150. ==> ob is being monitored. Internally selected based on MADIS QC flag values'
             write(lun_vis,'(a)')   'rmuse=-5000.==> this non-mesonet wind was was in the reject list' 
             write(lun_vis,'(a)')   'rmuse=-6000.==> this mesonet wind did not belong to any of the GSD uselists, and neither was it'
             write(lun_vis,'(a)')   '                in the reject list'
             write(lun_vis,'(a)')   'rmuse=-6100.==> this mesonet wind was in the reject list and on at least one of the GSD uselists'
             write(lun_vis,'(a)')   'rmuse=-6200.==> this mesonet wind was in the reject list and on neither one of the GSD uselists'
             write(lun_vis,'(a)')   'dtime is the hour relative to the valid analysis time. For example, dtime=-0.1' 
             write(lun_vis,'(a)')   '       means 0.1h (i.e. 6 minutes) before the valid analysis time'
             write(lun_vis,'(a)')   'rejectlist: list of sub-standard obs where ob was found. It can be (i) static (sta),'
             write(lun_vis,'(a)')   '            (ii) from weather forecast office (wfo), (iii) global based on MADIS QC stats (glb), or' 
             write(lun_vis,'(a)')   '            (iv) dynamic (dyn). Note that the ob can be in more than one reject list' 
             write(lun_vis,'(a)')   '                              '
             write(lun_vis,'(a)')   cheader//' rejectlist'
         endif

         if (ipblh > 0) then
             write(lun_pblh,'(a)')   'RTMA PBLH OBS' 
             write(lun_pblh,'(a)')   'UNITS of  oberr, ob, and guess is m. Ob is used only if rmuse=+1.0 or +2.0' 
             write(lun_pblh,'(a)')   'rmuse=-100. ==> user chose to monitor this ob and see how well it agrees with the guess' 
             write(lun_pblh,'(a)')   'rmuse=-150. ==> ob is being monitored. Internally selected based on MADIS QC flag values'
             write(lun_pblh,'(a)')   'rmuse=-5000.==> this non-mesonet wind was was in the reject list' 
             write(lun_pblh,'(a)')   'rmuse=-6000.==> this mesonet wind did not belong to any of the GSD uselists, and neither was it'
             write(lun_pblh,'(a)')   '                in the reject list'
             write(lun_pblh,'(a)')   'rmuse=-6100.==> this mesonet wind was in the reject list and on at least one of the GSD uselists'
             write(lun_pblh,'(a)')   'rmuse=-6200.==> this mesonet wind was in the reject list and on neither one of the GSD uselists'
             write(lun_pblh,'(a)')   'dtime is the hour relative to the valid analysis time. For example, dtime=-0.1' 
             write(lun_pblh,'(a)')   '       means 0.1h (i.e. 6 minutes) before the valid analysis time'
             write(lun_pblh,'(a)')   'rejectlist: list of sub-standard obs where ob was found. It can be (i) static (sta),'
             write(lun_pblh,'(a)')   '            (ii) from weather forecast office (wfo), (iii) global based on MADIS QC stats (glb), or' 
             write(lun_pblh,'(a)')   '            (iv) dynamic (dyn). Note that the ob can be in more than one reject list' 
             write(lun_pblh,'(a)')   '                              '
             write(lun_pblh,'(a)')   cheader//' rejectlist'
         endif

         if (idist > 0) then
             write(lun_dist,'(a)')   'RTMA CEELING HEIGHT OBS' 
             write(lun_dist,'(a)')   'UNITS of  oberr, ob, and guess is m. Ob is used only if rmuse=+1.0 or +2.0' 
             write(lun_dist,'(a)')   'rmuse=-100. ==> user chose to monitor this ob and see how well it agrees with the guess' 
             write(lun_dist,'(a)')   'rmuse=-150. ==> ob is being monitored. Internally selected based on MADIS QC flag values'
             write(lun_dist,'(a)')   'rmuse=-5000.==> this non-mesonet wind was was in the reject list' 
             write(lun_dist,'(a)')   'rmuse=-6000.==> this mesonet wind did not belong to any of the GSD uselists, and neither was it'
             write(lun_dist,'(a)')   '                in the reject list'
             write(lun_dist,'(a)')   'rmuse=-6100.==> this mesonet wind was in the reject list and on at least one of the GSD uselists'
             write(lun_dist,'(a)')   'rmuse=-6200.==> this mesonet wind was in the reject list and on neither one of the GSD uselists'
             write(lun_dist,'(a)')   'dtime is the hour relative to the valid analysis time. For example, dtime=-0.1' 
             write(lun_dist,'(a)')   '       means 0.1h (i.e. 6 minutes) before the valid analysis time'
             write(lun_dist,'(a)')   'rejectlist: list of sub-standard obs where ob was found. It can be (i) static (sta),'
             write(lun_dist,'(a)')   '            (ii) from weather forecast office (wfo), (iii) global based on MADIS QC stats (glb), or' 
             write(lun_dist,'(a)')   '            (iv) dynamic (dyn). Note that the ob can be in more than one reject list' 
             write(lun_dist,'(a)')   '                              '
             write(lun_dist,'(a)')   cheader//' rejectlist'
         endif

        
         if (iter==1) call init_auxlists(clun3,igust,ivis,ipblh,idist)


         filename='t_rejectlist'
         call readin_rjlist(filename,t_rjlist,nrejectmax,tlistexist,ntrjs)
         print*,'in get_ob_lists, tlistexist,ntrjs=',tlistexist,ntrjs

         filename='q_rejectlist'
         call readin_rjlist(filename,q_rjlist,nrejectmax,qlistexist,nqrjs)
         print*,'in get_ob_lists, qlistexist,nqrjs=',qlistexist,nqrjs

         filename='p_rejectlist'
         call readin_rjlist(filename,p_rjlist,nrejectmax,plistexist,nprjs)
         print*,'in get_ob_lists, plistexist,nprjs=',plistexist,nprjs

         filename='w_rejectlist'
         call readin_rjlist(filename,w_rjlist,nrejectmax,wlistexist,nwrjs)
         print*,'in get_ob_lists, wlistexist,nwrjs=',wlistexist,nwrjs

         cblank1=' '
         cblank15='               '


         ntot=0
         nmp(nmplower:nmpupper,1:nflds)=0

         call proj_info(cgrid,da8,alat18,elon18,elonv8,alatan8)
         print*,'in get_ob_lists: da8,alat18,elon18,elonv8,alatan8=', &
                 da8,alat18,elon18,elonv8,alatan8

         read(7,end=200) idate

100      continue
         read(7,end=200)otype,nchar,nreal,ii,mype

         print*,'in get_ob_lists, idate=',idate
         print*,'in get_ob_lists, otype,nchar,nreal,ii,mype=', & 
                    otype,nchar,nreal,ii,mype
  
         if (ii==0) then 
            read(7,end=200)
            if (r15887_diagfile_fmt) read(7,end=200)
            goto 100
         endif

         allocate(cdiagbuf(ii))
         allocate(cprvstg(ii))
         allocate(csprvstg(ii))
         allocate(rdiagbuf(nreal,ii))

         if (r15887_diagfile_fmt) then
             read(7) cdiagbuf,rdiagbuf
             read(7) cprvstg,csprvstg
           else
             read(7) cdiagbuf,rdiagbuf,cprvstg,csprvstg
         endif

        if (otype(1:3)=='  t') then
          lun=lun_t
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then 
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           clistorig=cblank15
           if (rmuse < -1.) call read_rejfileorig(itype,cstation,t_rjlist,nrejectmax,ntrjs,clistorig)
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           slm=rdiagbuf(20,i)
           hgt=rdiagbuf(21,i)
           ctvts='tv'
           if (qtflg .gt. 0.) ctvts='ts'
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8
           oberr2=1.e10
           if (oberr.gt.1.e-05) oberr2=1./oberr  
!          write(lun,124) cstation,itype,rlat,rlon,dtime,oberr2,ob,ob_model,rmuse,clistorig,ctvts
           caux18=clistorig//cblank1//ctvts   !use this to get around formatting pbs in write statement
           write(lun,124) cstation,itype,rlat,rlon,dtime,oberr2,ob,ob_model,rmuse,caux18

           if (iter==1) call add2auxlist(otype,cstation,cprovider,csubprovider,itype,shgt0,hgt0,hgt,slm)

           ifld=1
           ntot(ifld)=ntot(ifld)+1
           if (rmuse .lt. 0.) kuse=nint(rmuse-0.001)
           if (rmuse .ge. 0.) kuse=nint(rmuse+0.001)
           nmp(kuse,ifld)=nmp(kuse,ifld)+1
          enddo
        end if


        if (otype(1:3)=='  q') then
          lun=lun_q
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           clistorig=cblank15
           if (rmuse < -1.) call read_rejfileorig(itype,cstation,q_rjlist,nrejectmax,nqrjs,clistorig)
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           slm=rdiagbuf(21,i)
           hgt=rdiagbuf(22,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8
           oberr2=1.e10
           if (oberr.gt.1.e-05) oberr2=1./oberr  
           write(lun,123) cstation,itype,rlat,rlon,dtime,oberr2,ob,ob_model,rmuse,clistorig

           if (iter==1) call add2auxlist(otype,cstation,cprovider,csubprovider,itype,shgt0,hgt0,hgt,slm)

           ifld=2
           ntot(ifld)=ntot(ifld)+1
           if (rmuse .lt. 0.) kuse=nint(rmuse-0.001)
           if (rmuse .ge. 0.) kuse=nint(rmuse+0.001)
           nmp(kuse,ifld)=nmp(kuse,ifld)+1
          enddo
        end if

        if (otype(2:3)=='ps') then
          lun=lun_ps
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           clistorig=cblank15
           if (rmuse < -1.) call read_rejfileorig(itype,cstation,p_rjlist,nrejectmax,nprjs,clistorig)
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           slm=rdiagbuf(20,i)
           hgt=rdiagbuf(21,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8
           oberr2=1.e10
           if (oberr.gt.1.e-05) oberr2=(1./oberr)*100. !convert to Pa
           ob=ob*100. !convert to Pa
           ob_model=ob_model*100.! convert to Pa
           write(lun,123) cstation,itype,rlat,rlon,dtime,oberr2,ob,ob_model,rmuse,clistorig

           if (iter==1) call add2auxlist(otype,cstation,cprovider,csubprovider,itype,shgt0,hgt0,hgt,slm)

           ifld=3
           ntot(ifld)=ntot(ifld)+1
           if (rmuse .lt. 0.) kuse=nint(rmuse-0.001)
           if (rmuse .ge. 0.) kuse=nint(rmuse+0.001)
           nmp(kuse,ifld)=nmp(kuse,ifld)+1
          enddo
        end if


        if (otype(2:3)=='uv') then
          lun=lun_u
          lun2=lun_v
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           clistorig=cblank15
           if (rmuse < -1.) call read_rejfileorig(itype,cstation,w_rjlist,nrejectmax,nwrjs,clistorig)
           oberr=rdiagbuf(16,i)
           uob=rdiagbuf(17,i)
           dudiff=rdiagbuf(18,i)
           uob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           vob=rdiagbuf(20,i)
           dvdiff=rdiagbuf(21,i)
           vob_model=rdiagbuf(20,i)-rdiagbuf(21,i)
           slm=rdiagbuf(24,i)
           hgt=rdiagbuf(25,i)
           rfactor=rdiagbuf(23,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8
           oberr2=1.e10
           if (oberr.gt.1.e-05) oberr2=1./oberr  
           write(lun,123) cstation,itype,rlat,rlon,dtime,oberr2,uob,uob_model,rmuse,clistorig
           write(lun2,123) cstation,itype,rlat,rlon,dtime,oberr2,vob,vob_model,rmuse,clistorig

           if (iter==1) call add2auxlist(otype,cstation,cprovider,csubprovider,itype,shgt0,hgt0,hgt,slm)

           ifld=4
           ntot(ifld)=ntot(ifld)+1
           if (rmuse .lt. 0.) kuse=nint(rmuse-0.001)
           if (rmuse .ge. 0.) kuse=nint(rmuse+0.001)
           nmp(kuse,ifld)=nmp(kuse,ifld)+1
          enddo
        end if

        if (otype(1:3)=='spd') then
          lun=lun_spd
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           clistorig=cblank15
           if (rmuse < -1.) call read_rejfileorig(itype,cstation,w_rjlist,nrejectmax,nwrjs,clistorig)
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           rfactor=rdiagbuf(20,i)
           slm=rdiagbuf(21,i)
           hgt=rdiagbuf(22,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8
           oberr2=1.e10
           if (oberr.gt.1.e-05) oberr2=1./oberr  
           write(lun,123) cstation,itype,rlat,rlon,dtime,oberr2,ob,ob_model,rmuse,clistorig

           if (iter==1) call add2auxlist(otype,cstation,cprovider,csubprovider,itype,shgt0,hgt0,hgt,slm)

           ifld=6
           ntot(ifld)=ntot(ifld)+1
           if (rmuse .lt. 0.) kuse=nint(rmuse-0.001)
           if (rmuse .ge. 0.) kuse=nint(rmuse+0.001)
           nmp(kuse,ifld)=nmp(kuse,ifld)+1
          enddo
        end if

        if (otype(1:3)=='gst' .and. igust >0) then
          lun=lun_gust
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           clistorig=cblank15
           if (rmuse < -1.) call read_rejfileorig(itype,cstation,w_rjlist,nrejectmax,nwrjs,clistorig)
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           slm=rdiagbuf(21,i)
           hgt=rdiagbuf(22,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8
           oberr2=1.e10
           if (oberr.gt.1.e-05) oberr2=1./oberr  
           write(lun,123) cstation,itype,rlat,rlon,dtime,oberr2,ob,ob_model,rmuse,clistorig

           if (iter==1) call add2auxlist(otype,cstation,cprovider,csubprovider,itype,shgt0,hgt0,hgt,slm)

           ifld=7
           ntot(ifld)=ntot(ifld)+1
           if (rmuse .lt. 0.) kuse=nint(rmuse-0.001)
           if (rmuse .ge. 0.) kuse=nint(rmuse+0.001)
           nmp(kuse,ifld)=nmp(kuse,ifld)+1
          enddo
        end if

        if (otype(1:3)=='vis' .and. ivis >0) then
          lun=lun_vis
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           clistorig=cblank15
!          if (rmuse < -1.) call read_rejfileorig(itype,cstation,vis_rjlist,nrejectmax,nwrjs,clistorig)
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           slm=rdiagbuf(21,i)
           hgt=rdiagbuf(22,i)
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8
           oberr2=1.e10
           if (oberr.gt.1.e-05) oberr2=1./oberr  
           write(lun,123) cstation,itype,rlat,rlon,dtime,oberr2,ob,ob_model,rmuse,clistorig

           if (iter==1) call add2auxlist(otype,cstation,cprovider,csubprovider,itype,shgt0,hgt0,hgt,slm)

           ifld=8
           ntot(ifld)=ntot(ifld)+1
           if (rmuse .lt. 0.) kuse=nint(rmuse-0.001)
           if (rmuse .ge. 0.) kuse=nint(rmuse+0.001)
           nmp(kuse,ifld)=nmp(kuse,ifld)+1
          enddo
        end if

        if (otype(1:3)=='pbl' .and. ipblh >0) then
          lun=lun_pblh
          do i=1,ii
           cstation=cdiagbuf(i)
           cprovider=cprvstg(i)
           csubprovider=csprvstg(i)
           itype=nint(rdiagbuf(1,i))
           rlat=rdiagbuf(3,i)
           rlon=rdiagbuf(4,i)
           dtime=rdiagbuf(8,i)
           shgt0=rdiagbuf(5,i)
           hgt0=rdiagbuf(7,i)
           qtflg=rdiagbuf(10,i)
           if (rdiagbuf(11,i) .gt. 1.) then
               rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
             else
               rmuse=rdiagbuf(12,i)
           endif
           clistorig=cblank15
!          if (rmuse < -1.) call read_rejfileorig(itype,cstation,pblh_rjlist,nrejectmax,nwrjs,clistorig)
           oberr=rdiagbuf(16,i)
           ob=rdiagbuf(17,i)
           ddiff=rdiagbuf(18,i)
           ob_model=rdiagbuf(17,i)-rdiagbuf(18,i)
           slm=-999.00000
           hgt=-999.00000
           rlat8=rlat*1._8
           rlon8=rlon*1._8
           call latlon_to_grid(da8,alat18,elon18,elonv8,alatan8, &
                               rlat8,rlon8,cgrid,xx8,yy8)
           xx=xx8
           yy=yy8
           oberr2=1.e10
           if (oberr.gt.1.e-05) oberr2=1./oberr  
           write(lun,123) cstation,itype,rlat,rlon,dtime,oberr2,ob,ob_model,rmuse,clistorig

           if (iter==1) call add2auxlist(otype,cstation,cprovider,csubprovider,itype,shgt0,hgt0,hgt,slm)

           ifld=9
           ntot(ifld)=ntot(ifld)+1
           if (rmuse .lt. 0.) kuse=nint(rmuse-0.001)
           if (rmuse .ge. 0.) kuse=nint(rmuse+0.001)
           nmp(kuse,ifld)=nmp(kuse,ifld)+1
          enddo
        end if

!       add dist stuff here

        deallocate(cdiagbuf,rdiagbuf)
        deallocate(cprvstg,csprvstg)

        goto 100
200     continue

        !v-wind component
        ntot(5)=ntot(4)
        nmp(:,5)=nmp(:,4)

        do n=1,nflds
           if (n==7  .and. igust <= 0) cycle
           if (n==8  .and. ivis  <= 0) cycle
           if (n==9  .and. ipblh <= 0) cycle
           if (n==10 .and. idist <= 0) cycle

           if (n==1) lun=lun_t
           if (n==2) lun=lun_q
           if (n==3) lun=lun_ps
           if (n==4) lun=lun_u
           if (n==5) lun=lun_v
           if (n==6) lun=lun_spd
           if (n==7) lun=lun_gust
           if (n==8) lun=lun_vis
           if (n==9) lun=lun_pblh
           if (n==10) lun=lun_dist


           write (lun,'(a)') '=================================================================================================='
           write (lun,'(a)') '=================================================================================================='
           write (lun,'(1x,a,i8)') 'total number of obs                            :', ntot(n)
           write (lun,'(1x,a,i8)') 'number of obs assimilated                      :', nmp(1,n)+nmp(2,n)+nmp(3,n)
           percnumber=0.
           if (ntot(n) .gt. 0) percnumber=float((nmp(1,n)+nmp(2,n)+nmp(3,n)))/float(ntot(n))*100.
           write (lun,'(1x,a,3x,f8.5,a)') 'percentage of obs assimilated                  :',percnumber,'%'
!          write (lun,'(1x,a,i8)') 'number of obs rejected by gross error check    :', nmp(-1,n)
!          percnumber=0.
!          if (ntot(n) .gt. 0) percnumber=float(nmp(-1,n))/float(ntot(n))*100.
!          write (lun,'(1x,a,3x,f8.5,a)') 'percentage of obs rejected by gross error check:',percnumber,'%'
           write (lun,'(1x,a,i8)') 'number of obs with rmuse=-1.                   :', nmp(-1,n)
           write (lun,'(1x,a,i8)') 'number of obs with rmuse=-2.                   :', nmp(-2,n)
           write (lun,'(1x,a,i8,2x,a)') 'number of obs with rmuse=-3.                   :', nmp(-3,n), 'must be zero if # of outer-iterations is only 2'
           do m=-499,-4
              if ( nmp(m,n) .gt. 0) then
                  write(chvar10,'(i10)') nmp(m,n)
                  chvar10=adjustl(chvar10)
                  kk=len_trim(chvar10)
!                 write(chvar1,"(i1)") kk+1
                  chvar1='8'
                  cfmt="(1x,"//"'number of obs with rmuse='"//',f5.0,'//"'                 :'"//',i'//chvar1//")"
                  write (lun,cfmt) float(m),nmp(m,n)
              endif
           enddo

           if (n.ne.4 .and. n.ne.5)  then
              write (lun,'(1x,a,i8)') 'number of obs in the reject list               :', nmp(-5000,n)
           endif
           if (n.eq.4 .or. n.eq.5 .or. (n.eq.7 .and. igust.gt.0)) then
              write (lun,'(1x,a,i8)') 'number of non-mesonets in the reject list      :', nmp(-5000,n)
              write (lun,'(1x,a,i8)') 'number of mesonet winds in none of the uselists and also not in the reject list:', nmp(-6000,n)
              write (lun,'(1x,a,i8)') 'number of mesonets winds in the rejectlist and in at least one of the uselists :', nmp(-6100,n)
              write (lun,'(1x,a,9x,i8)') 'number of mesonets winds in the rejectlist but in none of the uselists:', nmp(-6200,n)
           endif
           write (lun,'(1x,a,i8)') 'sum of all obs                                 :', sum(nmp(nmplower:nmpupper,n))
        enddo


!123     format(a8,2x,i3,3(2x,f6.2),3(2x,f9.2)) 
!123     format(a8,2x,i3,2(2x,f6.2),4(2x,E11.4)) 
!123     format(a8,4x,i3,2x,2(2x,f6.2),3(2x,E11.4),2x,f5.1) 
!123     format(a8,4x,i3,4x,f6.2,4x,f6.2,2x,3(2x,E11.4),2x,f5.1) 
!123     format(a8,4x,i3,4x,f6.2,4x,f6.2,4x,f6.2,2x,3(2x,E11.4),2x,f5.1) 
!123     format(a8,4x,i3,4x,f6.2,4x,f6.2,4x,f6.2,2x,3(2x,E11.4),2x,f8.1) 
!124     format(a8,4x,i3,4x,f6.2,4x,f6.2,4x,f6.2,2x,3(2x,E11.4),2x,f8.1,7x,a2) 
 123     format(a8,4x,i3,4x,f6.2,4x,f6.2,4x,f6.2,2x,3(2x,E11.4),2x,f8.1,1x,15a) 
 124     format(a8,4x,i3,4x,f6.2,4x,f6.2,4x,f6.2,2x,3(2x,E11.4),2x,f8.1,1x,a18) 
!124     format(a8,4x,i3,4x,f6.2,4x,f6.2,4x,f6.2,2x,3(2x,E11.4),2x,f8.1,1x,15a,' ',a2) 
!124     format(a8,4x,i3,4x,f6.2,4x,f6.2,4x,f6.2,2x,3(2x,E11.4),2x,f8.1,1x,15a,1x,a2) 
!124     format(a8,4x,i3,4x,f6.2,4x,f6.2,4x,f6.2,2x,3(2x,E11.4),2x,f8.1,1x,17a,a2) 
        close(lun_t)
        close(lun_q)
        close(lun_ps)
        close(lun_u)
        close(lun_v)
        close(lun_spd)
        if (igust > 0) close(lun_gust)
        if (ivis  > 0) close(lun_vis)
        if (ipblh > 0) close(lun_pblh)
        if (idist > 0) close(lun_dist)
        close(7)

        if (iter==1) call destroy_auxlists(igust,ivis,ipblh,idist)

        return
        end
!***************************************************************
         subroutine readin_rjlist(filename,rjlist,nmax,fexist,nrjs)

         implicit none

         integer(4),intent(in)::nmax
         integer(4),intent(out)::nrjs
         character(80),intent(in)::filename
         character(90),intent(out)::rjlist(nmax)
         logical,intent(out)::fexist

         integer(4),parameter::meso_unit=61
         integer(4) ncount,m,n
         character(1),parameter::cblank=' '
         character(90) cstring
         character(80) filename2
         character(1) cvar
         integer(4) nrjs0

         do n=1,nmax
           do m=1,90
              rjlist(n)(m:m)=cblank
           enddo
         enddo

         nrjs=0

         fexist=.false.

         inquire(file=trim(filename),exist=fexist)
         if(fexist) then
             open (meso_unit,file=trim(filename),form='formatted')
             ncount=0
             do m=1,3
              read(meso_unit,*,end=151) cstring
             enddo
150          continue
             ncount=ncount+1
             read(meso_unit,*,end=151) rjlist(ncount)
             goto 150
151          continue
             nrjs=ncount-1 ; nrjs=max(nrjs,0)
         endif
         close(meso_unit)
         print*,'in readin_rjlist: filename,nrjs=',trim(filename),nrjs

         if(trim(filename)=='t_rejectlist' .or. & 
            trim(filename)=='q_rejectlist') then

!           cvar=trim(filename)(1:1)
            cvar=filename(1:1)

            do n=1,2
               if (n==1) filename2=cvar//'_day_rejectlist'
               if (n==2) filename2=cvar//'_night_rejectlist'

               inquire(file=trim(filename2),exist=fexist)
               if (fexist) then
                  open (meso_unit,file=trim(filename2),form='formatted')
                  ncount=nrjs ; nrjs0=nrjs
                  do m=1,3
                     read(meso_unit,*,end=251) cstring
                  enddo
250               continue
                  ncount=ncount+1
                  read(meso_unit,*,end=251) rjlist(ncount)
                  goto 250
251               continue
                  nrjs=ncount-1 ; nrjs=max(nrjs,nrjs0)
                  close(meso_unit)
                  print*,'in readin_rjlist: add diurnal lists:n,filename2,nrjs=',n,trim(filename2),nrjs
               endif
            enddo
         endif

         return
         end
!***************************************************************
        subroutine read_rejfileorig(kx,c_station_id,rjlist,nmax,nrjs,clistorig)
        implicit none

         character(8),intent(in)::c_station_id
         character(90),intent(in)::rjlist(nmax)
         character(15),intent(out)::clistorig
         integer(4),intent(in)::kx,nmax,nrjs

         character(8) ch8
         integer(4) m,nlen


         do m=1,nrjs
           ch8(1:8)=rjlist(m)(1:8)
            nlen=len_trim(ch8)
            if ((trim(c_station_id) == trim(ch8)) .or. &
                ((kx==188.or.kx==288.or.kx==195.or.kx==295) .and. c_station_id(1:nlen)==ch8(1:nlen))) then !handle wfo's mesonets which never end with
               clistorig(1:15)=rjlist(m)(70:84)                          !an "a" or "x" in the eight position following blanks
               exit
            endif
         enddo

        return
        end
!***************************************************************
!note:

!***************************************************************
         subroutine init_auxlists(clun3,igust,ivis,ipblh,idist)

         implicit none

         character(3),intent(in)::clun3
         integer(4),intent(in):: igust,ivis,ipblh,idist

         integer(4),parameter::lun_t=100
         integer(4),parameter::lun_q=110
         integer(4),parameter::lun_ps=120
         integer(4),parameter::lun_u=130
         integer(4),parameter::lun_v=140
         integer(4),parameter::lun_spd=150
         integer(4),parameter::lun_gust=160
         integer(4),parameter::lun_vis=170
         integer(4),parameter::lun_pblh=180
         integer(4),parameter::lun_dist=190

         character(96) cheader

         open(lun_t,file='t_obs.listing_iter_'//trim(clun3)//'_aux',form='formatted') !output file

         open(lun_q,file='q_obs.listing_iter_'//trim(clun3)//'_aux',form='formatted') !output file

         open(lun_ps,file='ps_obs.listing_iter_'//trim(clun3)//'_aux',form='formatted') !output file

         open(lun_u,file='u_obs.listing_iter_'//trim(clun3)//'_aux',form='formatted') !output file

         open(lun_v,file='v_obs.listing_iter_'//trim(clun3)//'_aux',form='formatted') !output file

         open(lun_spd,file='spd_obs.listing_iter_'//trim(clun3)//'_aux',form='formatted') !output file

         if (igust > 0) open(lun_gust,file='gust_obs.listing_iter_'//trim(clun3)//'_aux',form='formatted') !output file
         if (ivis > 0) open(lun_vis,file='vis_obs.listing_iter_'//trim(clun3)//'_aux',form='formatted') !output file
         if (ipblh > 0) open(lun_pblh,file='pblh_obs.listing_iter_'//trim(clun3)//'_aux',form='formatted') !output file
         if (idist > 0) open(lun_dist,file='dist_obs.listing_iter_'//trim(clun3)//'_aux',form='formatted') !output file


         cheader='stnname   obtype provider  subprovider         shgt0       hgt0        hgt       slm'

         write(lun_t,'(a)')   'RTMA TEMPERATURE OBS' 
         write(lun_t,'(a)')   'shgt0 ==> station height'
         write(lun_t,'(a)')   'hgt0  ==> observation elevation'
         write(lun_t,'(a)')   'hgt   ==> model terrain at ob (x,y) location'
         write(lun_t,'(a)')   'slm   ==> dominant surface type. this is the surface type of the grid point'
         write(lun_t,'(a)')   '          of the enclosing box nearest to the observation.'
         write(lun_t,'(a)')   '          slm is 0 for water and 1 for land. note: must subtract 3 if slm>=3'
         write(lun_t,'(a)')   '          values of slm>=3 are used to indicate that at least two of the'
         write(lun_t,'(a)')   '          grid points of the enclosing box are of different surface types.'
         write(lun_t,'(a)')   '                              '
         write(lun_t,'(a)')   trim(cheader)

         write(lun_q,'(a)')   'RTMA SPECIFIC HUMIDITY OBS' 
         write(lun_q,'(a)')   'shgt0 ==> station height'
         write(lun_q,'(a)')   'hgt0  ==> observation elevation'
         write(lun_q,'(a)')   'hgt   ==> model terrain at ob (x,y) location'
         write(lun_q,'(a)')   'slm   ==> dominant surface type. this is the surface type of the grid point'
         write(lun_q,'(a)')   '          of the enclosing box nearest to the observation.'
         write(lun_q,'(a)')   '          slm is 0 for water and 1 for land. note: must subtract 3 if slm>=3'
         write(lun_q,'(a)')   '          values of slm>=3 are used to indicate that at least two of the'
         write(lun_q,'(a)')   '          grid points of the enclosing box are of different surface types.'
         write(lun_q,'(a)')   '                              '
         write(lun_q,'(a)')   trim(cheader)

         write(lun_ps,'(a)')   'RTMA SURFACE PRESSURE OBS' 
         write(lun_ps,'(a)')   'shgt0 ==> station height'
         write(lun_ps,'(a)')   'hgt0  ==> observation elevation'
         write(lun_ps,'(a)')   'hgt   ==> model terrain at ob (x,y) location'
         write(lun_ps,'(a)')   'slm   ==> dominant surface type. this is the surface type of the grid point'
         write(lun_ps,'(a)')   '          of the enclosing box nearest to the observation.'
         write(lun_ps,'(a)')   '          slm is 0 for water and 1 for land. note: must subtract 3 if slm>=3'
         write(lun_ps,'(a)')   '          values of slm>=3 are used to indicate that at least two of the'
         write(lun_ps,'(a)')   '          grid points of the enclosing box are of different surface types.'
         write(lun_ps,'(a)')   '                              '
         write(lun_ps,'(a)')   trim(cheader)

         write(lun_u,'(a)')   'RTMA U-WIND OBS' 
         write(lun_u,'(a)')   'shgt0 ==> station height. but ob pressure in hPa for SATWND'
         write(lun_u,'(a)')   'hgt0  ==> observation elevation. for SATWND it is computed from the US  stand. atmos.'
         write(lun_u,'(a)')   'hgt   ==> model terrain at ob (x,y) location'
         write(lun_u,'(a)')   'slm   ==> dominant surface type. this is the surface type of the grid point'
         write(lun_u,'(a)')   '          of the enclosing box nearest to the observation.'
         write(lun_u,'(a)')   '          slm is 0 for water and 1 for land. note: must subtract 3 if slm>=3'
         write(lun_u,'(a)')   '          values of slm>=3 are used to indicate that at least two of the'
         write(lun_u,'(a)')   '          grid points of the enclosing box are of different surface types.'
         write(lun_u,'(a)')   '                              '
         write(lun_u,'(a)')   trim(cheader)

         write(lun_v,'(a)')   'RTMA V-WIND OBS' 
         write(lun_v,'(a)')   'shgt0 ==> station height. but ob pressure in hPa for SATWND'
         write(lun_v,'(a)')   'hgt0  ==> observation elevation. for SATWND it is computed from the US  stand. atmos.'
         write(lun_v,'(a)')   'hgt   ==> model terrain at ob (x,y) location'
         write(lun_v,'(a)')   'slm   ==> dominant surface type. this is the surface type of the grid point'
         write(lun_v,'(a)')   '          of the enclosing box nearest to the observation.'
         write(lun_v,'(a)')   '          slm is 0 for water and 1 for land. note: must subtract 3 if slm>=3'
         write(lun_v,'(a)')   '          values of slm>=3 are used to indicate that at least two of the'
         write(lun_v,'(a)')   '          grid points of the enclosing box are of different surface types.'
         write(lun_v,'(a)')   '                              '
         write(lun_v,'(a)')   trim(cheader)

         write(lun_spd,'(a)')   'RTMA WIND SPEED OBS' 
         write(lun_spd,'(a)')   'shgt0 ==> station height'
         write(lun_spd,'(a)')   'hgt0  ==> observation elevation'
         write(lun_spd,'(a)')   'hgt   ==> model terrain at ob (x,y) location'
         write(lun_spd,'(a)')   'slm   ==> dominant surface type. this is the surface type of the grid point'
         write(lun_spd,'(a)')   '          of the enclosing box nearest to the observation.'
         write(lun_spd,'(a)')   '          slm is 0 for water and 1 for land. note: must subtract 3 if slm>=3'
         write(lun_spd,'(a)')   '          values of slm>=3 are used to indicate that at least two of the'
         write(lun_spd,'(a)')   '          grid points of the enclosing box are of different surface types.'
         write(lun_spd,'(a)')   '                              '
         write(lun_spd,'(a)')   trim(cheader)

         if (igust > 0) then
             write(lun_gust,'(a)')   'RTMA GUST WIND OBS'
             write(lun_gust,'(a)')   'shgt0 ==> station height'
             write(lun_gust,'(a)')   'hgt0  ==> observation elevation'
             write(lun_gust,'(a)')   'hgt   ==> model terrain at ob (x,y) location'
             write(lun_gust,'(a)')   'slm   ==> dominant surface type. this is the surface type of the grid point'
             write(lun_gust,'(a)')   '          of the enclosing box nearest to the observation.'
             write(lun_gust,'(a)')   '          slm is 0 for water and 1 for land. note: must subtract 3 if slm>=3'
             write(lun_gust,'(a)')   '          values of slm>=3 are used to indicate that at least two of the'
             write(lun_gust,'(a)')   '          grid points of the enclosing box are of different surface types.'
             write(lun_gust,'(a)')   '                              '
             write(lun_gust,'(a)')   trim(cheader)
         endif

         if (ivis > 0) then
             write(lun_vis,'(a)')   'RTMA VISIBILITY OBS'
             write(lun_vis,'(a)')   'shgt0 ==> station height'
             write(lun_vis,'(a)')   'hgt0  ==> observation elevation'
             write(lun_vis,'(a)')   'hgt   ==> model terrain at ob (x,y) location'
             write(lun_vis,'(a)')   'slm   ==> dominant surface type. this is the surface type of the grid point'
             write(lun_vis,'(a)')   '          of the enclosing box nearest to the observation.'
             write(lun_vis,'(a)')   '          slm is 0 for water and 1 for land. note: must subtract 3 if slm>=3'
             write(lun_vis,'(a)')   '          values of slm>=3 are used to indicate that at least two of the'
             write(lun_vis,'(a)')   '          grid points of the enclosing box are of different surface types.'
             write(lun_vis,'(a)')   '                              '
             write(lun_vis,'(a)')   trim(cheader)
         endif

         if (ipblh > 0) then
             write(lun_pblh,'(a)')   'RTMA PBLH OBS'
             write(lun_pblh,'(a)')   'shgt0 ==> station height'
             write(lun_pblh,'(a)')   'hgt0  ==> observation elevation'
             write(lun_pblh,'(a)')   'hgt   ==> model terrain at ob (x,y) location'
             write(lun_pblh,'(a)')   'slm   ==> dominant surface type. this is the surface type of the grid point'
             write(lun_pblh,'(a)')   '          of the enclosing box nearest to the observation.'
             write(lun_pblh,'(a)')   '          slm is 0 for water and 1 for land. note: must subtract 3 if slm>=3'
             write(lun_pblh,'(a)')   '          values of slm>=3 are used to indicate that at least two of the'
             write(lun_pblh,'(a)')   '          grid points of the enclosing box are of different surface types.'
             write(lun_pblh,'(a)')   '                              '
             write(lun_pblh,'(a)')   trim(cheader)
         endif

         if (idist > 0) then
             write(lun_dist,'(a)')   'RTMA CEELING HEIGHT OBS'
             write(lun_dist,'(a)')   'shgt0 ==> station height'
             write(lun_dist,'(a)')   'hgt0  ==> observation elevation'
             write(lun_dist,'(a)')   'hgt   ==> model terrain at ob (x,y) location'
             write(lun_dist,'(a)')   'slm   ==> dominant surface type. this is the surface type of the grid point'
             write(lun_dist,'(a)')   '          of the enclosing box nearest to the observation.'
             write(lun_dist,'(a)')   '          slm is 0 for water and 1 for land. note: must subtract 3 if slm>=3'
             write(lun_dist,'(a)')   '          values of slm>=3 are used to indicate that at least two of the'
             write(lun_dist,'(a)')   '          grid points of the enclosing box are of different surface types.'
             write(lun_dist,'(a)')   '                              '
             write(lun_dist,'(a)')   trim(cheader)
         endif

         return
         end

!***************************************************************
         subroutine add2auxlist(otype,cstation,cprovider,csubprovider, & 
                                itype,shgt0,hgt0,hgt,slm)

         implicit none

         character(3),intent(in):: otype
         character(8),intent(in):: cstation
         character(8),intent(inout):: cprovider,csubprovider


         integer(4),parameter::lun_t=100
         integer(4),parameter::lun_q=110
         integer(4),parameter::lun_ps=120
         integer(4),parameter::lun_u=130
         integer(4),parameter::lun_v=140
         integer(4),parameter::lun_spd=150
         integer(4),parameter::lun_gust=160
         integer(4),parameter::lun_vis=170
         integer(4),parameter::lun_pblh=180
         integer(4),parameter::lun_dist=190


         integer(4),intent(in):: itype
         real(4),intent(in):: shgt0,hgt0,hgt,slm
         
         integer(4) lun,lun2

         if (otype(1:3)=='  t')  lun=lun_t
         if (otype(1:3)=='  q')  lun=lun_q
         if (otype(2:3)=='ps')   lun=lun_ps
         if (otype(2:3)=='uv')   lun=lun_u
         if (otype(1:3)=='spd')  lun=lun_spd
         if (otype(1:3)=='gst')  lun=lun_gust
         if (otype(1:3)=='vis')  lun=lun_vis
         if (otype(1:3)=='pbl')  lun=lun_pblh
         !  add dist case here

         lun2=lun_v

         if (cprovider(1:4)=='B7Hv' .or. cprovider(5:8)=='vH7B') then !this provider name comes with weird characters
             cprovider(1:4)='B7Hv'
             cprovider(5:8)='   '
         endif

         if (csubprovider(1:4)=='B7Hv' .or. csubprovider(5:8)=='vH7B') then
             csubprovider(1:4)='B7Hv'
             csubprovider(5:8)='   '
         endif
         

!        write(lun,121) trim(adjustl(cstation)),itype,trim(adjustl(cprovider)), & 
!                       trim(adjustl(csubprovider)),shgt0,hgt0,hgt,slm

         write(lun,121) cstation,itype,cprovider,csubprovider,shgt0,hgt0,hgt,slm

         if (otype(2:3)=='uv') then 
!         write(lun2,121) trim(adjustl(cstation)),itype,trim(adjustl(cprovider)), & 
!                         trim(adjustl(csubprovider)),shgt0,hgt0,hgt,slm
         write(lun2,121) cstation,itype,cprovider,csubprovider,shgt0,hgt0,hgt,slm
         endif

121      format(a8,4x,i3,4x,a8,4x,a8,4x,f10.3,2x,e10.3,2x,f10.3,2x,f7.3)
!!!121      format(a8,4x,i3,4x,a8,4x,a8,4x,3(f10.3,2x),f7.3)

         return
         end
!***************************************************************
         subroutine destroy_auxlists(igust,ivis,ipblh,idist)

         implicit none

         integer(4),intent(in):: igust,ivis,ipblh,idist

         integer(4),parameter::lun_t=100
         integer(4),parameter::lun_q=110
         integer(4),parameter::lun_ps=120
         integer(4),parameter::lun_u=130
         integer(4),parameter::lun_v=140
         integer(4),parameter::lun_spd=150
         integer(4),parameter::lun_gust=160
         integer(4),parameter::lun_vis=170
         integer(4),parameter::lun_pblh=180
         integer(4),parameter::lun_dist=190

         close(lun_t)
         close(lun_q)
         close(lun_ps)
         close(lun_u)
         close(lun_v)
         close(lun_spd)
         if (igust > 0) close(lun_gust)
         if (ivis  > 0) close(lun_vis)
         if (ipblh > 0) close(lun_pblh)
         if (idist > 0) close(lun_dist)

         return
         end
!***************************************************************
!        if (otype(3:3)=='t') then
!***************************************************************
!       rstation_id     = data(id,i)
!       cdiagbuf(ii)    = station_id         ! station id
! 
!       rdiagbuf(1,ii)  = ictype(ikx)        ! observation type
!       rdiagbuf(2,ii)  = icsubtype(ikx)     ! observation subtype
!   
!       rdiagbuf(3,ii)  = data(ilate,i)      ! observation latitude (degrees)
!       rdiagbuf(4,ii)  = data(ilone,i)      ! observation longitude (degrees)
!       rdiagbuf(5,ii)  = data(istnelv,i)    ! station elevation (meters)
!       rdiagbuf(6,ii)  = prest              ! observation pressure (hPa)
!       rdiagbuf(7,ii)  = data(iobshgt,i)    ! observation height (meters)
!       rdiagbuf(8,ii)  = dtime              ! obs time (hours relative to analysis time)

!       rdiagbuf(9,ii)  = data(iqc,i)        ! input prepbufr qc or event mark
!       rdiagbuf(10,ii) = data(iqt,i)        ! setup qc or event mark (currently qtflg only)
!       rdiagbuf(11,ii) = data(iuse,i)       ! read_prepbufr data usage flag
!       if(muse(i)) then
!          rdiagbuf(12,ii) = one             ! analysis usage flag (1=use, -1=not used)
!       else
!          rdiagbuf(12,ii) = -one
!       endif
!
!       err_input = data(ier2,i)
!       err_adjst = data(ier,i)
!       if (ratio_errors*error>tiny_r_kind) then
!          err_final = one/(ratio_errors*error)
!       else
!          err_final = huge_single
!       endif
!
!       errinv_input = huge_single
!       errinv_adjst = huge_single
!       errinv_final = huge_single
!       if (err_input>tiny_r_kind) errinv_input=one/err_input
!       if (err_adjst>tiny_r_kind) errinv_adjst=one/err_adjst
!       if (err_final>tiny_r_kind) errinv_final=one/err_final
!
!       rdiagbuf(13,ii) = rwgt               ! nonlinear qc relative weight
!       rdiagbuf(14,ii) = errinv_input       ! prepbufr inverse obs error (K**-1)
!       rdiagbuf(15,ii) = errinv_adjst       ! read_prepbufr inverse obs error (K**-1)
!       rdiagbuf(16,ii) = errinv_final       ! final inverse observation error (K**-1)
!
!       rdiagbuf(17,ii) = data(itob,i)       ! temperature observation (K)
!       rdiagbuf(18,ii) = ddiff              ! obs-ges used in analysis (K)
!       rdiagbuf(19,ii) = tob-tges           ! obs-ges w/o bias correction (K) (future slot)
!***************************************************************

 
!***************************************************************
!       if (otype(3:3)=='q') then
!***************************************************************
!       rstation_id     = data(id,i)
!       cdiagbuf(ii)    = station_id         ! station id
!
!       rdiagbuf(1,ii)  = ictype(ikx)        ! observation type
!       rdiagbuf(2,ii)  = icsubtype(ikx)     ! observation subtype
!   
!       rdiagbuf(3,ii)  = data(ilate,i)      ! observation latitude (degrees)
!       rdiagbuf(4,ii)  = data(ilone,i)      ! observation longitude (degrees)
!       rdiagbuf(5,ii)  = data(istnelv,i)    ! station elevation (meters)
!       rdiagbuf(6,ii)  = presq              ! observation pressure (hPa)
!       rdiagbuf(7,ii)  = data(iobshgt,i)    ! observation height (meters)
!       rdiagbuf(8,ii)  = dtime              ! obs time (hours relative to analysis time)
! 
!       rdiagbuf(9,ii)  = data(iqc,i)        ! input prepbufr qc or event mark
!       rdiagbuf(10,ii) = rmiss_single       ! setup qc or event mark 
!       rdiagbuf(11,ii) = data(iuse,i)       ! read_prepbufr data usage flag
!       if(muse(i)) then
!          rdiagbuf(12,ii) = one             ! analysis usage flag (1=use, -1=not used)
!       else
!          rdiagbuf(12,ii) = -one                    
!       endif
! 
!       err_input = data(ier2,i)*qsges            ! convert rh to q
!       err_adjst = data(ier,i)*qsges             ! convert rh to q
!       if (ratio_errors*error>tiny_r_kind) then
!          err_final = one/(ratio_errors*error)
!       else
!          err_final = huge_single
!       endif
!
!       errinv_input = huge_single
!       errinv_adjst = huge_single
!       errinv_final = huge_single
!       if (err_input>tiny_r_kind) errinv_input = one/err_input
!       if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst
!       if (err_final>tiny_r_kind) errinv_final = one/err_final
!
!       rdiagbuf(13,ii) = rwgt               ! nonlinear qc relative weight
!       rdiagbuf(14,ii) = errinv_input       ! prepbufr inverse observation error
!       rdiagbuf(15,ii) = errinv_adjst       ! read_prepbufr inverse obs error
!       rdiagbuf(16,ii) = errinv_final       ! final inverse observation error
!
!       rdiagbuf(17,ii) = data(iqob,i)       ! observation
!       rdiagbuf(18,ii) = ddiff              ! obs-ges used in analysis
!       rdiagbuf(19,ii) = qob-qges           ! obs-ges w/o bias correction (future slot)
!
!       rdiagbuf(20,ii) = qsges              ! guess saturation specific humidity
!       
!***************************************************************


!***************************************************************
!       if (otype(2:3)=='ps') then
!***************************************************************
!       rstation_id     = data(id,i)
!       cdiagbuf(ii)    = station_id         ! station id
!
!       rdiagbuf(1,ii)  = ictype(ikx)        ! observation type
!       rdiagbuf(2,ii)  = icsubtype(ikx)     ! observation subtype
!       rdiagbuf(3,ii)  = data(ilate,i)      ! observation latitude (degrees)
!       rdiagbuf(4,ii)  = data(ilone,i)      ! observation longitude (degrees)
!       rdiagbuf(5,ii)  = data(istnelv,i)    ! station elevation (meters)
!       rdiagbuf(6,ii)  = data(ipres,i)*r10  ! observation pressure (hPa)
!       rdiagbuf(7,ii)  = dhgt               ! observation height (meters)
!       rdiagbuf(8,ii)  = dtime              ! obs time (hours relative to analysis time)
!
!       rdiagbuf(9,ii)  = data(iqc,i)        ! input prepbufr qc or event mark
!       rdiagbuf(10,ii) = rmiss_single       ! setup qc or event mark
!       rdiagbuf(11,ii) = data(iuse,i)       ! read_prepbufr data usage flag
!       if(muse(i)) then
!          rdiagbuf(12,ii) = one             ! analysis usage flag (1=use, -1=not used)
!       else
!          rdiagbuf(12,ii) = -one                    
!       endif
!
!       pob      = pob*r10
!       pges     = pges*r10
!       pgesorig = pgesorig*r10
!
!       err_input = data(ier2,i)*r10   ! r10 converts cb to mb
!       err_adjst = data(ier,i)*r10
!       if (ratio_errors*error/r10>tiny_r_kind) then
!          err_final = r10/(ratio_errors*error)
!       else
!          err_final = huge_single
!       endif
!
!       errinv_input = huge_single
!       errinv_adjst = huge_single
!       errinv_final = huge_single
!       if (err_input>tiny_r_kind) errinv_input = one/err_input
!       if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst
!       if (err_final>tiny_r_kind) errinv_final = one/err_final
!
!       rdiagbuf(13,ii) = rwgt               ! nonlinear qc relative weight
!       rdiagbuf(14,ii) = errinv_input       ! prepbufr inverse obs error (hPa**-1)
!       rdiagbuf(15,ii) = errinv_adjst       ! read_prepbufr inverse obs error (hPa**-1)
!       rdiagbuf(16,ii) = errinv_final       ! final inverse observation error (hPa**-1)
!
!       rdiagbuf(17,ii) = pob                ! surface pressure observation (hPa)
!       rdiagbuf(18,ii) = pob-pges           ! obs-ges used in analysis (coverted to hPa)
!       rdiagbuf(19,ii) = pob-pgesorig       ! obs-ges w/o adjustment to guess surface pressure (hPa)
!***************************************************************


!***************************************************************
!       if (otype(2:3)=='uv') then
!***************************************************************
!       rstation_id     = data(id,i)
!       cdiagbuf(ii)    = station_id         ! station id
!
!       rdiagbuf(1,ii)  = ictype(ikx)        ! observation type
!       rdiagbuf(2,ii)  = icsubtype(ikx)     ! observation subtype
!   
!       rdiagbuf(3,ii)  = data(ilate,i)      ! observation latitude (degrees)
!       rdiagbuf(4,ii)  = data(ilone,i)      ! observation longitude (degrees)
!       rdiagbuf(5,ii)  = data(ielev,i)      ! station elevation (meters)
!       rdiagbuf(6,ii)  = presw              ! observation pressure (hPa)
!       rdiagbuf(7,ii)  = data(ihgt,i)       ! observation height (meters)
!       rdiagbuf(8,ii)  = dtime              ! obs time (hours relative to analysis time)
!
!       rdiagbuf(9,ii)  = data(iqc,i)        ! input prepbufr qc or event mark
!       rdiagbuf(10,ii) = rmiss_single       ! setup qc or event mark
!       rdiagbuf(11,ii) = data(iuse,i)       ! read_prepbufr data usage flag
!       if(muse(i)) then
!          rdiagbuf(12,ii) = one             ! analysis usage flag (1=use, -1=not used)
!       else
!          rdiagbuf(12,ii) = -one
!       endif
!
!       err_input = data(ier2,i)
!       err_adjst = data(ier,i)
!       if (ratio_errors*error>tiny_r_kind) then
!          err_final = one/(ratio_errors*error)
!       else
!          err_final = huge_single
!       endif
!
!       errinv_input = huge_single
!       errinv_adjst = huge_single
!       errinv_final = huge_single
!       if (err_input>tiny_r_kind) errinv_input = one/err_input
!       if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst
!       if (err_final>tiny_r_kind) errinv_final = one/err_final
!
!       rdiagbuf(13,ii) = rwgt               ! nonlinear qc relative weight
!       rdiagbuf(14,ii) = errinv_input       ! prepbufr inverse obs error (m/s)**-1
!       rdiagbuf(15,ii) = errinv_adjst       ! read_prepbufr inverse obs error (m/s)**-1
!       rdiagbuf(16,ii) = errinv_final       ! final inverse observation error (m/s)**-1
!
!       rdiagbuf(17,ii) = data(iuob,i)       ! u wind component observation (m/s)
!       rdiagbuf(18,ii) = dudiff             ! u obs-ges used in analysis (m/s)
!       rdiagbuf(19,ii) = uob-ugesin         ! u obs-ges w/o bias correction (m/s) (future slot)
!
!       rdiagbuf(20,ii) = data(ivob,i)       ! v wind component observation (m/s)
!       rdiagbuf(21,ii) = dvdiff             ! v obs-ges used in analysis (m/s)
!       rdiagbuf(22,ii) = vob-vgesin         ! v obs-ges w/o bias correction (m/s) (future slot)
!
!       rdiagbuf(23,ii) = factw              ! 10m wind reduction factor
!***************************************************************


!***************************************************************
!       if (otype(1:3)=='spd') then
!***************************************************************
!       rstation_id     = data(id,i)
!       cdiagbuf(ii)    = station_id         ! station id
!
!       rdiagbuf(1,ii)  = ictype(ikx)        ! observation type
!       rdiagbuf(2,ii)  = icsubtype(ikx)     ! observation subtype
!    
!       rdiagbuf(3,ii)  = data(ilate,i)      ! observation latitude (degrees)
!       rdiagbuf(4,ii)  = data(ilone,i)      ! observation longitude (degrees)
!       rdiagbuf(5,ii)  = data(istnelv,i)    ! station elevation (meters)
!       rdiagbuf(6,ii)  = presw              ! observation pressure (hPa)
!       rdiagbuf(7,ii)  = data(ihgt,i)       ! observation height (meters)
!       rdiagbuf(8,ii)  = dtime              ! obs time (hours relative to analysis time)
!
!       rdiagbuf(9,ii)  = data(iqc,i)        ! input prepbufr qc or event mark
!       rdiagbuf(10,ii) = rmiss_single       ! setup qc or event mark
!       rdiagbuf(11,ii) = data(iuse,i)       ! read_prepbufr data usage flag
!       if(muse(i)) then
!          rdiagbuf(12,ii) = one             ! analysis usage flag (1=use, -1=not used)
!       else
!          rdiagbuf(12,ii) = -one
!       endif
!
!       spdob0    = sqrt(data(iuob,i)*data(iuob,i)+data(ivob,i)*data(ivob,i))
!       err_input = data(ier2,i)
!       err_adjst = data(ier,i)
!       if (ratio_errors*error>tiny_r_kind) then
!          err_final = one/(ratio_errors*error)
!       else
!          err_final = huge_single
!       endif
!
!       errinv_input = huge_single
!       errinv_adjst = huge_single
!       errinv_final = huge_single
!       if (err_input>tiny_r_kind) errinv_input = one/err_input
!       if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst
!       if (err_final>tiny_r_kind) errinv_final = one/err_final
!
!       rdiagbuf(13,ii) = rwgt               ! nonlinear qc relative weight
!       rdiagbuf(14,ii) = errinv_input       ! prepbufr inverse obs error (m/s)**-1
!       rdiagbuf(15,ii) = errinv_adjst       ! read_prepbufr inverse obs error (m/s)**-1
!       rdiagbuf(16,ii) = errinv_final       ! final inverse observation error (m/s)**-1
!
!       rdiagbuf(17,ii) = spdob              ! wind speed observation (m/s)
!       rdiagbuf(18,ii) = ddiff              ! obs-ges used in analysis (m/s)
!       rdiagbuf(19,ii) = spdob0-spdges      ! obs-ges w/o bias correction (m/s) (future slot)
!
!       rdiagbuf(20,ii) = factw              ! 10m wind reduction factor
!***************************************************************
