subroutine read_tcps(nread,ndata,nodata,infile,obstype,lunout,sis,nobs)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_tcps                   read tcvitals ascii file
!   prgmmr: kleist             org: np20                date: 2009-02-02
!
! abstract:  This routine reads the ascii TC psmin data file 
!
!
! program history log:
!   2009-02-02  kleist
!   2010-09-08  treadon - add station_id and destroy_tcv_card; increase
!                         maxdat to 10; remove i_kind suffix from integer 
!                         constants; 
!   2013-01-26  parrish - change from grdcrd to grdcrd1
!   2015-10-01  guo     - consolidate use of ob location (in deg)
!
!   input argument list:
!     infile   - unit from which to read ascii file
!     obstype  - observation type to process
!     lunout   - unit to which to write data for further processing
!
!   output argument list:
!     nread    - number of bogus data read
!     ndata    - number of bogus data 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,i_kind,r_double
  use gridmod, only: nlat,nlon,rlats,rlons,regional,tll2xy
  use constants, only: deg2rad,zero,one_tenth,one
  use convinfo, only: nconvtype,ictype,icuse
  use obsmod, only: ianldate
  use tcv_mod, only: get_storminfo,numstorms,stormlat,stormlon,stormpsmin,stormdattim,&
       centerid,stormid,destroy_tcv_card,tcp_refps,tcp_width,tcp_ermin,tcp_ermax
  use gsi_4dvar, only: time_4dvar
  use mpimod, only:npe
  implicit none

! Declare passed variables
  character(len=*),intent(in  ) :: obstype,infile
  character(20)  ,intent(in   ) :: sis
  integer(i_kind),intent(in   ) :: lunout
  integer(i_kind),intent(inout) :: nread,ndata,nodata
  integer(i_kind),dimension(npe),intent(inout) :: nobs

! Declare local parameters
  real(r_kind),parameter:: r360=360.0_r_kind
  integer(i_kind),parameter:: maxobs=2e6
  integer(i_kind),parameter:: maxdat=10

! Declare local variables
  real(r_double) rstation_id
  character(8) station_id
  equivalence(rstation_id,station_id)

  real(r_kind) dlat,dlon,dlat_earth,dlon_earth
  real(r_kind) dlat_earth_deg,dlon_earth_deg
  real(r_kind),allocatable,dimension(:,:):: cdata_all
  real(r_kind) ohr,olat,olon,psob,pob,oberr,usage,toff
  real(r_kind) psdif,alpha

  integer(i_kind) i,k,lunin,nc
  integer(i_kind) ilat,ilon,ikx,nreal,nchanl,noutside,nmrecs

  logical endfile, outside

  data lunin / 10 /

!**************************************************************************
! Initialize variables
  nmrecs=0
  nreal=maxdat
  nchanl=0
  noutside=0
  ilon=2
  ilat=3
  endfile=.false.

  allocate(cdata_all(maxdat,maxobs))

  open(lunin,file=trim(infile),form='formatted',action='read')
  call datelen(10)

! Big loop over binary file

  call get_storminfo(lunin)

  write(6,*) 'READ_TCPS:  IANLDATE = ',ianldate

  do i=1,numstorms
     nmrecs=nmrecs+1
     nread=nread+1

! Set ikx here...pseudo-mslp tc_vitals obs are assumed to be type 111
     do nc=1,nconvtype
        if (ictype(nc)==112) ikx=nc
     end do

! Set usage variable
     usage = zero
     if(icuse(ikx) <= 0)usage=100._r_kind

     if (stormdattim(i)/=ianldate) then
        write(6,*) 'READ_TCPS:  IGNORE TC_VITALS ENTRY # ',i
        write(6,*) 'READ_TCPS:  MISMATCHED FROM ANALYSIS TIME, OBS / ANL DATES = ',stormdattim(i),ianldate
        cycle
     end if

! Set center and storm id (only used in diagnostic file)
  station_id = centerid(i)(1:4) // '_' // stormid(i)(1:3)

! Observation occurs at analysis time as per date check above
! Set observation lat, lon, mslp, and default obs-error
     call time_4dvar(ianldate,toff)
!    write(6,*)'READ_TCPS: bufr file date is ',ianldate
!    write(6,*)'READ_TCPS: time offset is ',toff,' hours.'
     ohr=toff
     olat=stormlat(i)
     olon=stormlon(i)
     psob=stormpsmin(i)      ! in mb
     psdif=tcp_refps-psob    ! in mb
     alpha=max(min(psdif/tcp_width,one),zero)
     oberr=tcp_ermin+(tcp_ermax-tcp_ermin)*alpha

! Make sure the psob is reasonable
     if ( (psob<850._r_kind) .or. (psob>1025._r_kind) )then
        usage=100._r_kind
     end if

     if (olon >= r360) olon=olon-r360
     if (olon < zero)  olon=olon+r360
     dlat_earth_deg = olat
     dlon_earth_deg = olon
     dlat_earth = olat * deg2rad
     dlon_earth = olon * deg2rad
     if(regional)then
        call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
        if (outside) then
           noutside=noutside+1
           cycle
        endif
     else 
        dlat = dlat_earth
        dlon = dlon_earth
        call grdcrd1(dlat,rlats,nlat,1)
        call grdcrd1(dlon,rlons,nlon,1)
     end if

! Extract observation.
     ndata=min(ndata+1,maxobs)
     nodata=min(nodata+1,maxobs)

! convert pressure (mb) to log(pres(cb))
     pob=one_tenth*psob
 
     cdata_all(1,ndata)=oberr*one_tenth       ! obs error converted to cb
     cdata_all(2,ndata)=dlon                  ! grid relative longitude
     cdata_all(3,ndata)=dlat                  ! grid relative latitude
     cdata_all(4,ndata)=pob                   ! pressure in cb 
     cdata_all(5,ndata)=toff                  ! obs time (analyis relative hour)
     cdata_all(6,ndata)=ikx                   ! obs type
     cdata_all(7,ndata)=dlon_earth_deg        ! earth relative longitude (degrees)
     cdata_all(8,ndata)=dlat_earth_deg        ! earth relative latitude (degrees)
     cdata_all(9,ndata)=usage                 ! usage parameter
     cdata_all(10,ndata)=rstation_id          ! storm name (centerid_stormid)

! End of loop over number of storms
  end do

  write(6,*) 'READ_TCPS:  NUMBER OF OBS READ IN = ', ndata
  write(6,*) 'READ_TCPS: # out of domain =', noutside

! Write observations to scratch file

  call count_obs(ndata,maxdat,ilat,ilon,cdata_all,nobs)
  write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
  write(lunout) ((cdata_all(k,i),k=1,maxdat),i=1,ndata)

! Close unit
  deallocate(cdata_all)
  call destroy_tcv_card
  close(lunin)

! End of routine
  return
end subroutine read_tcps