subroutine read_ahi(mype,val_img,ithin,rmesh,jsatid,gstime,&
     infile,lunout,obstype,nread,ndata,nodata,twind,sis, &
     mype_root,mype_sub,npe_sub,mpi_comm_sub,nobs,nrec_start,dval_use)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_ahi                    read himawari-8 ahi data
!   prgmmr: zaizhong ma     org: np23                date: 2014-11-19
!
! abstract:  This routine reads Himawari8 AHI radiance (brightness
!            temperature) files.  Optionally, the data are thinned to
!            a specified 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:
!   2014-11-19 zaizhong began developping with the proxy data from JMA
!   2014-12-12 zaizhong done the first version of this subroutine
!   2014-12-23 zaizhong cleaned up and finalized with the proxy data
!   2015-03-23 zaizhong cleaned up and finalized with the real sample data
!   2015-09-17 Thomas   add l4densvar and thin4d to data selection procedure
!   2016-03-11 j. guo   Fixed {dlat,dlon}_earth_deg in the obs data stream
!   2018-05-21 j.jin    Added time-thinning. Moved the checking of thin4d into satthin.F90.
!
!   input argument list:
!     mype     - mpi task id
!     val_img  - weighting factor applied to super obs
!     ithin    - flag to thin data
!     rmesh    - thinning mesh size (km)
!     jsatid   - satellite to read
!     gstime   - analysis time in minutes from reference date
!     infile   - unit from which to read BUFR data
!     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
!     nrec_start - first subset with useful information
!
!   output argument list:
!     nread    - number of BUFR GOES imager observations read
!     ndata    - number of BUFR GOES imager profiles retained for further processing
!     nodata   - number of BUFR GOES imager observations retained for further processing
!     nobs     - array of observations on each subdomain for each processor
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,r_double,i_kind
  use satthin, only: super_val,itxmax,makegrids,map2tgrid,destroygrids, &
      checkob,finalcheck,score_crit
  use satthin, only: radthin_time_info,tdiff2crit
  use obsmod,  only: time_window_max
  use gridmod, only: diagnostic_reg,regional,nlat,nlon,txy2ll,tll2xy,rlats,rlons
  use constants, only: deg2rad,zero,one,rad2deg,r60inv,r60
  use obsmod, only: bmiss
  use radinfo, only: iuse_rad,jpch_rad,nusis
  use gsi_4dvar, only: l4dvar,iwinbgn,winlen,l4densvar
  use deter_sfc_mod, only: deter_sfc
  use gsi_nstcouplermod, only: nst_gsi,nstinfo
  use gsi_nstcouplermod, only: gsi_nstcoupler_skindepth, gsi_nstcoupler_deter
  use file_utility, only : get_lun     
  use mpimod, only: npe
  implicit none

! Declare passed variables
  character(len=*),intent(in   ) :: infile,obstype,jsatid
  character(len=*),intent(in  ) :: sis
  integer(i_kind) ,intent(in   ) :: mype,lunout,ithin,nrec_start
  integer(i_kind),dimension(npe)  ,intent(inout) :: nobs
  integer(i_kind) ,intent(inout) :: ndata,nodata
  integer(i_kind) ,intent(inout) :: nread
  real(r_kind)    ,intent(in   ) :: rmesh,gstime,twind
  real(r_kind)    ,intent(inout) :: val_img
  integer(i_kind) ,intent(in   ) :: mype_root
  integer(i_kind) ,intent(in   ) :: mype_sub
  integer(i_kind) ,intent(in   ) :: npe_sub
  integer(i_kind) ,intent(in   ) :: mpi_comm_sub
  logical         ,intent(in)    :: dval_use

! Declare local parameters
  integer(i_kind),parameter:: nimghdr=11
  real(r_kind),parameter:: r360=360.0_r_kind
  real(r_kind),parameter:: r180=180.0_r_kind
  real(r_kind),parameter:: tbmin=50.0_r_kind
  real(r_kind),parameter:: tbmax=550.0_r_kind
  character(80),parameter:: hdrh8  = &            ! Himawari-8 AHI header
              'SAID YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH SAZA SOZA'

! Declare local variables
  logical outside,iuse,assim

  character(8) subset

  integer(i_kind) nchanl,ilath,ilonh,ilzah,iszah,irec,next,maxinfo,nchn
  integer(i_kind) nmind,lnbufr,idate,ilat,ilon
  integer(i_kind) ireadmg,ireadsb,iret,nreal,nele,itt
  integer(i_kind) itx,i,k,isflg,kidsat,n,iscan,idomsfc
  integer(i_kind) idate5(5)
  integer(i_kind),allocatable,dimension(:)::nrec

  real(r_kind) dg2ew,sstime,tdiff,t4dv,sfcr
  real(r_kind) dlon,dlat,crit1,dist1
  real(r_kind) dlon_earth,dlat_earth
  real(r_kind) dlon_earth_deg,dlat_earth_deg
  real(r_kind) pred
  real(r_kind),dimension(0:3):: sfcpct
  real(r_kind),dimension(0:3):: ts
  real(r_kind) :: tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10
  real(r_kind) :: zob,tref,dtw,dtc,tz_tr
  real(r_kind),allocatable,dimension(:,:):: data_all

  logical :: allchnmiss
  real(r_kind) rclrsky,rcldfrc
  real(r_double),dimension(nimghdr) :: hdrh8arr       !   Himawari8 AHI data
  real(r_double),allocatable,dimension(:,:) :: dataahi,dataahibt,dataahisd  !   Himawari8 AHI data for NCLDMNT,BT,SDTB
  real(r_kind),dimension(0:4):: rlndsea

!--start-- variables for AHI cloud detection
  real(r_kind)                        :: ts_coef0
  real(r_kind), dimension(4)          :: ts_coef
  integer(i_kind), dimension(2)       :: ts_ichan
  real(r_kind)                        :: seca, dbt_ts
  real(r_kind)                        :: dts_thresh = 330.0_r_kind
  integer(i_kind)                     :: qc_thresh = 1
  real(r_kind), dimension(2)          :: bt_ts
  !---regression sst from split window test
  real(r_kind)                        :: ts_reg
  !---difference between sst from regression and surface
  real(r_kind)                        :: sst_test
!-end--- variables for AHI cloud detection

  real(r_kind) cdist,disterr,disterrmax,dlon00,dlat00
  integer(i_kind) ntest
  real(r_kind)    :: ptime,timeinflat,crit0
  integer(i_kind) :: ithin_time,n_tbin,it_mesh

!**************************************************************************
! Initialize variables
  lnbufr = 10
  disterrmax=zero
  ntest=0
  dg2ew = r360*deg2rad

  ilon=3
  ilat=4

  if (nst_gsi > 0 ) then
     call gsi_nstcoupler_skindepth(obstype, zob)         ! get penetration depth (zob) for the obstype
  endif

  rlndsea(0) = zero
  rlndsea(1) = 15._r_kind
  rlndsea(2) = 10._r_kind
  rlndsea(3) = 15._r_kind
  rlndsea(4) = 30._r_kind

  nread=0
  ndata=0
  nodata=0

  nchn=12         ! total numer of channels
  nchanl=10       ! total number of channels with valid BT and SDTB
  ilath=8        ! the position of latitude in the header
  ilonh=9        ! the position of longitude in the header
  ilzah=10       ! satellite zenith angle
  iszah=11       ! solar zenith angle

! If all channels of a given sensor are set to monitor or not
! assimilate mode (iuse_rad<1), reset relative weight to zero.
! We do not want such observations affecting the relative
! weighting between observations within a given thinning group.

  assim=.false.
  search: do i=1,jpch_rad
     if ((nusis(i)==sis) .and. (iuse_rad(i)>0)) then
        assim=.true.
        exit search
     endif
  end do search
  if (.not.assim) val_img=zero

  call radthin_time_info(obstype, jsatid, sis, ptime, ithin_time)
  if( ptime > 0.0_r_kind) then
     n_tbin=nint(2*time_window_max/ptime)
  else
     n_tbin=1
  endif
! Make thinning grids
  call makegrids(rmesh,ithin,n_tbin=n_tbin)


! Open bufr file.
  open(lnbufr,file=trim(infile),form='unformatted')
  call openbf(lnbufr,'IN',lnbufr)
  call datelen(10)
  call readmg(lnbufr,subset,idate,iret)

! Check the data set
  if( iret/=0) then
     write(6,*) 'READ_AHI: SKIP PROCESSING OF AHI FILE'
     write(6,*) 'infile=', lnbufr, infile
     return
  endif

  allocate(dataahi(1,nchn))      ! NCLDMNT: 2 for ASR, not channel dependent; ncld for CSR, chn dependent
  allocate(dataahibt(1,nchn))    ! BT: channel dependent: all, clear, cloudy, low, middle and high clouds
  allocate(dataahisd(1,nchn))    ! SDTB: channel dependent: all, clear, cloudy, low, middle and high clouds

! Allocate arrays to hold all data for given satellite
  maxinfo=32
  maxinfo=maxinfo+nchanl
  if(dval_use) maxinfo = maxinfo + 2
  nreal = maxinfo + nstinfo
  nele  = nreal   + nchanl
  allocate(data_all(nele,itxmax),nrec(itxmax))

  call closbf(lnbufr)
  close(lnbufr)
  open(lnbufr,file=trim(infile),form='unformatted')
  call openbf(lnbufr,'IN',lnbufr)
  if(jsatid == 'himawari8') then
     kidsat = 173
  elseif (jsatid == 'himawari9') then
     kidsat = 174
  else
     write(6,*) 'READ_AHI: Unrecognized value for jsatid '//jsatid//': RETURNING'
     return
  end if


  next=0
  nrec=999999
  irec=0

!---regression coefficients trained in clear simulation
  !   using CRTM v2.1.3 AHI and ECMWF
  ts_coef    = (/1.16778_r_kind, -1.27133_r_kind, 0.416716_r_kind, &
                 2.16380_r_kind/)
  ts_coef0   = -51.0104_r_kind
  !---should be channels 11.2 and 12.38 microns
  ts_ichan   = (/7,8/) 
  !---threshold for difference in regression from tsavg
  !   COAT used 1.25 so we may need to make this smaller.
  dts_thresh = 2.00
  !---currently not used
  qc_thresh  = -1

! Big loop over bufr file
  read_msg: do while(IREADMG(lnbufr,subset,idate) >= 0)
     irec=irec+1
     if(irec < nrec_start) cycle read_msg
     next=next+1
     if(next == npe_sub)next=0
     if(next /= mype_sub)cycle
     read_loop: do while (IREADSB(lnbufr) == 0)

!       Read through each reacord
        call ufbint(lnbufr,hdrh8arr,nimghdr,1,iret,hdrh8)
        if(hdrh8arr(1) /= kidsat) cycle read_loop

!       only the first 10 chns are valid
        call ufbrep(lnbufr,dataahibt,1,nchn,iret,'TMBRST')
        nread=nread+nchanl

        allchnmiss=.true.
        do n=1,nchanl
           if(dataahibt(1,n)<500.0_r_kind)  then
              allchnmiss=.false.
           end if
        end do
        if(allchnmiss) cycle read_loop

!      first step QC filter out data with SAZA>60degree,check use 60 or 65
        if (hdrh8arr(ilzah) >65.0_r_kind .or. hdrh8arr(iszah) >= r180) then
          print*, 'SAZA & Satellite azimuth',hdrh8arr(ilzah),hdrh8arr(iszah)
          cycle read_loop
        end if

!       Convert obs location from degrees to radians
        if (hdrh8arr(ilonh)>=r360) hdrh8arr(ilonh)=hdrh8arr(ilonh)-r360
        if (hdrh8arr(ilonh)< zero) hdrh8arr(ilonh)=hdrh8arr(ilonh)+r360

        dlon_earth_deg = hdrh8arr(ilonh)
        dlat_earth_deg = hdrh8arr(ilath)
        dlon_earth=hdrh8arr(ilonh)*deg2rad
        dlat_earth=hdrh8arr(ilath)*deg2rad

!       If regional, map obs lat,lon to rotated grid.
        if(regional)then

!          Convert to rotated coordinate.  dlon centered on 180 (pi), 
!          so always positive for limited area
           call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)

           if(diagnostic_reg) then
              call txy2ll(dlon,dlat,dlon00,dlat00)
              ntest=ntest+1
              cdist=sin(dlat_earth)*sin(dlat00)+cos(dlat_earth)*cos(dlat00)* &
                   (sin(dlon_earth)*sin(dlon00)+cos(dlon_earth)*cos(dlon00))
              cdist=max(-one,min(cdist,one))
              disterr=acos(cdist)*rad2deg
              disterrmax=max(disterrmax,disterr)
           end if

!          Check to see if in domain.  outside=.true. if dlon_earth,
!          dlat_earth outside domain, =.false. if inside
           if(outside) cycle read_loop

!       Global case
        else
           dlon=dlon_earth
           dlat=dlat_earth
           call grdcrd1(dlat,rlats,nlat,1)
           call grdcrd1(dlon,rlons,nlon,1)
        endif

!       Compare relative obs time with window.  If obs 
!       falls outside of window, don't use this obs
        idate5(1) = hdrh8arr(2)     !year
        idate5(2) = hdrh8arr(3)     ! month
        idate5(3) = hdrh8arr(4)     ! day
        idate5(4) = hdrh8arr(5)     ! hours
        idate5(5) = hdrh8arr(6)     ! minutes
        call w3fs21(idate5,nmind)
        t4dv = (real((nmind-iwinbgn),r_kind) + real(hdrh8arr(7),r_kind)*r60inv)*r60inv
        sstime = real(nmind,r_kind) + real(hdrh8arr(7),r_kind)*r60inv
        tdiff=(sstime-gstime)*r60inv
        if (l4dvar.or.l4densvar) then
           if (t4dv<zero .OR. t4dv>winlen) cycle read_loop
        else
           if (abs(tdiff)>twind) cycle read_loop
        endif

        crit0 = 0.01_r_kind
        timeinflat=6.0_r_kind
        call tdiff2crit(tdiff,ptime,ithin_time,timeinflat,crit0,crit1,it_mesh)
        call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis,it_mesh=it_mesh)
        if(.not. iuse)cycle read_loop

           rclrsky=bmiss
           rcldfrc=bmiss

           call ufbrep(lnbufr,dataahi,1,nchn,iret,'NCLDMNT')
!          dataahi(1,2) is high-peaking water vapor channel
!          for AHI CSR, clear-sky percentage are the same for all the channels
           if(dataahi(1,2)>= zero .and. dataahi(1,2) <= 100.0_r_kind ) then
              rclrsky=dataahi(1,2)
!             first QC filter out data with less clear sky fraction
              if ( rclrsky < 70.0_r_kind ) cycle read_loop
           end if

        call ufbrep(lnbufr,dataahisd,1,nchn,iret,'SDTB')
 
!       toss data if SDTB>1.3 
        do i=1,nchanl
           if(i==2 .or. i==3 .or. i==4) then   ! 3 water-vapor channels
              if(dataahisd(1,i)>1.3_r_kind) then
                 cycle read_loop
              end if
           end if
        end do

!       Locate the observation on the analysis grid.  Get sst and land/sea/ice
!       mask.  

!     isflg    - surface flag
!                0 sea
!                1 land
!                2 sea ice
!                3 snow
!                4 mixed                         


        call deter_sfc(dlat,dlon,dlat_earth,dlon_earth,t4dv,isflg,idomsfc,sfcpct, &
            ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr)

!       if (isflg >= 1) cycle read_loop   !!!test ocean only

!       Set common predictor parameters

!! comment the following line to add the clear filter for Himawari-8 AHI !!
!start by ZZMA

        !---secant of the sensor zenith angle  
        seca = 1.0_r_kind / COS( hdrh8arr(ilzah) * deg2rad )
        !---brightness temperature of channels used for regression
        bt_ts   = dataahibt(1,ts_ichan)
         
        !---difference in BTs water-window or window-water 
        !  (depends on channel selection)
        dbt_ts  = bt_ts(2)-bt_ts(1)

        !---calculate regression surface temperature using 
        !   empirical relation
        ts_reg = &
             ts_coef0 + ts_coef(1)*bt_ts(2) + ts_coef(2)*dbt_ts + &
             ts_coef(3)*dbt_ts*dbt_ts + ts_coef(4) * seca
        !---automatically reject freezing regression temperatures
        if ( ts_reg <= 273.00_r_kind) ts_reg = -10.
        
        !---two options with the split window test
        ! 1.) throw out observations with large SST_reg - SST_detersfc
        ! 2.) use thinning routines from satthin.F90 module :: checkob () 
        !     to pick "best" observation in thinning box

        !---Option 1.)
        ! the following line will through out any observation 
        ! with delta.ts > dts_threshold
        !  tsavg from deter_sfc
        sst_test = tsavg-ts_reg
!       if (abs(sst_test) >= dts_thresh) cycle read_loop

        !---Option 2.)
        !---or we can do this --use sathin module to select best pixels
        !   in thinning box
        !---larger values are bad
        ! sst_test = max(zero,sst_test)
        !---or i would prefer -- larger values still bad for whatever reason
        sst_test = ABS(sst_test)
        pred  = 15._r_kind*sst_test
        
!end by ZZMA

        crit1=crit1+rlndsea(isflg)
        call checkob(dist1,crit1,itx,iuse)
        if(.not. iuse)cycle read_loop

!       use NCLDMNT from chn7 (10.8 micron) as a QC predictor
!       add SDTB from chn7 as QC predictor
        pred=10.0_r_kind-dataahi(1,7)/10.0_r_kind+dataahisd(1,7)*10.0_r_kind
        
        crit1 = crit1 + pred
!       Compute "score" for observation.  All scores>=0.0.  Lowest score is "best"
        call finalcheck(dist1,crit1,itx,iuse)
        if(.not. iuse) cycle read_loop

!       Map obs to grids
        iscan = nint(hdrh8arr(ilzah))+1.001_r_kind ! integer scan position
!
!       interpolate NSST variables to Obs. location and get dtw, dtc, tz_tr
!
        if ( nst_gsi > 0 ) then
           tref  = ts(0)
           dtw   = zero
           dtc   = zero
           tz_tr = one
           if ( sfcpct(0) > zero ) then
              call gsi_nstcoupler_deter(dlat_earth,dlon_earth,t4dv,zob,tref,dtw,dtc,tz_tr)
           endif
        endif

!       Transfer information to work array
        data_all( 1,itx) = hdrh8arr(1)                ! satellite id
        data_all( 2,itx) = t4dv                       ! analysis relative time
        data_all( 3,itx) = dlon                       ! grid relative longitude
        data_all( 4,itx) = dlat                       ! grid relative latitude
        data_all( 5,itx) = hdrh8arr(ilzah)*deg2rad    ! satellite zenith angle (radians)
        data_all( 6,itx) = bmiss                      ! satellite azimuth angle (radians)
        data_all( 7,itx) = rclrsky                    ! clear sky amount
        data_all( 8,itx) = iscan                      ! integer scan position
        data_all( 9,itx) = hdrh8arr(iszah)            ! solar zenith angle
        data_all(10,itx) = bmiss                      ! solar azimuth angle
        data_all(11,itx) = sfcpct(0)                  ! sea percentage of
        data_all(12,itx) = sfcpct(1)                  ! land percentage
        data_all(13,itx) = sfcpct(2)                  ! sea ice percentage
        data_all(14,itx) = sfcpct(3)                  ! snow percentage
        data_all(15,itx)= ts(0)                       ! ocean skin temperature
        data_all(16,itx)= ts(1)                       ! land skin temperature
        data_all(17,itx)= ts(2)                       ! ice skin temperature
        data_all(18,itx)= ts(3)                       ! snow skin temperature
        data_all(19,itx)= tsavg                       ! average skin temperature
        data_all(20,itx)= vty                         ! vegetation type
        data_all(21,itx)= vfr                         ! vegetation fraction
        data_all(22,itx)= sty                         ! soil type
        data_all(23,itx)= stp                         ! soil temperature
        data_all(24,itx)= sm                          ! soil moisture
        data_all(25,itx)= sn                          ! snow depth
        data_all(26,itx)= zz                          ! surface height
        data_all(27,itx)= idomsfc + 0.001_r_kind      ! dominate surface type
        data_all(28,itx)= sfcr                        ! surface roughness
        data_all(29,itx)= ff10                        ! ten meter wind factor
        data_all(30,itx)= dlon_earth_deg              ! earth relative longitude (degrees)
        data_all(31,itx)= dlat_earth_deg              ! earth relative latitude (degrees)
        data_all(32,itx) = rcldfrc                    ! total cloud fraction from AHIASR
        do k=1,nchanl
           data_all(32+k,itx) = dataahisd(1,k)      ! BT standard deviation from AHICSR
        end do

        if(dval_use)then
           data_all(maxinfo-1,itx) = val_img
           data_all(maxinfo,itx) = itt
        end if

        if ( nst_gsi > 0 ) then
           data_all(maxinfo+1,itx) = tref         ! foundation temperature
           data_all(maxinfo+2,itx) = dtw          ! dt_warm at zob
           data_all(maxinfo+3,itx) = dtc          ! dt_cool at zob
           data_all(maxinfo+4,itx) = tz_tr        ! d(Tz)/d(Tr)
        endif

!       Transfer observation location and other data to local arrays

        do k=1,nchanl
           data_all(k+nreal,itx)=dataahibt(1,k)   ! test only for AHI channels:7-16
        end do

        nrec(itx)=irec

     enddo read_loop
  enddo read_msg
  call closbf(lnbufr)
  close(lnbufr)

  call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,&
     nele,itxmax,nread,ndata,data_all,score_crit,nrec)

! If no observations read, jump to end of routine.
  if (mype_sub==mype_root.and.ndata>0) then
     do n=1,ndata
        do k=1,nchanl
           if(data_all(k+nreal,n) > tbmin .and. &
              data_all(k+nreal,n) < tbmax)nodata=nodata+1
        end do
     end do
     if(dval_use .and. assim)then
       do n=1,ndata
          itt=nint(data_all(maxinfo,n))
          super_val(itt)=super_val(itt)+val_img
       end do
     end if

! Write final set of "best" observations to output file
     call count_obs(ndata,nele,ilat,ilon,data_all,nobs)
     write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
     write(lunout) ((data_all(k,n),k=1,nele),n=1,ndata)
  end if

! Deallocate local arrays
  deallocate(data_all,nrec)
  deallocate(dataahi,dataahibt,dataahisd)

! Deallocate satthin arrays
  call destroygrids

  if(diagnostic_reg.and.ntest>0) write(6,*)'READ_AHI:  ',&
     'mype,ntest,disterrmax=',mype,ntest,disterrmax

  return
end subroutine read_ahi