program prepbufr_encode_selected_obs
!
! read all observations out from prepbufr. 
! read bufr table from prepbufr file
! write the selected obs to the new file.

 implicit none

 ! Define mnemonic (coming from GSI read_prepbufr.f90 )
 integer, parameter :: mxmn=35, mxlv=250
 character(80):: hdstr='SID XOB YOB DHR TYP ELV SAID T29'
 character(80):: obstr='POB QOB TOB ZOB UOB VOB PWO CAT PRSS'
 character(80):: qcstr='PQM QQM TQM ZQM WQM NUL PWQ     '
 character(80):: oestr='POE QOE TOE NUL WOE NUL PWE     '
  
 real(8) :: hdr(mxmn),obs(mxmn,mxlv),qcf(mxmn,mxlv),oer(mxmn,mxlv)
 real(8) :: rstation_id
 real    :: lat,lon
 
 integer        :: ireadmg,ireadsb
 integer        :: unit_in=10,unit_table=24,unit_out=20,idate,nmsg,ntb
 integer        :: i,k,iret, nlvl,report_type

 character(8)   :: c_sid, msgtype


 equivalence(rstation_id,c_sid) ! convert rstation_id to character c_sid

 open(unit_table,file='prepobs_prep.bufrtable')
 open(unit_in,file='prepbufr',form='unformatted',status='old')
 open(unit_out,file='prepbufr_selected_obs',action='write',form='unformatted')


 call openbf(unit_in,'IN',unit_in)
 call dxdump(unit_in,unit_table) ! dump bufr table from existing prepbufr file
 call openbf(unit_out,'OUT',unit_table) ! connect the bufr table to bufr lib,  
                                        ! it will write bufr table to new file. 

 call datelen(10)
 nmsg=0
 msg_report: do while (ireadmg(unit_in,msgtype,idate) == 0) ! Read message type from existing file
    call openmb(unit_out,msgtype,idate) ! Opens new message type, will write it to new file
    nmsg=nmsg+1
    ntb = 0
    sb_report: do while (ireadsb(unit_in) == 0) ! Read data subset in the message type of existing file
       ntb = ntb+1
       call ufbint(unit_in,hdr,mxmn,1   ,iret,hdstr)
       call ufbint(unit_in,obs,mxmn,mxlv,iret,obstr)
       call ufbint(unit_in,oer,mxmn,mxlv,iret,oestr)
       call ufbint(unit_in,qcf,mxmn,mxlv,iret,qcstr)
       nlvl = iret ! Actually returned number of pressure level
       lon = hdr(2)
       lat = hdr(3)
       report_type = hdr(5) 

! add conditions to select data, only write the selected data to the new file.
       if ( report_type .eq. 120 .and. (lat < 5.0 .or.  lat > 20.) ) then
          write(6,*)  'msgtype=',msgtype,'   cycle time =',idate
          write(6,*) iret,c_sid,(hdr(i),i=1,8)
          do k=1,iret
            write(6,'(i3,a10,9f14.1)') k,'obs=',(obs(i,k),i=1,9)
            write(6,'(i3,a10,9f14.1)') k,'oer=',(oer(i,k),i=1,7)
            write(6,'(i3,a10,9f14.1)') k,'qcf=',(qcf(i,k),i=1,7)
          enddo

! Now write the selected obs to new file
         call ufbint(unit_out,hdr,mxmn,1   ,iret,hdstr)
         call ufbint(unit_out,obs,mxmn,nlvl,iret,obstr)
         call ufbint(unit_out,oer,mxmn,nlvl,iret,oestr)
         call ufbint(unit_out,qcf,mxmn,nlvl,iret,qcstr)
         call writsb(unit_out)  
       end if 

    enddo sb_report
    call closmg(unit_out) 
 enddo msg_report

 call closbf(unit_in)
 call closbf(unit_out)

end program