subroutine read_no2(nread,ndata,nodata,jsatid,infile,gstime,lunout, &
           obstype,twind,sis,ithin,rmesh)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_no2                    read no2 data
!   prgmmr: yang             org: np23                date: 1998-05-15
!
! abstract:  This routine reads omi no2 observations.  Both layer
!            and total column values are read in.  The routine has
!            the ability to read both IEEE and BUFR format SBUV/2
!            no2 data files. OMI and GOME data is optionally thinned
!            to a specific resolution using simple quality control checks.
!
!            When running the gsi in regional mode, the code only
!            retains those observations that fall within the regional
!            domain
!
! program history log:
!   1998-05-15  yang, weiyu
!   1999-08-24  derber, j., treadon, r., yang, w., first frozen mpp version
!   2004-06-16  treadon - update documentation
!   2004-07-29  treadon - add only to module use, add intent in/out
!   2004-09-17  todling - fixed intent of jsatid
!   2004-12-02  todling - compilation in OSF1 forces big_endian for bufr files;
!                         need to force little_endian for ieee files
!   2004-12-22  kokron  - change cpp tokens to add support for ifort compiler
!                         efc does not have a convert option so it should use
!                         the other 'open'
!   2005-03-14  treadon - define numeric constants to r_kind precision
!   2005-05-12  wu - add OMI total ozone 
!   2005-06-27  guo     - bug fix: hour read from header was incorrect
!   2005-09-08  derber - modify to use input group time window
!   2005-09-19  treadon - add check on NOAA-17 sbuv data (toss bad data)
!   2005-10-17  treadon - add grid and earth relative obs location to output file
!   2005-10-18  treadon - remove array obs_load and call to sumload
!   2005-12-23  treadon - bound longitude to be less than 360.0
!   2006-01-26  treadon - remove ieee sbuv option
!   2006-02-03  derber  - modify for new obs control and obs count
!   2007-03-01  tremolet - measure time from beginning of assimilation window
!   2007-07-10  zhou    - modify to read version 8 SBUV/2 BUFR data(keep 
!                         option to read version 6 data), also add 
!                         total ozone and ozone profile quality control.
!   2007-09-11  h.liu - add kidsat for nimbus-7, n09, n11, n14
!   2007-10-16  zhou    - organize ozone flag control for all satellites
!   2008-04-16  h.liu   - thin OMI and read in GOME data
!   2008-05-27  safford - rm unused vars and uses
!   2008-05-30  treadon - accept  version8 poq=7 obs for further processing
!   2008-06-01  treadon - adjust logic to correctly handle zero length BUFR files
!   2008-06-03  treadon - add use_poq7 flag
!   2008-09-08  lueken  - merged ed's changes into q1fy09 code
!   2009-01-20  sienkiewicz - merge in changes for MLS ozone
!   2009-04-21  derber  - add ithin to call to makegrids
!   2009-3-05   h.liu   - read in OMI bufr, QC GOME2 and OMI
!   2009-7-02   h.liu   - toss the OMI data with AFBO=3 (c-pair correction) and clean up codes
!   2010-05-26  treadon - add timedif=zero for l4dvar (used in thinning)
!   2010-06-02  sienkiewicz - care for closing bufr other than for o3lev
!   2011-07-04  todling  - fixes to run either single or double precision
!   2011-08-01  lueken  - replaced F90 with f90 (no machine logic)
!   2012-10-12  h.liu  - read in MLS v2 Near Real Time (NRT) and v2.2 standard bufr data
!   2013-01-17  h.liu  - read in MLS v3 Near Real Time (NRT) 
!   2013-01-26  parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS)
!
!   input argument list:
!     obstype  - observation type to process
!     jsatid   - satellite id to read
!     infile   - unit from which to read no2 data
!     gstime   - analysis time in minutes from reference date
!     lunout   - unit to which to write data for further processing
!     obstype  - observation type to process
!     twind    - input group time window (hours)
!     sis      - satellite/instrument/sensor indicator
!     ithin    - flag to thin data
!     rmesh    - thinning mesh size (km)
!
!   output argument list:
!     nread    - number of sbuv/omi no2 observations read
!     ndata    - number of sbuv/omi no2 profiles retained for further processing
!     nodata   - number of sbuv/omi no2 observations retained for further processing
!
! remarks:
!   NCEP stopped producing IEEE format sbuv ozone files in April 2004.  
!   Hence, the IEEE portion of this routine no future application.  It 
!   is retained in the GSI package for use with retrospective runs.  The
!   IEEE portion of this routine may be removed from the GSI at a later date.
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,r_double,i_kind
  use satthin, only: makegrids,map2tgrid,destroygrids, &
      finalcheck,itxmax
  use gridmod, only: nlat,nlon,regional,tll2xy,rlats,rlons
  use constants, only: deg2rad,zero,rad2deg,one_tenth,r60inv,two
  use obsmod, only: iadate,&
      cloudfilter !JDE
  use convinfo, only: nconvtype, &
      icuse,ictype,ioctype
  use gsi_4dvar, only: l4dvar,iwinbgn,winlen
  use qcmod, only: use_poq7
  use mpimod, only : mype
  use omimod, only : nlevscatwt,presscatwt,no2pe
  implicit none

! Declare passed variables
  character(len=*),intent(in   ) :: obstype,infile,jsatid
  character(len=*),intent(in   ) :: sis
  integer(i_kind) ,intent(in   ) :: lunout,ithin
  integer(i_kind) ,intent(inout) :: nread
  integer(i_kind) ,intent(inout) :: ndata,nodata
  real(r_kind)    ,intent(in   ) :: gstime,twind,rmesh

! Declare local parameters
  real(r_kind),parameter:: r6   = 6.0_r_kind
  real(r_kind),parameter:: r76  = 76.0_r_kind
  real(r_kind),parameter:: r84  = 84.0_r_kind

  real(r_kind),parameter:: r360 = 360.0_r_kind
  real(r_kind),parameter:: rmiss = -9999.9_r_kind
!  real(r_kind),parameter:: badno2 = 10000.0_r_kind
  real(r_kind),parameter:: badno2 = 10000.0_r_kind

! Declare local variables
  logical outside,version6,version8,iuse
  
  character(2) version
  character(8) subset
  character(49) no2str
  character(42) no2str2
  character(23) tropomino2str2
  character(40) no2str3
  character(30) tropomino2str3

  integer(i_kind) maxobs,no2dat,nlno2
  integer(i_kind) idate,jdate,ksatid,kk,iy,iret,im,ihh,idd,lunin
  integer(i_kind) nmind,i
  integer(i_kind) imin,jskip(4)
  integer(i_kind) nmrecs,k,ilat,ilon,nreal,nchanl
! integer(i_kind) ithin,kidsat
  integer(i_kind) kidsat
  integer(i_kind) idate5(5)
  integer(i_kind) JULIAN,IDAYYR,IDAYWK
  integer(i_kind) itype, ikx
  integer(i_kind) isnd, ilev, iflg, mflg
  integer(i_kind) decimal,binary(14)


  integer(i_kind) itx,itt,ipoq7

  real(r_kind) tdiff,sstime,slons,slats,dlon,dlat,t4dv,timedif,crit1,dist1
  real(r_kind) slons0,slats0,rsat,solzen,solzenp,dlat_earth,dlon_earth
  real(r_kind) rsec, ppmv, prec, pres, pob, obserr, usage

! maximum number of observations set to 
  real(r_kind),allocatable,dimension(:,:):: no2out
  real(r_double) toq,poq
!  real(r_double),dimension(10):: hdrno2
  real(r_double),dimension(10):: hdrno2o
  real(r_double),dimension(4) :: hdrno22
  real(r_double),dimension(5) :: hdrno23
!  integer(i_kind),parameter :: nlevscatwt=35
  real(r_double) :: scatwt(nlevscatwt)
!  real(r_double) :: presscatwt(nlevscatwt),scatwt(nlevscatwt)
  real(r_kind),allocatable,dimension(:):: usage1
  integer(i_kind),allocatable,dimension(:):: ipos

!  real(r_double) totno2,totno2std
  real(r_double) slantno2,slantno2std
  real(r_double) no2trop,no2tropstd,amftrop,amftropstd,troppres,amftot
  real(r_double) trplay
  real(r_double) testread
  integer(i_kind) :: k0
  integer(i_kind) :: no2fill
  integer(i_kind) :: iskip,iskip2,iskipvcqf,itryread
  integer, parameter :: nlevavgkernel=34 ! for tropomi
  integer (i_kind) :: ipreskern,ikern,ino2trop
  real(r_double) :: avgkern(nlevavgkernel),pressavgkern(nlevavgkernel)
  logical :: first
  logical :: mask, tmask
  data first/.true./
  save first


  data no2str &
       / 'SAID CLAT CLON YEAR MNTH DAYS HOUR MINU SECO SOZA' /
! since 2009020412, the omi bufr contains fovn
  data no2str2 &
       / 'CLDMNT VZAN VCDQF XTRKQF' /
  data no2str3 &
       / 'NO2TRP NO2TRPSD AMFTRP AMFTRPSD TROPPRES'/
  data tropomino2str2 &
       / 'CLDMNT VZAN QUAL' /


  data lunin / 10 /

!**************************************************************************
! Set constants.  Initialize variables
  write(6,*)mype,'obstype',obstype,'readno2'
  write(6,*) 'JDE read_no2.f90 badno2: ',badno2
  call flush(6)
  jskip=0
  itryread=0
  iskip=0
  iskip2=0
  iskipvcqf=0
  rsat=999._r_kind
  no2fill=255
  maxobs=1e6
  ilon=3
  ilat=4
  ipoq7=0
  if (use_poq7) ipoq7=7


! Separately process  omi no2
  write(6,*)mype,'obstype read_no2',obstype,' infile ',infile
  call flush(6) 
  if ( obstype == 'tropomino2') then

     ! JDE Open file for writing excluded lat/lons
     open(unit=15, file='lat_exclude', access='SEQUENTIAL', form='UNFORMATTED', status='UNKNOWN') !JDE write out data
     open(unit=16, file='lon_exclude', access='SEQUENTIAL', form='UNFORMATTED', status='UNKNOWN') !JDE write out data
     open(unit=17, file='lat_all', access='SEQUENTIAL', form='UNFORMATTED', status='UNKNOWN') !JDE write out data
     open(unit=18, file='lon_all', access='SEQUENTIAL', form='UNFORMATTED', status='UNKNOWN') !JDE write out data
     open(unit=19, file='ll_mask', access='SEQUENTIAL', form='UNFORMATTED', status='UNKNOWN') !JDE write out data
     open(unit=20, file='qa_value', access='SEQUENTIAL', form='UNFORMATTED', status='UNKNOWN') !JDE write out data
     open(unit=21, file='time_mask', access='SEQUENTIAL', form='UNFORMATTED', status='UNKNOWN') !JDE write out data
     write(6,*)'JDE read_no2.f90 Open files for writing excluded lat/lon'
     write(6,*)'JDE read_no2.f90 r_kind is = ',r_kind

     if( mype == 1 ) write(6,*)'JDE obstype test passed. ',obstype
     call flush(6)
     ipreskern=17
     ikern=ipreskern+nlevavgkernel
     ino2trop=ikern+nlevavgkernel
!    Make thinning grids
     call makegrids(rmesh,ithin)
     write(6,*)'tropomino2 rmesh',rmesh,'ithin',ithin
     call flush(6)

     nmrecs=0
     write(6,*)'infile',infile
     open(lunin,file=infile,form='unformatted')
     call openbf(lunin,'IN',lunin)
     call datelen(10)
     call readmg(lunin,subset,idate,iret)
     write(6,*)'idate',idate
     if (subset == 'NC008013') then
        write(6,*)'READ_NO2:  TROPOMI data type, subset=',subset
     else
        write(6,*)'READ_NO2 :  *** WARNING: unknown no2 data type, subset=',subset
        write(6,*)' infile=',infile, ', lunin=',lunin, ', obstype=',obstype,', jsatid=',jsatid
        goto 170
     endif

!    Set dependent variables and allocate arraysn
!     nreal=14
!     nreal=18
!     nreal=17
!     nreal=17+nlevscatwt
     nreal=ino2trop-1
!     nreal=16
     nlno2=0
     nchanl=1
     no2dat=nreal+nchanl
     if( mype == 1 ) write(6,*)'JDE no2out before allocate: ',no2dat, itxmax
     allocate (no2out(no2dat,itxmax))
     do k=1,itxmax
        do i=1,no2dat
           no2out(i,k)=rmiss
        end do
     end do

     iy=0
     im=0
     idd=0
     ihh=0
     if(iret/=0) goto 160

1130 continue
     call readsb(lunin,iret)
!     write(6,*)mype,'130 loop iret',iret,'lunin',lunin
     if (iret/=0) then ! if no more subsets in current message then
        call readmg(lunin,subset,jdate,iret)
        if (iret/=0) goto 150 ! if no more messages in BUFR goto 150
        goto 1130
     endif

!    extract header information
     itryread=itryread+1
     call ufbint(lunin,hdrno22,3,1,iret,tropomino2str2)
     !if (hdrno22(3) < 100000000000.000) write(6,*)mype,'JDE iret tropomino2str2 ',iret
     !if (hdrno22(3) < 100000000000.000) write(6,*)mype,'hdrno22',hdrno22
     call ufbint(lunin,hdrno2o,10,1,iret,no2str)
     !if (hdrno22(3) < 100000000000.000) write(6,*)mype,'JDE iret no2str ',iret
     call ufbint(lunin,testread,1,1,iret,'VCDQF')
     !if (hdrno22(3) < 100000000000.000) write(6,*)mype,'VCDQF ',testread,iret
     call ufbint(lunin,testread,1,1,iret,'XTRKQF')
     !if (hdrno22(3) < 100000000000.000) write(6,*)mype,'XTRKQF ',testread,iret
     call ufbint(lunin,testread,1,1,iret,'CLDMNT')
     !if (hdrno22(3) < 100000000000.000) write(6,*)mype,'CLDMNT ',testread,iret
     call ufbint(lunin,testread,1,1,iret,'NO2')
     !if (hdrno22(3) < 100000000000.000) write(6,*)mype,'NO2 ',testread,iret
     rsat = hdrno2o(1); ksatid=rsat
     !if (hdrno22(3) < 100000000000.000) write(6,*)'jsatid ',jsatid,'rsat',rsat
     if(jsatid == 'aura'.or.jsatid == 'TROPOMI')kidsat = 785
     !if (hdrno22(3) < 100000000000.000) write(6,*)mype,'JDE ksatid, kidsat ',ksatid,kidsat
     if (ksatid /= kidsat) go to 1130
     flush(6)


     nmrecs=nmrecs+nlno2+1

!    Convert observation location to radians
     slats0= hdrno2o(2)
     slons0= hdrno2o(3)
!     write(6,*)mype,'slat',slats0,slons0
!     call flush(6)
     if(abs(slats0)>90._r_kind .or. abs(slons0)>r360) go to 1130  
     if(slons0< zero) slons0=slons0+r360
     if(slons0==r360) slons0=zero
     dlat_earth = slats0 * deg2rad
     dlon_earth = slons0 * deg2rad

     if(regional)then
        call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
        write(17) slats0
        write(18) slons0
        write(20) hdrno22(3)
        if(outside)then !JDE block
            mask=.TRUE.
        else
            mask=.FALSE.
        endif
        write(19)mask ! TRUE -- this is masked
        if(outside)jskip(1)=jskip(1)+1
        if(outside) write(15) slats0! JDE write outside lats/lons to file
        !write(6,*),'JDE read_no2.f90 slats0 = ',slats0
        if(outside) write(16) slons0! JDE write outside lats/lons to file
        if(outside) go to 1130
        !write(6,*)'JDE latitude = ',slats0
     else
        dlat = dlat_earth
        dlon = dlon_earth
        call grdcrd1(dlat,rlats,nlat,1)
        call grdcrd1(dlon,rlons,nlon,1)
     endif

! convert observation time to relative time
     idate5(1) = hdrno2o(4)  !year
     idate5(2) = hdrno2o(5)  !month
     idate5(3) = hdrno2o(6)  !day
     idate5(4) = hdrno2o(7)  !hour
     idate5(5) = hdrno2o(8)  !minute
     call w3fs21(idate5,nmind)
!     write(6,*)mype,'nmind',nmind,'idate5',idate5

     t4dv=real((nmind-iwinbgn),r_kind)*r60inv
!     write(6,*)mype,'ajl t4dv',t4dv,'nmind',nmind,iwinbgn,'ldvar',l4dvar
!     call flush(6)
     if (l4dvar) then
        if (t4dv<zero .OR. t4dv>winlen) go to 1130
     else
        sstime=real(nmind,r_kind)
        tdiff=(sstime-gstime)*r60inv
!        write(6,*)mype,'ajl tdiff',tdiff,'twind',twind
        if(abs(tdiff) > twind) jskip(2)=jskip(2)+1
        tmask=.FALSE.
        if(abs(tdiff) > twind) tmask=.TRUE. !JDE block
        write(21)tmask ! TRUE -- this is masked
        if(abs(tdiff) > twind) go to 1130
!        write(6,*)'twind',twind,'tdiff',tdiff,'hr',hdrno2o(7)+hdrno2o(8)/60.,'lat',slats0,slons0
!        call flush(6)
     end if

!    extract total no2
!     call ufbint(lunin,totno2,1,1,iret,'NO2')
!     write(6,*)mype,'totno2',totno2
!     call flush(6)
!    also NO2 precision
!     call ufbint(lunin,totno2std,1,1,iret,'NO2STD')
!    extrack slant column no2
     call ufbint(lunin,slantno2,1,1,iret,'SLNO2')
     call ufbint(lunin,slantno2std,1,1,iret,'SLNO2STD')
     call ufbint(lunin,no2trop,1,1,iret,'NO2TRP')
     call ufbint(lunin,no2tropstd,1,1,iret,'NO2TRPSD')
     call ufbint(lunin,amftrop,1,1,iret,'AMFTRP')
     call ufbint(lunin,amftot,1,1,iret,'AMFTOT')
     call ufbint(lunin,trplay,1,1,iret,'TRPLAY')
!    check to see what badno2 would be
!     write(6,*)mype,'totno2',totno2,'basno2',badno2
!     call flush(6)
     if (no2trop > badno2 .or. no2trop < 0.0 ) jskip(3)=jskip(3)+1
     !if (no2trop > badno2 .or. no2trop < 0.0 ) write(6,*)'JDE read_no2.f90 jskip no2trop val: ',no2trop
     if (no2trop > badno2 .or. no2trop < 0.0 ) goto 1130
! JDE skip this block, it is redundant, check performed in BUFR converter
! JDE     if(hdrno22(3)<=.5)then
! JDE!      bad flag set
! JDE!       write(6,*)mype,'skip  **** hdrno22(3)',hdrno22(3)
! JDE       iskipvcqf=iskipvcqf+1
! JDE       write(6,*)'JDE read_no2.f90, bad QA val = ',hdrno22(3)
! JDE       go to 1130
! JDE     endif

!    ajl add cloud check
     
     ! JDE replace cloud check with "cloudfilter" car
     !JDE if (hdrno22(1) > 30. ) jskip(4)=jskip(4)+1
     !JDE if (hdrno22(1) > 30. ) write(6,*)'cloud ',hdrno22(1)
     !JDE if (hdrno22(1) > 30. ) go to 1130

     if (hdrno22(1) > cloudfilter ) jskip(4)=jskip(4)+1
     if (hdrno22(1) > cloudfilter ) write(6,*)'cloud ',hdrno22(1)
     if (hdrno22(1) > cloudfilter ) go to 1130

!     write(6,*)mype,'got though',hdrno22(1)
!       write(6,*) 'NO2TRP NO2TRPSD AMFTRP AMFTRPSD TROPPRES'
!     write(6,*)'hdrno23',hdrno23


!    thin OMI data

     if (l4dvar) then 
        timedif = zero 
     else 
        timedif = r6*abs(tdiff)        ! range:  0 to 18 
     endif 
     crit1 = 0.01_r_kind+timedif
     call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis)
     if(.not.iuse)then
           iskip2=iskip2+1
           write(6,*)'out of bounds? ',dlat_earth,dlon_earth,'timedif',timedif
     endif
     if(.not. iuse)go to 1130
 
     call finalcheck(dist1,crit1,itx,iuse)
     if(.not.iuse)then
          iskip2=iskip2+1
!        write(6,*)'no2 finalcheck',iskip2,dlat_earth,dlon_earth,'itx',itx,'crit1',crit1,'dist1',dist1
     endif
     if(.not. iuse)go to 1130
   
     ndata=ndata+1
     nodata=ndata
!     write(6,*)mype,'ajl ndata',ndata
!     call flush(6)



     no2out(1,itx)=rsat
     no2out(2,itx)=t4dv
     no2out(3,itx)=dlon               ! grid relative longitude
     no2out(4,itx)=dlat               ! grid relative latitude
     no2out(5,itx)=dlon_earth*rad2deg ! earth relative longitude (degrees)
     no2out(6,itx)=dlat_earth*rad2deg ! earth relative latitude (degrees)
!    this could be vcdqflags
     no2out(7,itx)=hdrno22(3)         ! vertical column no2  quality code
     no2out(8,itx)=hdrno2o(10)         ! solar zenith angle
     no2out(9,itx)=amftot            ! air mass factor total atm
!    this could be xtrackflgs
     no2out(10,itx)=hdrno22(1)        !  cloud amount
     no2out(11,itx)=hdrno22(2)        !  vzan
!     no2out(12,itx)=totno2
     no2out(12,itx)=no2tropstd
     no2out(13,itx)=slantno2
     no2out(14,itx)=slantno2std
     no2out(15,itx)=amftrop           ! air mass factor troposphere
     no2out(17,itx)=trplay
!     write(6,*)mype,'slat',slats0,slons0
!     write(6,*)'no2out',no2out(1:17,itx)
     call ufbrep(lunin,pressavgkern,1,nlevavgkernel,iret,'PRLNO2')
!     write(6,*)'ret prlno2 ',iret
     call ufbrep(lunin,avgkern,1,nlevavgkernel,iret,'AVKERN')
!     write(6,*)'iret avkern',iret,'itryread',itryread
!     call flush(6)
     do k=1,nlevavgkernel
       no2out(ipreskern+k,itx)=pressavgkern(k)
       no2out(ikern+k,itx)=avgkern(k)
!       write(6,*)'pres',k,pressavgkern(k),avgkern(k)
     end do
!     no2out(ikern,itx)=0.0 ! top pressure
     no2out(ino2trop,itx)=no2trop
!     write(6,*)'no2trop',no2trop
!     call flush(6)
!    last should be no2trop now
!     write(500+mype,*)hdrno22(3),'ratio error/no2trop',no2tropstd/no2trop

!     if(mype.eq.0.and.itx<10)then
!       write(6,*)itx,'rsat',rsat,'t4dv',t4dv,'dlon',no2out(5:6,itx)
!       write(6,*)'colno2',no2out(12,itx),'sza',no2out(8,itx)
!       write(6,*)'vdqf',no2out(7,itx),'xtrackqf',no2out(9,itx)
!       write(6,*)'std',no2out(13,itx),'cloud',no2out(10,itx)
!       write(6,*)'slant',no2out(14:15,itx)
!       do k=1,nlevscatwt
!         print *,'psct',presscatwt(k),'wt',no2out(17+k,itx)
!       enddo
!     endif

!    End of loop over observations
     go to 1130
  elseif ( obstype == 'omino2') then
     if( mype == 1 ) write(6,*)'JDE obstype: ',obstype

!    Make thinning grids
     call makegrids(rmesh,ithin)
     write(6,*)'omino2 rmesh',rmesh,'ithin',ithin
     call flush(6)

     nmrecs=0
     write(6,*)'infile',infile
     open(lunin,file=infile,form='unformatted')
     call openbf(lunin,'IN',lunin)
     call datelen(10)
     call readmg(lunin,subset,idate,iret)
     write(6,*)'idate',idate
     if (subset == 'NC008013') then
        write(6,*)'READ_NO2:  OMI data type, subset=',subset
     else
        write(6,*)'READ_NO2 :  *** WARNING: unknown no2 data type, subset=',subset
        write(6,*)' infile=',infile, ', lunin=',lunin, ', obstype=',obstype,', jsatid=',jsatid
        goto 170
     endif

!    Set dependent variables and allocate arraysn
!     nreal=14
!     nreal=18
!     nreal=17
     nreal=17+nlevscatwt
!     nreal=16
     nlno2=0
     nchanl=1
     no2dat=nreal+nchanl
     allocate (no2out(no2dat,itxmax))
     do k=1,itxmax
        do i=1,no2dat
           no2out(i,k)=rmiss
        end do
     end do

     iy=0
     im=0
     idd=0
     ihh=0
     if(iret/=0) goto 160

130  continue
     call readsb(lunin,iret)
!     write(6,*)mype,'130 loop iret',iret,'lunin',lunin
     if (iret/=0) then
        call readmg(lunin,subset,jdate,iret)
        if (iret/=0) goto 150
        goto 130
     endif

!    extract header information
     itryread=itryread+1
     call ufbint(lunin,hdrno2o,10,1,iret,no2str)
!     write(6,*)mype,'ajl iret no2str ',iret
     call ufbint(lunin,hdrno22,4,1,iret,no2str2)
     call ufbint(lunin,hdrno23,5,1,iret,no2str3)
!     write(6,*)mype,'ajl iret no2str2 ',iret
!     write(6,*)mype,'hdrno22',hdrno22
     call ufbint(lunin,testread,1,1,iret,'VCDQF')
!     write(6,*)mype,'VCDQF',testread,iret
!     call flush(6)
     call ufbint(lunin,testread,1,1,iret,'XTRKQF')
!     write(6,*)mype,'XTRKQF',testread,iret
!     call flush(6)
     call ufbint(lunin,testread,1,1,iret,'CLDMNT')
!     write(6,*)mype,'CLDMNT',testread,iret
!     call flush(6)
     call ufbint(lunin,testread,1,1,iret,'NO2')
!     write(6,*)mype,'NO2',testread,iret
     rsat = hdrno2o(1); ksatid=rsat

     if(jsatid == 'aura')kidsat = 785
!     write(6,*)mype,'ajl ksatid',ksatid,kidsat
     if (ksatid /= kidsat) go to 130


     nmrecs=nmrecs+nlno2+1

!    Convert observation location to radians
     slats0= hdrno2o(2)
     slons0= hdrno2o(3)
!     write(6,*)mype,'slat',slats0,slons0
!     call flush(6)
     if(abs(slats0)>90._r_kind .or. abs(slons0)>r360) go to 130  
     if(slons0< zero) slons0=slons0+r360
     if(slons0==r360) slons0=zero
     dlat_earth = slats0 * deg2rad
     dlon_earth = slons0 * deg2rad

     if(regional)then
        call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
        if(outside)jskip(1)=jskip(1)+1
        if(outside) go to 130
     else
        dlat = dlat_earth
        dlon = dlon_earth
        call grdcrd1(dlat,rlats,nlat,1)
        call grdcrd1(dlon,rlons,nlon,1)
     endif

! convert observation time to relative time
     idate5(1) = hdrno2o(4)  !year
     idate5(2) = hdrno2o(5)  !month
     idate5(3) = hdrno2o(6)  !day
     idate5(4) = hdrno2o(7)  !hour
     idate5(5) = hdrno2o(8)  !minute
     call w3fs21(idate5,nmind)
!     write(6,*)mype,'nmind',nmind,'idate5',idate5

     t4dv=real((nmind-iwinbgn),r_kind)*r60inv
!     write(6,*)mype,'ajl t4dv',t4dv,'nmind',nmind,iwinbgn,'ldvar',l4dvar
!     call flush(6)
     if (l4dvar) then
        if (t4dv<zero .OR. t4dv>winlen) go to 130
     else
        sstime=real(nmind,r_kind)
        tdiff=(sstime-gstime)*r60inv
!        write(6,*)mype,'ajl tdiff',tdiff,'twind',twind
        if(abs(tdiff) > twind) jskip(2)=jskip(2)+1
        tmask=.FALSE.
        if(abs(tdiff) > twind) tmask=.TRUE. !JDE block
        write(21)tmask ! TRUE -- this is masked
        if(abs(tdiff) > twind) go to 130
!        write(6,*)'twind',twind,'tdiff',tdiff,'hr',hdrno2o(7)+hdrno2o(8)/60.,'lat',slats0,slons0
!        call flush(6)
     end if

!    extract total no2
!     call ufbint(lunin,totno2,1,1,iret,'NO2')
!     write(6,*)mype,'totno2',totno2
!     call flush(6)
!    also NO2 precision
!     call ufbint(lunin,totno2std,1,1,iret,'NO2STD')
!    extrack slant column no2
     call ufbint(lunin,slantno2,1,1,iret,'SLNO2')
     call ufbint(lunin,slantno2std,1,1,iret,'SLNO2STD')
     call ufbint(lunin,no2trop,1,1,iret,'NO2TRP')
     call ufbint(lunin,no2tropstd,1,1,iret,'NO2TRPSD')
     call ufbint(lunin,amftrop,1,1,iret,'AMFTRP')
     call ufbint(lunin,amftropstd,1,1,iret,'AMFTRPSD')
     call ufbint(lunin,troppres,1,1,iret,'TROPPRES')
!    check to see what badno2 would be
!     write(6,*)mype,'totno2',totno2,'basno2',badno2
!     call flush(6)
     if (no2trop > badno2 .or. no2trop < 0.0 ) jskip(3)=jskip(3)+1
     !if (no2trop > badno2 .or. no2trop < 0.0 ) write(6,*)'JDE read_no2.f90 jskip no2trop val: ',no2trop
     if (no2trop > badno2 .or. no2trop < 0.0 ) goto 130
!    XTRACK gives row anonm?
     if (nint(hdrno22(4)).ne.0.and.nint(hdrno22(4)).ne.no2fill)then
!       write(6,*)'skip **** hdrno22(4)',hdrno22(4),'nofill',no2fill
!       call flush(6)
       iskip=iskip+1
       go to 130
     endif
!    Bit 10 in TOQF represents row anomaly. 
!     decimal=int(hdrno22(6))
!     call dec2bin(decimal,binary,14)
!     if (binary(10) == 1 ) then
!        goto 130
!     endif
     if(mod(nint(hdrno22(3)),2).eq.1)then
!      bad flag set
!       write(6,*)mype,'skip  **** hdrno22(3)',hdrno22(3)
       iskipvcqf=iskipvcqf+1
       go to 130
     endif

!    ajl add cloud check
     
     ! JDE replace cloud check with "cloudfilter" var
     !JDE if (hdrno22(1) > 30. ) jskip(4)=jskip(4)+1
     !JDE if (hdrno22(1) > 30. ) go to 130

     if (hdrno22(1) > cloudfilter ) jskip(4)=jskip(4)+1
     if (hdrno22(1) > cloudfilter ) go to 130

!     write(6,*)mype,'got though',hdrno22(1)
!       write(6,*) 'NO2TRP NO2TRPSD AMFTRP AMFTRPSD TROPPRES'
!     write(6,*)'hdrno23',hdrno23


!    thin OMI data

     if (l4dvar) then 
        timedif = zero 
     else 
        timedif = r6*abs(tdiff)        ! range:  0 to 18 
     endif 
     crit1 = 0.01_r_kind+timedif
     call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis)
     if(.not.iuse)then
           iskip2=iskip2+1
           write(6,*)'out of bounds? ',dlat_earth,dlon_earth,'timedif',timedif
     endif
     if(.not. iuse)go to 130
 
     call finalcheck(dist1,crit1,itx,iuse)
     if(.not.iuse)then
          iskip2=iskip2+1
!        write(6,*)'no2 finalcheck',iskip2,dlat_earth,dlon_earth,'itx',itx,'crit1',crit1,'dist1',dist1
     endif
     if(.not. iuse)go to 130
   
     ndata=ndata+1
     nodata=ndata
!     write(6,*)mype,'ajl ndata',ndata
!     call flush(6)



     no2out(1,itx)=rsat
     no2out(2,itx)=t4dv
     no2out(3,itx)=dlon               ! grid relative longitude
     no2out(4,itx)=dlat               ! grid relative latitude
     no2out(5,itx)=dlon_earth*rad2deg ! earth relative longitude (degrees)
     no2out(6,itx)=dlat_earth*rad2deg ! earth relative latitude (degrees)
!    this could be vcdqflags
     no2out(7,itx)=hdrno22(3)         ! vertical column no2  quality code
     no2out(8,itx)=hdrno2o(10)         ! solar zenith angle
     no2out(9,itx)=hdrno22(4)         ! xtrack qf contains row anomaly flag
!    this could be xtrackflgs
     no2out(10,itx)=hdrno22(1)        !  cloud amount
     no2out(11,itx)=hdrno22(4)        !  vzan
!     no2out(12,itx)=totno2
     no2out(12,itx)=no2tropstd
     no2out(13,itx)=slantno2
     no2out(14,itx)=slantno2std
     no2out(15,itx)=amftrop
     no2out(16,itx)=amftropstd
     no2out(17,itx)=troppres
!     write(6,*)'nlevscatwt',nlevscatwt,shape(presscatwt),shape(scatwt)
     if(first)then
       call ufbrep(lunin,presscatwt,1,nlevscatwt,iret,'PRLNO2')
!       write(6,*)'iret prlcno2',iret,'itx',itx
       first=.false.
!       write(6,*)'r_kind',r_kind,'r_double',r_double
!       do k=1,nlevscatwt
!         write(6,*)'presscatwt',presscatwt(k)
!         if(presscatwt(k)<0.0.or.presscatwt(k)>2000.)then
!           write(6,*)'presscatwt',k,presscatwt(k)
!           first=.true.
!         endif
!       end do
       if(.not.first)then
!          print *,'did presscatwt',presscatwt,'mype',mype
         no2pe=mype
       endif
     endif
     call ufbrep(lunin,scatwt,1,nlevscatwt,iret,'SCATWT')
!     write(6,*)'iret scatwt',iret
     do k=1,nlevscatwt
       no2out(17+k,itx)=scatwt(k)
     end do
     no2out(17+nlevscatwt+1,itx)=no2trop
!    last should be no2trop now
!     write(500+mype,*)hdrno22(3),'ratio error/no2trop',no2tropstd/no2trop

!     if(mype.eq.0.and.itx<10)then
!       write(6,*)itx,'rsat',rsat,'t4dv',t4dv,'dlon',no2out(5:6,itx)
!       write(6,*)'colno2',no2out(12,itx),'sza',no2out(8,itx)
!       write(6,*)'vdqf',no2out(7,itx),'xtrackqf',no2out(9,itx)
!       write(6,*)'std',no2out(13,itx),'cloud',no2out(10,itx)
!       write(6,*)'slant',no2out(14:15,itx)
!       do k=1,nlevscatwt
!         print *,'psct',presscatwt(k),'wt',no2out(17+k,itx)
!       enddo
!     endif

!    End of loop over observations
     go to 130

! End of OMI block
  endif


! Jump here when eof detected
150 continue
  write(6,*)mype,'did read_no2 ajl ****** itxmax ',itxmax,' no2dat',no2dat

! If gome or omi data, compress no2out array to thinned data
  if (obstype=='omino2'.or.obstype=='tropomino2') then
     if( mype == 1 ) write(6,*)'JDE no2out thinning if called'
     kk=0
     do k=1,itxmax
        if (no2out(1,k)>zero) then
           kk=kk+1
           do i=1,no2dat
              no2out(i,kk)=no2out(i,k)
           end do
        endif
     end do
     ndata=kk
     nodata=ndata
  endif
  write(6,*)'omino2 ndata',ndata,' iskip2 ',iskip2

! Write header record and data to output file for further processing
  write(6,*)'itryread no2 ',itryread
  write(6,*)mype,'output no2 lunout ',lunout,'nreal ',nreal,nchanl,'ilat,ilon ',ilat,ilon,'ndata ',ndata
  call flush(6)
  if( mype == 1 ) write(6,*)'JDE write to obs_intpu.0001 start'
  write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
  call flush(lunout)
  write(lunout) ((no2out(k,i),k=1,no2dat),i=1,ndata)
  if( mype == 1 ) write(6,*)'JDE write to obs_intpu.0001 end'
  call flush(lunout)
  nread=nmrecs

! Deallocate local arrays
160 continue
  deallocate(no2out)

! Close unit to input data file
170 continue
  if (obstype .ne. 'o3lev') call closbf(lunin)
  close(lunin)

! Deallocate satthin arrays
  if (obstype == 'omino2' .or. obstype == 'gome'.or.obstype=='tropomino2')call destroygrids
  if (iskip+iskip2.ne.0)write(6,*)'iskip no2',iskip,iskip2,' ndata ',ndata
  do k=1,4
    if(jskip(k).ne.0)then
      write(6,*)k,' jskip ',jskip(k)
    endif
  end do

  ! JDE close lats/lons files
  close(unit=15)!status='SAVE') !JDE write out data
  close(unit=16)!status='SAVE') !JDE write out data
  close(unit=17)!status='SAVE') !JDE write out data
  close(unit=18)!status='SAVE') !JDE write out data
  close(unit=19)!status='SAVE') !JDE write out data
  close(unit=20)!status='SAVE') !JDE write out data
  close(unit=21)!status='SAVE') !JDE write out data

  write(6,*)'JDE read_no2.f90 itryread  = ',itryread
  write(6,*)'JDE read_no2.f90 iskip     = ',iskip
  write(6,*)'JDE read_no2.f90 jskip     = ',jskip
  write(6,*)'JDE read_no2.f90 iskip2    = ',iskip2
  write(6,*)'JDE read_no2.f90 iskipvcqf = ',iskipvcqf
  write(6,*)'-----------------JDE END OF read_no2.f90 ----------------------'
  write(6,*)'---------------------------------------------------------------'

  return
  
end subroutine read_no2
