module sfcobsqc
!$$$ module documentation block
!           .      .    .                                       .
! module:   sfcobsqc
!   prgmmr: pondeca
!
! abstract: contains subroutines for the qc of surface obs based
!           on (i) the provider uselist for mesonet winds, (ii) 
!           station uselist for mesonet winds, and (iii) rejectlists 
!           for any ob type (u and v-wind, wind speed, temperatature, 
!           moisture, surface pressure). the code inquires for the 
!           existence of thsese lists in the gsi working directorty 
!           and applies them if present.
!           
!
! program history log:
!   2007-10-19  pondeca
!   2011-02-15  zhu - add get_gustqm
!   2014-04-10  pondeca - add reject lists for td, mxtm, mitm, pmsl
!   2014-05-07  pondeca - add reject list for howv
!   2014-07-11  carley - add reject list for lcbas and tcamt
!   2014-10-01  Xue - add GSD surface data uselist
!   2015-07-10  pondeca - add reject list for cldch
!   2018-03-14  pondeca - add station accept list for mesonet visibility
!
! subroutines included:
!   sub init_rjlists
!   sub get_usagerj
!   sub get_gustqm
!   readin_rjlists
!   sub destroy_rjlists
!
! variable definitions:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use kinds, only: i_kind,r_kind,r_single,r_double

  implicit none

  private

  character(16),allocatable,dimension(:)::cprovider
  character(5),allocatable,dimension(:)::csta_winduse
  character(80),allocatable,dimension(:)::w_rjlist,t_rjlist,p_rjlist,q_rjlist
  character(80),allocatable,dimension(:)::td_rjlist,mxtm_rjlist,mitm_rjlist,pmsl_rjlist,howv_rjlist, &
                                          lcbas_rjlist,tcamt_rjlist,cldch_rjlist
  character(80),allocatable,dimension(:)::t_day_rjlist,t_night_rjlist
  character(80),allocatable,dimension(:)::q_day_rjlist,q_night_rjlist
  character(8),allocatable,dimension(:,:)::csta_windbin
  character(80),allocatable,dimension(:)::csta_visuse

  integer(i_kind) sfcuselist_nt_use
  character(8),allocatable,dimension(:)::sfcuselist_use_id
  character(1),allocatable,dimension(:)::w_use_sfcuselist
  character(1),allocatable,dimension(:)::t_use_sfcuselist
  character(1),allocatable,dimension(:)::td_use_sfcuselist

  integer(i_kind) nprov,nwrjs,ntrjs,nprjs,nqrjs,nsta_mesowind_use,nsta_mesovis_use
  integer(i_kind) ntdrjs,nmxtmrjs,nmitmrjs,npmslrjs,nhowvrjs,&
                  nlcbasrjs,ntcamtrjs,ncldchrjs
  integer(i_kind) ntrjs_day,ntrjs_night
  integer(i_kind) nqrjs_day,nqrjs_night
  integer(i_kind) nbins
  integer(i_kind),allocatable,dimension(:)::nwbaccpts

  logical gsdsfclistexist
  logical gsdsfcproviderlistexist
  logical listexist
  logical wlistexist
  logical tlistexist
  logical plistexist
  logical qlistexist
  logical tdlistexist
  logical mxtmlistexist
  logical mitmlistexist
  logical pmsllistexist
  logical howvlistexist
  logical lcbaslistexist
  logical tcamtlistexist
  logical cldchlistexist
  logical listexist2
  logical t_day_listexist
  logical t_night_listexist
  logical q_day_listexist
  logical q_night_listexist
  logical wbinlistexist
  logical vis_uselistexist

  public init_rjlists
  public get_usagerj
  public get_gustqm
  public readin_rjlists
  public get_sunangle
  public get_wbinid
  public destroy_rjlists

  public init_gsd_sfcuselist
  public apply_gsd_sfcuselist
  public destroy_gsd_sfcuselist

  logical :: verbose = .false.
contains

subroutine init_gsd_sfcuselist
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_gsd_sfcuselist
!   prgmmr:
!
! abstract: read in GSD surface observation uselist 
!
! program history log:
!   2014-10-01  Xue - initial code
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use mpimod, only: mype

  implicit none

  integer(i_kind) use_unit
  character(150) cstring
  character(80) clistname

  integer(i_kind), parameter::nmax=60000_i_kind

  data use_unit / 777_i_kind /
!**************************************************************************
  if (mype==0) verbose =.false.
  sfcuselist_nt_use= 0
  allocate(sfcuselist_use_id(nmax))
  allocate(w_use_sfcuselist(nmax))
  allocate(t_use_sfcuselist(nmax))
  allocate(td_use_sfcuselist(nmax))

  gsdsfclistexist=.false.

  inquire(file='gsd_sfcobs_uselist.txt',exist=gsdsfclistexist)
  if(gsdsfclistexist) then
    open (use_unit,file='gsd_sfcobs_uselist.txt',form='formatted')

    read_use: do
      read(use_unit,'(a150)',end=7745) cstring
      if(cstring(1:1) == ';') cycle read_use    ! skip comments marked as ;

      sfcuselist_nt_use=sfcuselist_nt_use+1
      sfcuselist_use_id(sfcuselist_nt_use)= adjustl(cstring(1:10))
      w_use_sfcuselist(sfcuselist_nt_use)= adjustl(cstring(11:12))
      t_use_sfcuselist(sfcuselist_nt_use)= adjustl(cstring(13:14))
      td_use_sfcuselist(sfcuselist_nt_use)= adjustl(cstring(15:16))
      if(verbose) print*,'sfcuselist_use_id=',sfcuselist_nt_use,&
                                    sfcuselist_use_id(sfcuselist_nt_use),&
                                ",",w_use_sfcuselist(sfcuselist_nt_use),&
                                ",",t_use_sfcuselist(sfcuselist_nt_use),&
                                ",",t_use_sfcuselist(sfcuselist_nt_use)

    end do read_use
7745 continue
  endif
  close(use_unit)

!==> Read mesonet provider names from the uselist
  clistname='gsd_sfcobs_provider.txt'
  call readin_rjlists(clistname,gsdsfcproviderlistexist,cprovider,500,nprov)
  if(verbose)&
    print*,'mesonetproviderlist: provider,nprov=',gsdsfcproviderlistexist,nprov


end subroutine init_gsd_sfcuselist

subroutine destroy_gsd_sfcuselist
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    destroy_gsd_sfcuselist
!   prgmmr:
!
! abstract:
!
! program history log:
!   2015-02-05  Hu - added subprogram doc block
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
  implicit none

  deallocate(sfcuselist_use_id)
  deallocate(w_use_sfcuselist)
  deallocate(t_use_sfcuselist)
  deallocate(td_use_sfcuselist)

end subroutine destroy_gsd_sfcuselist

subroutine apply_gsd_sfcuselist(kx,obstype,c_station_id,c_prvstg,c_sprvstg, &
                       usage_rj)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    apply_gsd_sfcuselist
!   prgmmr:
!
! abstract: use GSD surface observation uselist  to decide
!           which surface observation should be used in the analysis
!
! program history log:
!   2014-10-01  Xue - initial code
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block



  use constants, only: zero_single
  implicit none

  integer(i_kind),intent(in   ) :: kx
  character(10)  ,intent(in   ) :: obstype
  character(8)   ,intent(in   ) :: c_station_id
  character(8)   ,intent(in   ) :: c_prvstg,c_sprvstg
  real(r_kind)   ,intent(inout) :: usage_rj

! Declare local variables
  integer(i_kind) m,nlen
  character(8)  ch8
  real(r_kind) usage_rj0

! Declare local parameters
  real(r_kind),parameter:: r6    = 6.0_r_kind
  real(r_kind),parameter:: r5000 = 5000._r_kind
  real(r_kind),parameter:: r5100 = 5100._r_kind
  real(r_kind),parameter:: r6000 = 6000._r_kind
  real(r_kind),parameter:: r6100 = 6100._r_kind
  real(r_kind),parameter:: r6200 = 6200._r_kind

  if (usage_rj >= r6) return  

  usage_rj0=usage_rj
  usage_rj=r6000

  if(gsdsfcproviderlistexist) then
     do m=1,nprov
        if (c_prvstg(1:8) == cprovider(m)(1:8) .and. &
           (c_sprvstg(1:8) == cprovider(m)(9:16) .or. &
            cprovider(m)(9:16) == 'allsprvs') ) then
              usage_rj=usage_rj0
              exit
        endif
     enddo
  endif

  if(.not.gsdsfclistexist) return 
  if (usage_rj==usage_rj0) return  ! station is in provider list, use it
                                   ! if not (usage_rj=r6000), check uselist

  if (kx<200 .and. (obstype=='t' .or. obstype=='q')) then  !<==mass obs
! assume everything is reject first, get back those obs exist in uselist 
    usage_rj= r5000                                           
    if(obstype=='t') then
        do m=1,sfcuselist_nt_use
           ch8(1:8)=sfcuselist_use_id(m)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then
             if (t_use_sfcuselist(m)=='0') then  ! put it to reject list
                usage_rj=r5000
             elseif(t_use_sfcuselist(m)=='1') then ! use it original usage value
                usage_rj= usage_rj0                                          
             endif
             exit
           endif
        enddo
     elseif(obstype=='q') then
        do m=1,sfcuselist_nt_use
           ch8(1:8)=sfcuselist_use_id(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then
              if (td_use_sfcuselist(m)=='0') then  ! put it to reject list
                 usage_rj=r5000
             elseif(td_use_sfcuselist(m)=='1') then ! use it original usage value
                usage_rj= usage_rj0                                          
             endif
             exit
           endif
        enddo
    endif
   endif ! kx < 200

   if (kx>=200) then  ! wind vector obs
     if(obstype=='uv') then
        usage_rj= r5000                                           
        do m=1,sfcuselist_nt_use
           ch8(1:8)=sfcuselist_use_id(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==288.or.kx==295).and. c_station_id(1:nlen)==ch8(1:nlen))) then
              if (w_use_sfcuselist(m)=='0') then  ! put it to reject list
                 usage_rj=r5000
             elseif(w_use_sfcuselist(m)=='1') then ! use it original usage value
                usage_rj= usage_rj0                                          
             endif
             exit
           endif
        enddo
     endif
   end if

end subroutine apply_gsd_sfcuselist

subroutine init_rjlists
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_rjlists
!   prgmmr:
!
! abstract: initialize qc lists 
!
! program history log:
!   2009-10-01  lueken - added subprogram doc block
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use mpimod, only: mype

  implicit none

! Declare passed variables

! Declare local variables
  integer(i_kind) meso_unit,ncount
  integer(i_kind) ibin
  real(r_single) aa1,aa2
  character(80) cstring
  character(80) clistname
  character(8) ach8,bch8
  character(16) ch16

  integer(i_kind), parameter::nmax=100000
  integer(i_kind) nmax2

  data meso_unit / 20 /
!**************************************************************************
  if (mype==0) verbose =.true.

  allocate(cprovider(500))
  allocate(w_rjlist(nmax))
  allocate(t_rjlist(nmax))
  allocate(p_rjlist(nmax))
  allocate(q_rjlist(nmax))
  allocate(td_rjlist(nmax))
  allocate(mxtm_rjlist(nmax))
  allocate(mitm_rjlist(nmax))
  allocate(pmsl_rjlist(nmax))
  allocate(howv_rjlist(nmax))
  allocate(lcbas_rjlist(nmax))
  allocate(tcamt_rjlist(nmax))
  allocate(cldch_rjlist(nmax))
  allocate(csta_winduse(nmax))
  allocate(t_day_rjlist(nmax))
  allocate(t_night_rjlist(nmax))
  allocate(q_day_rjlist(nmax))
  allocate(q_night_rjlist(nmax))
  allocate(csta_visuse(nmax))

!==> Read mesonet provider names from the uselist
 clistname='mesonetuselist'
 call readin_rjlists(clistname,listexist,cprovider,500,nprov)
 if(verbose)&
 print*,'mesonetuselist: listexist,nprov=',listexist,nprov


 
!==> Read in station names from the reject list for wind
 clistname='w_rejectlist'
 call readin_rjlists(clistname,wlistexist,w_rjlist,nmax,nwrjs)
 if(verbose)&
 print*,'w_rejectlist: wlistexist,nwrjs=',wlistexist,nwrjs


 
!==> Read in station names from the reject lists for temperature
 clistname='t_rejectlist'
 call readin_rjlists(clistname,tlistexist,t_rjlist,nmax,ntrjs)
 if(verbose)&
 print*,'t_rejectlist: tlistexist,ntrjs=',tlistexist,ntrjs


 
 clistname='t_day_rejectlist'
 call readin_rjlists(clistname,t_day_listexist,t_day_rjlist,nmax,ntrjs_day)
 if(verbose)&
 print*,'t_day_rejectlist: t_day_listexist,ntrjs_day=',t_day_listexist,ntrjs_day


 
 clistname='t_night_rejectlist'
 call readin_rjlists(clistname,t_night_listexist,t_night_rjlist,nmax,ntrjs_night)
 if(verbose)&
 print*,'t_night_rejectlist: t_night_listexist,ntrjs_night=',t_night_listexist,ntrjs_night


 
!==> Read in station names from the reject list for surface pressure
 clistname='p_rejectlist'
 call readin_rjlists(clistname,plistexist,p_rjlist,nmax,nprjs)
 if(verbose)&
 print*,'p_rejectlist: plistexist,nprjs=',plistexist,nprjs


 
!==> Read in station names from the reject lists for specific humidity
 clistname='q_rejectlist'
 call readin_rjlists(clistname,qlistexist,q_rjlist,nmax,nqrjs)
 if(verbose)&
 print*,'q_rejectlist: qlistexist,nqrjs=',qlistexist,nqrjs


 
 clistname='q_day_rejectlist'
 call readin_rjlists(clistname,q_day_listexist,q_day_rjlist,nmax,nqrjs_day)
 if(verbose)&
 print*,'q_day_rejectlist: q_day_listexist,nqrjs_day=',q_day_listexist,nqrjs_day

  

 clistname='q_night_rejectlist'
 call readin_rjlists(clistname,q_night_listexist,q_night_rjlist,nmax,nqrjs_night)
 if(verbose)&
 print*,'q_night_rejectlist: q_night_listexist,nqrjs_night=',q_night_listexist,nqrjs_night



 clistname='td_rejectlist'
 call readin_rjlists(clistname,tdlistexist,td_rjlist,nmax,ntdrjs)
 if(verbose)&
 print*,'td_rejectlist: tdlistexist,ntdrjs=',tdlistexist,ntdrjs

 clistname='mxtm_rejectlist'
 call readin_rjlists(clistname,mxtmlistexist,mxtm_rjlist,nmax,nmxtmrjs)
 if(verbose)&
 print*,'mxtm_rejectlist: mxtmlistexist,nmxtmrjs=',mxtmlistexist,nmxtmrjs

 clistname='mitm_rejectlist'
 call readin_rjlists(clistname,mitmlistexist,mitm_rjlist,nmax,nmitmrjs)
 if(verbose)&
 print*,'mitm_rejectlist: mitmlistexist,nmitmrjs=',mitmlistexist,nmitmrjs


 clistname='pmsl_rejectlist'
 call readin_rjlists(clistname,pmsllistexist,pmsl_rjlist,nmax,npmslrjs)
 if(verbose)&
 print*,'pmsl_rejectlist: pmsllistexist,npmslrjs=',pmsllistexist,npmslrjs


 clistname='howv_rejectlist'
 call readin_rjlists(clistname,howvlistexist,howv_rjlist,nmax,nhowvrjs)
 if(verbose)&
 print*,'howv_rejectlist: howvlistexist,nhowvrjs=',howvlistexist,nhowvrjs
 

 clistname='lcbas_rejectlist'
 call readin_rjlists(clistname,lcbaslistexist,lcbas_rjlist,nmax,nlcbasrjs)
 if(verbose)&
 print*,'lcbas_rejectlist: lcbaslistexist,nlcbasrjs=',lcbaslistexist,nlcbasrjs

 clistname='tcamt_rejectlist'
 call readin_rjlists(clistname,tcamtlistexist,tcamt_rjlist,nmax,ntcamtrjs)
 if(verbose)&
 print*,'tcamt_rejectlist: tcamtlistexist,ntcamtrjs=',tcamtlistexist,ntcamtrjs

 clistname='cldch_rejectlist'
 call readin_rjlists(clistname,cldchlistexist,cldch_rjlist,nmax,ncldchrjs)
 if(verbose)&
 print*,'cldch_rejectlist: cldchlistexist,ncldchrjs=',cldchlistexist,ncldchrjs

!==> Read in 'good' mesonet station names from the station uselist
 inquire(file='mesonet_stnuselist',exist=listexist2)
 if(listexist2) then
    open (meso_unit,file='mesonet_stnuselist',form='formatted')
    ncount=0
    do
       ncount=ncount+1
       read(meso_unit,'(a5,a80)',end=181) csta_winduse(ncount),cstring
    end do
181 continue
    nsta_mesowind_use=ncount-1
    if(verbose)&
    print*,'mesonet_stnuselist: nsta_mesowind_use=',nsta_mesowind_use
 endif
 close(meso_unit)

!==> Read in wind direction stratified wind accept lists
 inquire(file='wbinuselist',exist=wbinlistexist)
 if(verbose)&
 print*,'wdirbinuselist: wbinuselist=',wbinlistexist

 nbins=0
 if(wbinlistexist) then
    open (meso_unit,file='wbinuselist',form='formatted')
    read(meso_unit,'(a16,i2)',end=191) ch16,nbins
    allocate(nwbaccpts(max(1,nbins)))
    nwbaccpts(:)=0
    do
       read(meso_unit,'(a8,3x,a8,3x,f7.4,2x,f9.4,3x,i2)',end=191) ach8,bch8,aa1,aa2,ibin
       nwbaccpts(ibin)=nwbaccpts(ibin)+1
    end do
191 continue
    if(verbose)then
       print*,'wdirbinuselist: number of bins=',nbins
       print*,'wdirbinuselist: (nwbaccpts(ibin),ibin=1,nbins)=',(nwbaccpts(ibin),ibin=1,nbins)
    endif

    nmax2=maxval(nwbaccpts(:))
    if(verbose)&
    print*,'wdirbinuselist: maximum number of obs in a single bin=',nmax2

    allocate(csta_windbin(max(1,nmax2),max(1,nbins)))

    rewind(meso_unit)
    read(meso_unit,'(a16,i2)',end=193) ch16,nbins
    nwbaccpts(:)=0
    do
       read(meso_unit,'(a8,3x,a8,3x,f7.4,2x,f9.4,3x,i2)',end=193) ach8,bch8,aa1,aa2,ibin
       nwbaccpts(ibin)=nwbaccpts(ibin)+1
       csta_windbin(nwbaccpts(ibin),ibin)=ach8
    end do
193 continue
 endif
 close(meso_unit)

!==> Read in 'good' mesonet station names from the station uselist for visibility
 inquire(file='mesonet_stnuselist_for_vis',exist=vis_uselistexist)
 if(vis_uselistexist) then
    open (meso_unit,file='mesonet_stnuselist_for_vis',form='formatted')
    ncount=0
    do
      ncount=ncount+1
      read(meso_unit,'(a5,a80)',end=194) csta_visuse(ncount),cstring
    end do
194 continue
    nsta_mesovis_use=ncount-1
!    if(verbose)&
    print*,'mesonet_stnuselist_for_vis: nsta_mesovis_use=',nsta_mesovis_use
 endif
 close(meso_unit)

end subroutine init_rjlists

subroutine get_usagerj(kx,obstype,c_station_id,c_prvstg,c_sprvstg, &
                       dlon,dlat,idate,dtime,udbl,vdbl,usage_rj)

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    get_usagerg
!   prgmmr:
!
! abstract: determine the usage value of read_prepbufr for surface obs. the following
!           is done: (i) if incoming usage value is >=100. then do nothing, since
!           read_prepbufr has already flagged this ob and assigned a specific usage 
!           value to it. (ii) use usage=500. for temperature, moisture, or surface pressure
!           obs which are found in the rejectlist. (iii) 
!
! program history log:
!   2009-10-01  lueken - added subprogram doc block
!
!   input argument list:
!    kx
!    obstype
!    c_station_id
!    c_prvstg,c_sprvstg
!    dlon
!    dlat
!    idate
!    udbl
!    vdbl
!
!   output argument list:
!    usage_rj
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use constants, only: zero_single
  use gridmod, only: twodvar_regional,tll2xy
  use ndfdgrids, only: valley_adjustment

  implicit none

  integer(i_kind),intent(in   ) :: kx
  integer(i_kind),intent(in   ) :: idate
  character(10)  ,intent(in   ) :: obstype
  character(8)   ,intent(in   ) :: c_station_id
  character(8)   ,intent(in   ) :: c_prvstg,c_sprvstg
  real(r_kind)   ,intent(in   ) :: dlon
  real(r_kind)   ,intent(in   ) :: dlat
  real(r_kind)   ,intent(in   ) :: dtime
  real(r_double) ,intent(in   ) :: udbl,vdbl
  real(r_kind)   ,intent(inout) :: usage_rj

! Declare local variables
  integer(i_kind) m,nlen
  integer(i_kind) ibin
  character(8)  ch8
  real(r_kind) usage_rj0
  real(r_single) sunangle
  real(r_kind) xob,yob
  logical outside

! Declare local parameters
  real(r_kind),parameter:: r6    = 6.0_r_kind
  real(r_kind),parameter:: r5000 = 5000._r_kind
  real(r_kind),parameter:: r5100 = 5100._r_kind
  real(r_kind),parameter:: r6000 = 6000._r_kind
  real(r_kind),parameter:: r6100 = 6100._r_kind
  real(r_kind),parameter:: r6200 = 6200._r_kind

  if (usage_rj >= r6) return

  usage_rj0=usage_rj

  if (kx<200) then  !<==mass obs

    if(obstype=='t') then

      if(tlistexist ) then
        do m=1,ntrjs
           ch8(1:8)=t_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then !handle wfo's mesonets which never end with
              usage_rj=r5000                                           !an "a" or "x" in the eight position following blanks
              exit
           endif
        enddo
      endif

      if((t_day_listexist .or. t_night_listexist) .and. usage_rj==usage_rj0) then
           call get_sunangle(idate,dtime,dlon,dlat,sunangle)
           if (sunangle > zero_single) then
              do m=1,ntrjs_day
                 ch8(1:8)=t_day_rjlist(m)(1:8)
                 nlen=len_trim(ch8)
                 if ((trim(c_station_id) == trim(ch8)) .or. &
                     ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then
                    usage_rj=r5100
                    exit
                 endif
              enddo
             else 
              do m=1,ntrjs_night
                 ch8(1:8)=t_night_rjlist(m)(1:8)
                 nlen=len_trim(ch8)
                 if ((trim(c_station_id) == trim(ch8)) .or. &
                     ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then
                    usage_rj=r5100
                    exit
                 endif
              enddo
           endif 
      endif

     elseif(obstype=='q') then

      if(qlistexist ) then
        do m=1,nqrjs
           ch8(1:8)=q_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then
              usage_rj=r5000
              exit
           endif
        enddo
      endif

      if((q_day_listexist .or. q_night_listexist) .and. usage_rj==usage_rj0) then
           call get_sunangle(idate,dtime,dlon,dlat,sunangle)
           if (sunangle > zero_single) then
              do m=1,nqrjs_day
                 ch8(1:8)=q_day_rjlist(m)(1:8)
                 nlen=len_trim(ch8)
                 if ((trim(c_station_id) == trim(ch8)) .or. &
                     ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then
                    usage_rj=r5100
                    exit
                 endif
              enddo
             else 
              do m=1,nqrjs_night
                 ch8(1:8)=q_night_rjlist(m)(1:8)
                 nlen=len_trim(ch8)
                 if ((trim(c_station_id) == trim(ch8)) .or. &
                     ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then
                    usage_rj=r5100
                    exit
                 endif
              enddo
           endif 
      endif

     elseif(obstype=='ps' .and. plistexist ) then
        do m=1,nprjs
           ch8(1:8)=p_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then
              usage_rj=r5000
              exit
           endif
        enddo

     elseif(obstype=='td2m' .and. tdlistexist ) then
        do m=1,ntdrjs
           ch8(1:8)=td_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then
              usage_rj=r5000
              exit
           endif
        enddo

     elseif(obstype=='mxtm' .and. mxtmlistexist) then
        do m=1,nmxtmrjs
           ch8(1:8)=mxtm_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then !handle wfo's mesonets which never end with
              usage_rj=r5000                                           !an "a" or "x" in the eight position following blanks
              exit
           endif
        enddo

     elseif(obstype=='mitm' .and. mitmlistexist) then
        do m=1,nmitmrjs
           ch8(1:8)=mitm_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then !handle wfo's mesonets which never end with
              usage_rj=r5000                                           !an "a" or "x" in the eight position following blanks
              exit
           endif
        enddo

     elseif(obstype=='pmsl' .and. pmsllistexist) then
        do m=1,npmslrjs
           ch8(1:8)=pmsl_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then !handle wfo's mesonets which never end with
              usage_rj=r5000                                           !an "a" or "x" in the eight position following blanks
              exit
           endif
        enddo

     elseif(obstype=='howv' .and. howvlistexist) then
        do m=1,nhowvrjs
           ch8(1:8)=howv_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then !handle wfo's mesonets which never end with
              usage_rj=r5000                                           !an "a" or "x" in the eight position following blanks
              exit
           endif
        enddo

     elseif(obstype=='lcbas' .and. lcbaslistexist) then
        do m=1,nlcbasrjs
           ch8(1:8)=lcbas_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then !handle wfo's mesonets which never end with
              usage_rj=r5000                                           !an "a" or "x" in the eight position following blanks
              exit
           endif
        enddo

     elseif(obstype=='tcamt' .and. tcamtlistexist) then
        do m=1,ntcamtrjs
           ch8(1:8)=tcamt_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then !handle wfo's mesonets which never end with
              usage_rj=r5000                                           !an "a" or "x" in the eight position following blanks
              exit
           endif
        enddo

     elseif(obstype=='cldch' .and. cldchlistexist) then
        do m=1,ncldchrjs
           ch8(1:8)=cldch_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
               ((kx==188.or.kx==195).and.c_station_id(1:nlen)==ch8(1:nlen))) then !handle wfo's mesonets which never end with
              usage_rj=r5000                                           !an "a" or "x" in the eight position following blanks
              exit
           endif
        enddo

     end if
  endif

  if (kx>=200) then !<==wind obs

     if (kx==288.or.kx==295) then
        usage_rj=r6000
        if (listexist) then !note that uselists must precede the rejectlist
           do m=1,nprov
              if (cprovider(m)(1:7)=='allprvs' .or. & 
                 (c_prvstg(1:8) == cprovider(m)(1:8) .and. (c_sprvstg(1:8) == cprovider(m)(9:16)  &
                                                      .or. cprovider(m)(9:16) == 'allsprvs')) ) then
                 usage_rj=usage_rj0
                 exit
              endif
           enddo
        endif
        if (listexist2 .and. usage_rj/=usage_rj0) then
           do m=1,nsta_mesowind_use
              if (c_station_id(1:5) == csta_winduse(m)(1:5)) then
                 usage_rj=usage_rj0
                 exit
              endif
           enddo
        endif
        if(wbinlistexist .and. usage_rj/=usage_rj0) then
          call get_wbinid(udbl,vdbl,nbins,ibin)
          if(ibin <= nbins ) then
            do m=1,nwbaccpts(ibin)
              ch8(1:8)=csta_windbin(m,ibin)(1:8)
              nlen=len_trim(ch8)
              if (c_station_id(1:nlen)==ch8(1:nlen)) then
                usage_rj=usage_rj0
                exit
              endif
            enddo
          endif
        endif
     endif

     if( (obstype=='uv' .or. obstype=='wspd10m' .or. obstype=='uwnd10m' .or.  obstype=='vwnd10m') .and. wlistexist ) then
        do m=1,nwrjs
           ch8(1:8)=w_rjlist(m)(1:8)
           nlen=len_trim(ch8)
           if ((trim(c_station_id) == trim(ch8)) .or. &
              ((kx==288.or.kx==295).and. c_station_id(1:nlen)==ch8(1:nlen))) then
              if (kx/=288.and.kx/=295) then
                 usage_rj=r5000
              else
                 if (usage_rj==usage_rj0)    usage_rj=r6100 !ob is in at least one of the above three uselists
                 if (usage_rj==r6000)        usage_rj=r6200 !ob is in none of the above three uselists
              endif
              exit
           endif
        enddo
     endif

  end if

!==> station uselist for mesonet visibility

  if (obstype=='vis' .and.( kx==188.or.kx==288.or.kx==195.or.kx==295) )then
     usage_rj=r6000
     if (vis_uselistexist .and. usage_rj/=usage_rj0) then !note that usage_rj could differ from usage_rj0 after the rejectlist application
        do m=1,nsta_mesovis_use                          !which happens to be currently unavailable for vis
           nlen=len_trim(csta_visuse(m))
           if (c_station_id(1:nlen) == csta_visuse(m)(1:nlen)) then
              usage_rj=usage_rj0
              exit
           endif
        enddo
     endif
  end if

  if (twodvar_regional) then
     call tll2xy(dlon,dlat,xob,yob,outside)
     if ((obstype=='t' .or. obstype=='q') .and. .not.outside) call valley_adjustment(xob,yob,usage_rj)
  endif
end subroutine get_usagerj


subroutine get_gustqm(kx,c_station_id,c_prvstg,c_sprvstg,gustqm)
!$$$  module documentation block
! abstract: determine the qm for gust 1) for gust that is not on the uselist qm=9
!                                     2) for gust that is on the rejectlist qm=3 or 15
! program history log:
!   2009-01-29  zhu
!
  implicit none

  integer(i_kind) kx
  integer(i_kind) nlen
  integer(i_kind) gustqm
  character(8),intent(in)::  c_station_id
  character(8),intent(in)::  c_prvstg,c_sprvstg

! Declare local variables
  integer(i_kind) m
  character(8)  ch8

  gustqm=9
  if (listexist) then
     do m=1,nprov
        if (cprovider(m)(1:7)=='allprvs' .or. &
           (c_prvstg(1:8) == cprovider(m)(1:8) .and. (c_sprvstg(1:8) == cprovider(m)(9:16)  &
                                               .or. cprovider(m)(9:16) == 'allsprvs')) ) then
           gustqm=0
           exit
         endif
     enddo
  endif
  if (listexist2) then
     do m=1,nsta_mesowind_use
        if (c_station_id(1:5) == csta_winduse(m)(1:5)) then
           gustqm=0
           exit
         endif
     enddo
  endif

  if(wlistexist ) then
     do m=1,nwrjs
        ch8(1:8)=w_rjlist(m)(1:8)
        nlen=len_trim(ch8)
        if ((trim(c_station_id) == trim(ch8)) .or. &
            ((kx==288.or.kx==295).and.c_station_id(1:nlen)==ch8(1:nlen))) then
            if (gustqm==0) then
               gustqm=3
            else
               gustqm=15
            end if
            exit
        endif
     enddo
  endif
end subroutine get_gustqm


subroutine readin_rjlists(clistname,fexist,clist,ndim,ncount)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    readin_rjlists
!   prgmmr:
!
! abstract: read in accept or reject list
!
! program history log:
!   2012-17-02  pondeca
!
!   input argument list:
!   clistname
!   ndim
!
!   output argument list:
!   clist
!   ncount
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  implicit none

! Declare passed variables
  integer(i_kind),intent(in):: ndim
  character(80),intent(in):: clistname

  integer(i_kind),intent(out):: ncount
  character(*),intent(out):: clist(ndim)
  logical,intent(out)::fexist

! Declare local variables
  integer(i_kind) meso_unit,n,m
  character(80) cstring

  data meso_unit / 20 /
!**************************************************************************
 ncount=0

 inquire(file=trim(clistname),exist=fexist)
 if(fexist) then
    open (meso_unit,file=trim(clistname),form='formatted')
    do m=1,3
       read(meso_unit,*,end=131) cstring
    enddo
    n=0
    do 
       n=n+1
       read(meso_unit,*,end=131) clist(n)
    end do
131 continue
    ncount=n-1

    close(meso_unit)
 endif
end subroutine readin_rjlists


subroutine destroy_rjlists
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    destroy_rjlists
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-10-01  lueken - added subprogram doc block
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
  implicit none

  deallocate(cprovider)
  deallocate(w_rjlist)
  deallocate(t_rjlist)
  deallocate(p_rjlist)
  deallocate(q_rjlist)
  deallocate(csta_winduse)
  deallocate(t_day_rjlist)
  deallocate(t_night_rjlist)
  deallocate(q_day_rjlist)
  deallocate(q_night_rjlist)
  deallocate(td_rjlist)
  deallocate(mxtm_rjlist)
  deallocate(mitm_rjlist)
  deallocate(pmsl_rjlist)
  deallocate(howv_rjlist)
  deallocate(lcbas_rjlist)
  deallocate(tcamt_rjlist)
  deallocate(cldch_rjlist)
  if(wbinlistexist) deallocate(nwbaccpts)
  if(wbinlistexist) deallocate(csta_windbin)
  deallocate(csta_visuse)
end subroutine destroy_rjlists


subroutine get_sunangle(idate,dtime,dlon,dlat,sunangle) 

  use constants, only: deg2rad
  implicit none

  integer(i_kind),intent(in)::idate
  real(r_kind),intent(in)::dtime,dlon,dlat
  real(r_single),intent(out)::sunangle

  real(r_single),parameter:: s24 = 24._r_single
  real(r_single),parameter:: s180 = 180._r_single
  real(r_single),parameter:: s360 = 360._r_single

  integer(i_kind) iyear,imonth,iday,ihour
  integer(i_kind) iw3jdn
  integer(i_kind) jday2,iyear3,imonth3,iday3,idaywk3,idayyr3
  integer(i_kind) ndays
  real(r_single) rlat,rlon
  real(r_single) time
  character(10) date

  write(date,'( i10)') idate
  read (date,'(i4,3i2)') iyear,imonth,iday,ihour
  jday2=iw3jdn(iyear,imonth,iday)

  time=real(ihour,kind=r_single)+real(dtime,kind=r_single)

  if (time.gt.s24) then
     time = time-s24
     jday2=jday2+1
   elseif (time.lt.0) then
     time=time+s24
     jday2=jday2-1
  endif

  call w3fs26(jday2,iyear3,imonth3,iday3,idaywk3,idayyr3)
  !account for leap year
  if(mod(iyear3,4).eq.0) then
     ndays=366
  else
     ndays=365
  endif

  rlon=real(dlon/deg2rad,kind=r_single) ; if (rlon > s180) rlon=rlon-s360
  rlat=real(dlat/deg2rad,kind=r_single)

  call calcsun(idayyr3,ndays,time,rlat,rlon,sunangle)

end subroutine get_sunangle

subroutine get_wbinid(udbl,vdbl,nbins,ibin)

  use constants, only: zero, tiny_r_kind

  implicit none

  integer(i_kind),intent(in   ):: nbins
  real(r_double) ,intent(in   ):: udbl,vdbl
  integer(i_kind),intent(out  ):: ibin

! Declare local variables
  integer(i_kind) n
  real(r_kind) binwidth
  real(r_kind) ue,ve,wdir

  real(r_kind),parameter::r360=360._r_kind

  ue=real(udbl,kind=r_kind)
  ve=real(vdbl,kind=r_kind)

  binwidth=r360/real(max(1,nbins),kind=r_kind)

  call getwdir(ue,ve,wdir)

  if (abs(wdir)<=tiny_r_kind .or. abs(wdir-r360)<=tiny_r_kind) then  
     if (abs(ue)<=tiny_r_kind .and. abs(ve)<=tiny_r_kind) then
        ibin=nbins+1 !don't use if wind ob is calm
     else
        ibin=nbins
     endif
   else
     do n=1,nbins
        if ( wdir >= float(n-1)*binwidth .and. wdir < float(n)*binwidth ) then 
            ibin=n
            exit
         endif
      enddo
   endif

end subroutine get_wbinid

end module sfcobsqc

subroutine calcsun (idayr,idaysy,baltm,xlon,ylat,solar)

  implicit none

  integer(4),intent(in)::idayr,idaysy
  real(4),intent(in)::baltm,xlon,ylat
  real(4),intent(out)::solar
  real(4)::dangl,sdgl,cdgl,tsnoon,afy,safy,cafy,ssod,alon,sha,rlat

  dangl = 6.2831853 * (real(idayr) - 79.)/real(idaysy)
  sdgl = sin(dangl)
  cdgl = cos(dangl)
  !solar noon in utc
  tsnoon = -.030*sdgl-.120*cdgl+.330*sdgl*cdgl+.0016*sdgl**2-.0008
  !afy = angular fraction of year
  afy=6.2831853*(real(idayr)-1. + (baltm/24.))/real(idaysy)
  safy=sin(afy)
  cafy=cos(afy)
  !ssod = sin of dolar declination
  ssod = 0.3978492*sin(4.88578+afy+(.033420*safy)-(0.001388*cafy)+(.000696*safy*cafy)+(.000028*(safy**2-cafy**2)))
  alon=mod(720.+360.-xlon,360.)
  !sha = solar hour angle
  sha=0.0174532925*((15.*(tsnoon+baltm+36.))-alon)
  rlat=0.0174532925*ylat
  solar=57.29578*asin((ssod*sin(rlat))+sqrt(1.0-ssod**2)*cos(rlat)*cos(sha))
  return

end subroutine calcsun