subroutine read_wcpbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
     prsl_full,nobs,nrec_start)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:  read_wcpbufr:  read obs from wcpbufr file
!   prgmmr: parrish          org: np22                date: 1990-10-07
!
! abstract:  This routine reads retrieved hydrometeor (water content) data in the wcpbufr file.  
!            Specific observation types read by this routine include: 
!            solid-water content path and liquid-water content path
!            derived from Hurricane GPROF (see Wu et al. 2016, Brown et al. 2016)
!            (they are called integrated solid-water content and integrated liquid-water content 
!             in Wu et al. 2016) 
!
!            When running the gsi in regional mode, the code only
!            retains those observations that fall within the regional
!            domain
!
! program history log:
!   2017-12-18  T.-C. Wu - adapted from read_prepbufr

!   input argument list:
!     infile   - unit from which to read BUFR data
!     obstype  - observation type to process
!     lunout   - unit to which to write data for further processing
!     prsl_full- 3d pressure on full domain grid
!     nrec_start - number of subsets without useful information
!
!   output argument list:
!     nread    - number of type "obstype" observations read
!     nodata   - number of individual "obstype" observations read
!     ndata    - number of type "obstype" observations retained for further processing
!     twindin  - input group time window (hours)
!     sis      - satellite/instrument/sensor indicator
!     nobs     - array of observations on each subdomain for each processor
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_single,r_kind,r_double,i_kind
  use constants, only: zero,one_tenth,one,deg2rad,half,&
      rad2deg,tiny_r_kind,huge_r_kind,huge_i_kind,&
      r60inv,r2000
  use gridmod, only: diagnostic_reg,regional,nlon,nlat,nsig,&
      tll2xy,txy2ll, rlats,rlons
  use convinfo, only: nconvtype,ctwind, &
      ncmiter,ncgroup,ncnumgrp,icuse,ictype,icsubtype,ioctype, &
      ithin_conv,rmesh_conv,pmesh_conv
  use converr,only: etabl
  use obsmod, only: iadate, offtime_data, oberrflg
  use gsi_4dvar, only: l4dvar,l4densvar,time_4dvar,winlen,thin4d
  use convthin, only: make3grids,map3grids,del3grids,use_all
  use mpimod, only: npe

  implicit none

! Declare passed variables
  character(len=*)                      ,intent(in   ) :: infile,obstype
  character(len=20)                     ,intent(in   ) :: sis
  integer(i_kind)                       ,intent(in   ) :: lunout,nrec_start
  integer(i_kind)                       ,intent(inout) :: nread,ndata,nodata
  integer(i_kind),dimension(npe)        ,intent(inout) :: nobs
  real(r_kind)                          ,intent(in   ) :: twindin
  real(r_kind),dimension(nlat,nlon,nsig),intent(in   ) :: prsl_full

! Declare local parameters
  real(r_kind),parameter:: r6   = 6.0_r_kind
  real(r_kind),parameter:: r90  = 90.0_r_kind
  real(r_kind),parameter:: r360 = 360.0_r_kind
  real(r_kind),parameter:: r1200= 1200.0_r_kind
  real(r_kind),parameter:: convert= 1.0e-3_r_kind ! from g m^-2 to kg m^-2

! Declare local variables
  logical swcpob, lwcpob
  logical outside
  logical luse,ithinp
  logical,allocatable,dimension(:,:):: lmsg           ! set true when convinfo entry id found in a message

  character(40) hdstr,qcstr,oestr,levstr,hdstr2
  character(80) obstr
  character(10) date
  character(8) subset
  character(8) c_station_id
  character(1) sidchr(8)

  integer(i_kind) ireadmg,ireadsb,icntpnt,icntpnt2,icount,iiout
  integer(i_kind) lunin,i,maxobs,nmsgmax,mxtb
  integer(i_kind) kk,klon1,klat1,klonp1,klatp1
  integer(i_kind) nc,nx,ntread,itx,ii,ncsave
  integer(i_kind) ihh,idd,idate,iret,im,iy,k,levs
  integer(i_kind) kx,nreal,nchanl,ilat,ilon,ithin
  integer(i_kind) qm, swcpq, lwcpq
  integer(i_kind) nlevp         ! vertical level for thinning
  integer(i_kind) ntmp,iout
  integer(i_kind) pflag,irec
  integer(i_kind) ntest,nvtest,iosub,ixsub,isubsub,iobsub
  integer(i_kind) kl,k1,k2
  integer(i_kind) itypex
  integer(i_kind) minobs,minan
  integer(i_kind) ntb,ntmatch,ncx
  integer(i_kind) nmsg                ! message index
  integer(i_kind),dimension(5):: idate5
  integer(i_kind),dimension(255):: pqm
  integer(i_kind),dimension(nconvtype)::ntxall
  integer(i_kind),dimension(nconvtype+1)::ntx
  integer(i_kind),allocatable,dimension(:):: isort,iloc,nrep
  integer(i_kind),allocatable,dimension(:,:):: tab
  real(r_kind) time,timex,timeobs,toff,t4dv,zeps
  real(r_kind) rmesh,ediff,usage
  real(r_kind) dx,dy,dx1,dy1,w00,w10,w01,w11
  real(r_kind) dlnpob,ppb
  real(r_kind) swcpoe, swcpmerr, lwcpoe, lwcpmerr
  real(r_kind) dlat,dlon,dlat_earth,dlon_earth
  real(r_kind) dlat_earth_deg,dlon_earth_deg
  real(r_kind) stnelev
  real(r_kind) cdist,disterr,disterrmax,rlon00,rlat00
  real(r_kind) vdisterrmax
  real(r_kind) del, swcperrmin, lwcperrmin
  real(r_kind) crit1,timedif,xmesh,pmesh
  real(r_kind) time_correction
  real(r_kind) perrmin
  real(r_kind),dimension(nsig):: presl
  real(r_kind),dimension(nsig-1):: dpres
  real(r_kind),dimension(255)::plevs
  real(r_kind),allocatable,dimension(:):: presl_thin
  real(r_kind),allocatable,dimension(:,:):: cdata_all,cdata_out

  real(r_double) rstation_id,qcmark_huge
  real(r_double),dimension(8):: hdr
  real(r_double),dimension(4,255):: qcmark,obserr
  real(r_double),dimension(5,255):: obsdat
  real(r_double),dimension(1,255):: levdat
  equivalence(rstation_id,c_station_id)
  equivalence(rstation_id,sidchr)

!  data statements
  data hdstr  /'SID XOB YOB DHR TYP ELV SAID T29'/
  data hdstr2 /'TYP SAID T29 SID'/
  data obstr  /'POB ZOB CWIO CWLO PRSS' /
  data qcstr  /'PQM ZQM CWIQ CWLQ     '/
  data oestr  /'POE NUL CWIE CWLE     '/
  data levstr  /'POB'/

  data lunin / 15 /
  data ithin / -9 /
  data rmesh / -99.999_r_kind /

  
!------------------------------------------------------------------------
! Initialize variables

  vdisterrmax=zero
  pflag=0                  !  dparrish debug compile run flags pflag as not defined ???????????

  swcpob = obstype == 'swcp'
  lwcpob = obstype == 'lwcp'
  if(swcpob) then
     nreal=16
  else if(lwcpob) then
     nreal=16
  else 
     write(6,*) ' illegal obs type in READ_WCPBUFR ',obstype
     call stop2(94)
  end if

  qcmark_huge = huge_i_kind

  perrmin=0.3_r_kind
  swcperrmin=one
  lwcperrmin=one
  
!------------------------------------------------------------------------
  ntread=1
  ntmatch=0
  ntx(ntread)=0
  ntxall=0
  do nc=1,nconvtype
     if(trim(ioctype(nc)) == trim(obstype))then
          ntmatch=ntmatch+1
          ntxall(ntmatch)=nc
     end if
     if(trim(ioctype(nc)) == trim(obstype) .and. abs(icuse(nc)) <= 1)then
           ithin=ithin_conv(nc)
           if(ithin > 0)then
              ntread=ntread+1
              ntx(ntread)=nc
           end if
     end if
  end do
  if(ntmatch == 0)then
     write(6,*) ' no matching obstype found in obsinfo ',obstype
     return
  end if

!! get message and subset counts

  call getcount_bufr(infile,nmsgmax,mxtb)
  allocate(lmsg(nmsgmax,ntread),tab(mxtb,3),nrep(nmsgmax))

  lmsg = .false.
  maxobs=0
  tab=0
  nmsg=0
  nrep=0
  ntb = 0
  irec = 0

! Open, then read date from bufr data
  call closbf(lunin)
  open(lunin,file=trim(infile),form='unformatted')
  call openbf(lunin,'IN',lunin)
  call datelen(10)

  msg_report: do while (ireadmg(lunin,subset,idate) == 0)
     irec = irec + 1
     if(irec < nrec_start) cycle msg_report

!    Time offset
     if(nmsg == 0) call time_4dvar(idate,toff)
     nmsg=nmsg+1
     if (nmsg>nmsgmax) then
        write(6,*)'READ_WCPBUFR: messages exceed maximum ',nmsgmax
        call stop2(50)
     endif
     loop_report: do while (ireadsb(lunin) == 0)
        ntb = ntb+1
        nrep(nmsg)=nrep(nmsg)+1
        if (ntb>mxtb) then
           write(6,*)'READ_WCPBUFR: reports exceed maximum ',mxtb
           call stop2(50)
        endif

!       Extract type information
        call ufbint(lunin,hdr,4,1,iret,hdstr2)
        kx=hdr(1)

! temporary specify iobsub until put in bufr file
        iobsub = 0                                                  

!  Match ob to proper convinfo type
        ncsave=0
        matchloop:do ncx=1,ntmatch
           nc=ntxall(ncx)
           if (kx /= ictype(nc))cycle 

!  Find convtype which match ob type and subtype
           if(icsubtype(nc) == iobsub) then
              ncsave=nc
              exit matchloop
           else
!  Find convtype which match ob type and subtype group (isubtype == ?*)
!       where ? specifies the group and icsubtype = ?0)
              ixsub=icsubtype(nc)/10
              iosub=iobsub/10
              isubsub=icsubtype(nc)-ixsub*10
              if(ixsub == iosub .and. isubsub == 0) then
                 ncsave=nc
!  Find convtype which match ob type and subtype is all remaining 
!       (icsubtype(nc) = 0)
              else if (ncsave == 0 .and. icsubtype(nc) == 0) then
                 ncsave=nc
              end if
           end if
        end do matchloop

!  Save information for next read
        if(ncsave /= 0) then

           call ufbint(lunin,levdat,1,255,levs,levstr)
           maxobs=maxobs+max(1,levs)
           nx=1
           if(ithin_conv(ncsave) > 0)then
              do ii=2,ntread
                 if(ntx(ii) == ncsave)nx=ii
              end do
           end if
           tab(ntb,1)=ncsave
           tab(ntb,2)=nx
           tab(ntb,3)=levs
           lmsg(nmsg,nx) = .true.
        end if

     end do loop_report
  enddo msg_report
  if (nmsg==0) then
     write(6,*)'READ_WCPBUFR: no messages/reports '
     call closbf(lunin)
     close(lunin)
     return
  end if
  write(6,*)'READ_WCPBUFR: messages/reports = ',nmsg,'/',ntb,' ntread = ',ntread
!------------------------------------------------------------------------

! loop over convinfo file entries; operate on matches
  
  allocate(cdata_all(nreal,maxobs),isort(maxobs))
  isort = 0
  cdata_all=zero
  nread=0
  ntest=0
  nvtest=0
  nchanl=0
  ilon=2
  ilat=3
  loop_convinfo: do nx=1, ntread

     use_all = .true.
     ithin=0
     if(nx > 1) then
        nc=ntx(nx)
        ithin=ithin_conv(nc)
        if (ithin > 0 ) then
           rmesh=rmesh_conv(nc)
           pmesh=pmesh_conv(nc)
           use_all = .false.
           if(pmesh > zero) then
              pflag=1
              nlevp=r1200/pmesh
           else
              pflag=0
              nlevp=nsig
           endif
           xmesh=rmesh

           call make3grids(xmesh,nlevp)

           if (.not.use_all) then
              allocate(presl_thin(nlevp))
              if (pflag==1) then
                 do k=1,nlevp
                    presl_thin(k)=(r1200-(k-1)*pmesh)*one_tenth
                 enddo
              endif
           endif
     
           write(6,*)'READ_WCPBUFR: obstype,ictype(nc),rmesh,pflag,nlevp,pmesh=',&
              trim(ioctype(nc)),ictype(nc),rmesh,pflag,nlevp,pmesh
        endif
     endif
       

     call closbf(lunin)
     open(lunin,file=infile,form='unformatted')
     call openbf(lunin,'IN',lunin)
     call datelen(10)

!    Big loop over prepbufr file	

     ntb = 0
     nmsg = 0
     icntpnt=0
     icntpnt2=0
     disterrmax=-9999.0_r_kind
     irec = 0
     loop_msg: do while (ireadmg(lunin,subset,idate)== 0)
        irec = irec + 1
        if(irec < nrec_start) cycle loop_msg

        nmsg = nmsg+1
        if(.not.lmsg(nmsg,nx)) then
           do i=ntb+1,ntb+nrep(nmsg)
              icntpnt2=icntpnt2+tab(i,3)
           end do
           ntb=ntb+nrep(nmsg)
           cycle loop_msg ! no useable reports this mesage, skip ahead report count
        end if 

        loop_readsb: do while(ireadsb(lunin) == 0)
!          use msg lookup table to decide which messages to skip
!          use report id lookup table to only process matching reports
           ntb = ntb+1
           if(icntpnt < icntpnt2)icntpnt=icntpnt2
           icntpnt2=icntpnt2+tab(ntb,3)
           nc=tab(ntb,1)
           if(nc <= 0 .or. tab(ntb,2) /= nx) cycle loop_readsb
                 
!          Extract type, date, and location information
           call ufbint(lunin,hdr,8,1,iret,hdstr)
           kx=hdr(5)

              if(abs(hdr(3))>r90 .or. abs(hdr(2))>r360) cycle loop_readsb
              if(hdr(2)== r360)hdr(2)=hdr(2)-r360
              if(hdr(2) < zero)hdr(2)=hdr(2)+r360
              dlon_earth_deg=hdr(2)
              dlat_earth_deg=hdr(3)
              dlon_earth=hdr(2)*deg2rad
              dlat_earth=hdr(3)*deg2rad

              if(regional)then
                 call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)    ! convert to rotated coordinate
                 if(diagnostic_reg) then
                    call txy2ll(dlon,dlat,rlon00,rlat00)
                    ntest=ntest+1
                    cdist=sin(dlat_earth)*sin(rlat00)+cos(dlat_earth)*cos(rlat00)* &
                         (sin(dlon_earth)*sin(rlon00)+cos(dlon_earth)*cos(rlon00))
                    cdist=max(-one,min(cdist,one))
                    disterr=acos(cdist)*rad2deg
                    disterrmax=max(disterrmax,disterr)
                 end if
                 if(outside) cycle loop_readsb   ! check to see if outside regional domain
              else
                 dlat = dlat_earth
                 dlon = dlon_earth
                 call grdcrd1(dlat,rlats,nlat,1)
                 call grdcrd1(dlon,rlons,nlon,1)
              endif

!------------------------------------------------------------------------

           if(offtime_data) then
 
!             in time correction for observations to account for analysis
!                      time being different from obs file time.
              write(date,'( i10)') idate
              read (date,'(i4,3i2)') iy,im,idd,ihh
              idate5(1)=iy
              idate5(2)=im
              idate5(3)=idd
              idate5(4)=ihh
              idate5(5)=0
              call w3fs21(idate5,minobs)    !  obs ref time in minutes relative to historic date
              idate5(1)=iadate(1)
              idate5(2)=iadate(2)
              idate5(3)=iadate(3)
              idate5(4)=iadate(4)
              idate5(5)=0
              call w3fs21(idate5,minan)    !  analysis ref time in minutes relative to historic date
 
!             Add obs reference time, then subtract analysis time to get obs time relative to analysis
 
              time_correction=float(minobs-minan)*r60inv

           else
              time_correction=zero
           end if

              timeobs=real(real(hdr(4),r_single),r_double)
              t4dv=timeobs + toff
              zeps=1.0e-8_r_kind
              if (t4dv<zero  .and.t4dv>      -zeps) t4dv=zero
              if (t4dv>winlen.and.t4dv<winlen+zeps) t4dv=winlen
              t4dv=t4dv + time_correction
              time=timeobs + time_correction

 
              if (l4dvar.or.l4densvar) then
                 if (t4dv<zero.OR.t4dv>winlen) cycle loop_readsb ! outside time window
              else
                 if((real(abs(time)) > real(ctwind(nc)) .or. real(abs(time)) > real(twindin)))cycle loop_readsb ! outside time window
              endif

              timex=time
     
!          Extract data information on levels
           call ufbint(lunin,obsdat,5,255,levs,obstr)
           call ufbint(lunin,qcmark,4,255,levs,qcstr)
           call ufbint(lunin,obserr,4,255,levs,oestr)

           if(oberrflg)then

!             Set lower limits for observation errors
              swcperrmin=one_tenth
              lwcperrmin=one_tenth
                do k=1,levs
                   itypex=kx
                   ppb=obsdat(1,k)
                   ppb=max(zero,min(ppb,r2000))
                   if(ppb>=etabl(itypex,1,1)) k1=1
                   do kl=1,32
                      if(ppb>=etabl(itypex,kl+1,1).and.ppb<=etabl(itypex,kl,1)) k1=kl
                   end do
                   if(ppb<=etabl(itypex,33,1)) k1=5
                   k2=k1+1
                   ediff = etabl(itypex,k2,1)-etabl(itypex,k1,1)
                   if (abs(ediff) > tiny_r_kind) then
                      del = (ppb-etabl(itypex,k1,1))/ediff
                   else
                      del = huge_r_kind
                   endif
                   del=max(zero,min(del,one))
                   obserr(1,k)=(one-del)*etabl(itypex,k1,5)+del*etabl(itypex,k2,5)
                   obserr(1,k)=max(obserr(1,k),perrmin)
                   obserr(3,k)=max(obserr(3,k),swcperrmin)
                   obserr(4,k)=max(obserr(3,k),lwcperrmin)
                enddo
           endif        ! endif for oberrflg

           nread=nread+1

!          Set station ID
           rstation_id=hdr(1)

!          Loop over levels
           do k=1,levs
              do i=1,4
                 qcmark(i,k) = min(qcmark(i,k),qcmark_huge)
              end do

!              if (kx == id_bias_ps) then
!                 plevs(k)=one_tenth*obsdat(1,k)+conv_bias_ps   ! convert mb to cb
!              else
                 plevs(k)=one_tenth*obsdat(1,k)   ! convert mb to cb
!              endif
              pqm(k)=nint(qcmark(1,k))
           end do

           stnelev=hdr(6)
           ithin=ithin_conv(nc)
           ithinp = ithin > 0 .and. pflag /= 0
           if(levs > 1 .or. ithinp)then
!             Interpolate guess pressure profile to observation location
              klon1= int(dlon);  klat1= int(dlat)
              dx   = dlon-klon1; dy   = dlat-klat1
              dx1  = one-dx;     dy1  = one-dy
              w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy
 
              klat1=min(max(1,klat1),nlat); klon1=min(max(0,klon1),nlon)
              if (klon1==0) klon1=nlon
              klatp1=min(nlat,klat1+1); klonp1=klon1+1
              if (klonp1==nlon+1) klonp1=1
              do kk=1,nsig
                 presl(kk)=w00*prsl_full(klat1 ,klon1 ,kk) +  &
                           w10*prsl_full(klatp1,klon1 ,kk) + &
                           w01*prsl_full(klat1 ,klonp1,kk) + &
                           w11*prsl_full(klatp1,klonp1,kk)
              end do

!             Compute depth of guess pressure layersat observation location
              if (levs > 1) then
                 do kk=1,nsig-1
                    dpres(kk)=presl(kk)-presl(kk+1)
                 end do
              endif
           end if
           LOOP_K_LEVS: do k=1,levs
              icntpnt=icntpnt+1

!             Extract quality marks
              if(swcpob) then
                 swcpq=nint(qcmark(3,k))
                 qm=swcpq
              else if(lwcpob) then
                 lwcpq=nint(qcmark(4,k))
                 qm=lwcpq
             end if

!             Check qc marks to see if obs should be processed or skipped

              if(qm > 15 .or. qm < 0) cycle loop_k_levs

!             Special block for data thinning - if requested
              if (ithin > 0) then
                 ntmp=ndata  ! counting moved to map3gridS
           
!                Set data quality index for thinning
                 if (thin4d) then
                    timedif = zero
                 else
                    timedif=abs(t4dv-toff)
                 endif
                    crit1 = timedif/r6+half

                 if (pflag==0) then
                    do kk=1,nsig
                       presl_thin(kk)=presl(kk)
                    end do
                 endif

                 call map3grids(-1,pflag,presl_thin,nlevp,dlat_earth,dlon_earth,&
                    plevs(k),crit1,ndata,iout,icntpnt,iiout,luse,.false.,.false.)

                 if (.not. luse) then
                    if(k==levs) then
                       cycle loop_readsb
                    else
                       cycle LOOP_K_LEVS
                    endif
                 endif
                 if(iiout > 0) isort(iiout)=0
                 if(ndata >  ntmp)then
                    nodata=nodata+1
                 end if
                 isort(icntpnt)=iout

              else
                 ndata=ndata+1
                 nodata=nodata+1
                 iout=ndata
                 isort(icntpnt)=iout
              endif

              if(ndata > maxobs) then
                 write(6,*)'READ_WCPBUFR:  ***WARNING*** ndata > maxobs for ',obstype
                 ndata = maxobs
              end if

!             Set usage variable              
              usage = zero


              if(icuse(nc) <= 0)usage=100._r_kind
              if(qm == 15 .or. qm == 12 .or. qm == 9)usage=100._r_kind
              if(plevs(k) < 0.0001_r_kind) then
                 write(*,*) 'warning: obs pressure is too small:',kx,k,plevs(k)
                 cycle
              endif

              if(ncnumgrp(nc)>0 )then                 ! default cross validation on
                 if(mod(ndata+1,ncnumgrp(nc))== ncgroup(nc)-1)usage=ncmiter(nc)
              end if

!             Extract pressure level and quality marks
              dlnpob=log(plevs(k))  ! ln(pressure in cb)
 
!             solid-water content path (Hurricane GPROF: TMI and GMI)
              if(swcpob) then

                 swcpoe=obserr(3,k)*convert
                 swcpmerr=swcpoe
                 cdata_all(1,iout)=swcpoe                  ! swcp error
                 cdata_all(2,iout)=dlon                    ! grid relative longitude
                 cdata_all(3,iout)=dlat                    ! grid relative latitude
                 cdata_all(4,iout)=obsdat(3,k)*convert     ! swcp obs
                 cdata_all(5,iout)=rstation_id             ! station id
                 cdata_all(6,iout)=t4dv                    ! time
                 cdata_all(7,iout)=nc                      ! type
                 cdata_all(8,iout)=swcpmerr                ! swcp max error
                 cdata_all(9,iout)=swcpq                   ! quality mark
                 cdata_all(10,iout)=swcpoe                 ! original obs error
                 cdata_all(11,iout)=usage                  ! usage parameter
                 cdata_all(12,iout)=dlon_earth_deg         ! earth relative longitude (degrees)
                 cdata_all(13,iout)=dlat_earth_deg         ! earth relative latitude (degrees)
                 cdata_all(14,iout)=stnelev                ! station elevation (m)
                 cdata_all(15,iout)=obsdat(1,k)            ! observation pressure (hPa)
                 cdata_all(16,iout)=obsdat(2,k)            ! observation height (m)
 
!             liquid-water content path (Hurricane GPROF: TMI and GMI)
              else if(lwcpob) then

                 lwcpoe=obserr(4,k)*convert
                 lwcpmerr=lwcpoe
                 cdata_all(1,iout)=lwcpoe                  ! lwcp error
                 cdata_all(2,iout)=dlon                    ! grid relative longitude
                 cdata_all(3,iout)=dlat                    ! grid relative latitude
                 cdata_all(4,iout)=obsdat(4,k)*convert     ! lwcp obs
                 cdata_all(5,iout)=rstation_id             ! station id
                 cdata_all(6,iout)=t4dv                    ! time
                 cdata_all(7,iout)=nc                      ! type
                 cdata_all(8,iout)=lwcpmerr                ! lwcp max error
                 cdata_all(9,iout)=lwcpq                   ! quality mark
                 cdata_all(10,iout)=lwcpoe                 ! original obs error
                 cdata_all(11,iout)=usage                  ! usage parameter
                 cdata_all(12,iout)=dlon_earth_deg         ! earth relative longitude (degrees)
                 cdata_all(13,iout)=dlat_earth_deg         ! earth relative latitude (degrees)
                 cdata_all(14,iout)=stnelev                ! station elevation (m)
                 cdata_all(15,iout)=obsdat(1,k)            ! observation pressure (hPa)
                 cdata_all(16,iout)=obsdat(2,k)            ! observation height (m)

              end if

           end do  LOOP_K_LEVS   !    End k loop over levs
        end do loop_readsb   !   End of bufr read loop
     enddo loop_msg

!    Close unit to bufr file
     call closbf(lunin)

!    Deallocate arrays used for thinning data
     if (.not.use_all) then
        deallocate(presl_thin)
        call del3grids
     endif

! Normal exit

  enddo loop_convinfo! loops over convinfo entry matches
  deallocate(lmsg,tab,nrep)

! Write header record and data to output file for further processing
  allocate(iloc(ndata))
  icount=0
  do i=1,maxobs
     if(isort(i) > 0)then
       icount=icount+1
       iloc(icount)=isort(i)
     end if
  end do
  if(ndata /= icount)then
     write(6,*) ' WCPBUFR: mix up in read_wcpbufr ,ndata,icount ',ndata,icount
     call stop2(50)
  end if
  allocate(cdata_out(nreal,ndata))
  do i=1,ndata
     itx=iloc(i)
     do k=1,nreal
        cdata_out(k,i)=cdata_all(k,itx)
     end do
  end do
  deallocate(iloc,isort,cdata_all)

  call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs)
  write(lunout) obstype,sis,nreal,nchanl,ilat,ilon,ndata
  write(lunout) cdata_out

  deallocate(cdata_out)

900 continue
  if(diagnostic_reg .and. ntest>0) write(6,*)'READ_WCPBUFR:  ',&
     'ntest,disterrmax=',ntest,disterrmax
  if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_WCPBUFR:  ',&
     'nvtest,vdisterrmax=',ntest,vdisterrmax

  if (ndata == 0) then 
     call closbf(lunin)
     write(6,*)'READ_WCPBUFR:  closbf(',lunin,')'
  endif

  close(lunin)

! End of routine
  return

end subroutine read_wcpbufr