program read_gps
! subprogram: read_gps                   read in and reformat gps data

  use kinds, only: r_kind,i_kind,r_double
  implicit none

! Declare passed variables
  character(len=80) :: infile

! Declare local parameters  
  integer(i_kind),parameter:: maxlevs=500

! Declare local variables
  
  character(10) nemo
  character(80) hdr1a
  character,dimension(8):: subset
  
  integer(i_kind) lnbufr,i,k,m,maxobs,ireadmg,ireadsb,said,ptid
  integer(i_kind) idate
  integer(i_kind) iret,levs,levsr,nreps_ROSEQ1
  integer(i_kind),parameter:: mxib=31
!  integer(i_kind) ibit(mxib),nib
!  logical lone
  logical ref_obs

  real(r_kind) qfro
 
  integer(i_kind),parameter:: n1ahdr=10
  real(r_double),dimension(n1ahdr):: bfr1ahdr
  real(r_double),dimension(50,maxlevs):: data1b
  real(r_double),dimension(50,maxlevs):: data2a
  real(r_double),dimension(maxlevs):: nreps_this_ROSEQ2
 
  data lnbufr/10/
  data hdr1a / 'YEAR MNTH DAYS HOUR MINU PCCF ELRC SAID PTID GEODU' / 
  data nemo /'QFRO'/
  
!***********************************************************************************

  infile='gpsbufr'
  ref_obs=.true.

! Open file for input, then read bufr data
  open(lnbufr,file=infile,form='unformatted')
  call openbf(lnbufr,'IN',lnbufr)
  call datelen(10)
  call readmg(lnbufr,subset,idate,iret)
  if (iret/=0) goto 1010


! Big loop over the bufr file

  do while(ireadmg(lnbufr,subset,idate)==0)
     read_loop:  do while(ireadsb(lnbufr)==0)

! Read/decode data in subset (profile)

! Extract header information
        call ufbint(lnbufr,bfr1ahdr,n1ahdr,1,iret,hdr1a)
        call ufbint(lnbufr,qfro,1,1,iret,nemo)

!        if (said == 4) then ! Gras
!           call upftbv(lnbufr,nemo,qfro,mxib,ibit,nib)
!
!           if(lone) then
!              write(6,*)'READ_GPS:  bad profile said=',said,'ptid=',ptid,&
!                   ' SKIP this report'
!              cycle read_loop
!           endif
!        endif


! Read bending angle information
! Get the number of occurences of sequence ROSEQ2 in this subset
! (will also be the number of replications of sequence ROSEQ1), nreps_ROSEQ1
! Also determine the number of replications of sequence ROSEQ2 nested inside
! each replication of ROSEQ1,
! nreps_this_ROSEQ2(1:nreps_ROSEQ1) - currently = 3 frequencies (L1, L2, zero)

        call ufbint(lnbufr,nreps_this_ROSEQ2,1,maxlevs,nreps_ROSEQ1,'{ROSEQ2}')

! Store entire contents of ROSEQ1 sequence (including contents of nested ROSEQ2 sequence)
! in array data1b
        call ufbseq(lnbufr,data1b,50,maxlevs,levs,'ROSEQ1') 

        if(levs.ne.nreps_ROSEQ1) then
           write(6,*) 'READ_GPS:  **WARNING** said,ptid=',said,ptid,&
                ' mismatch between sequence of ROSEQ1 and ROSEQ2 occurence',levs,nreps_ROSEQ1, &
                ' SKIP this report'
           cycle read_loop
        endif

! Check we have the same number of levels for ref and bending angle
! when ref_obs on to get lat/lon information

        call ufbseq(lnbufr,data2a,50,maxlevs,levsr,'ROSEQ3') ! refractivity
        if ((ref_obs).and.(levs/=levsr)) then
           write(6,*) 'READ_GPS:  **WARNING** said,ptid=',said,ptid,&
                ' with gps_bnd levs=',levs,&
                ' and gps_ref levsr=',levsr,&
                ' SKIP this report'
           cycle read_loop
        endif

     enddo read_loop        ! subsets
  enddo                     ! messages

! Close unit to input file
1010 continue
  call closbf(lnbufr)

end program read_gps