c==========================================================================
c
c NAME: grid_ch4
c
c PURPOSE:  Grids HIRS cloud height scan line data and bins into grid. Currently
c           grids cloud top pressure and temperature and emissivity. Separates
c           ascending and descending nodes
c
c INPUTS:
c     day_file.........File name of daily clear-sky file
c     jday.............julian day of year
c     iyr4.............4 digit year
c     nlon.............number of longitude grid points
c     nlat.............number of latitude grid points
c     max_cnt..........Maximum number of observations in time array
c     iparm............Parameter to grid
c     s_range(2).......Range of scan points (1-28)
c     gres.............grid resolution
c
c OUTPUT:
c     sum..............sum of the values
c     sumsq............sum squared of values
c     obs_time.........array of observation times in julian days from Jan 1, 1969.
c     icnt.............number of points.
c
c EXTERNALS:
c            hirs_grid.inc
c            recl_parms.inc
c            binit.f
c            date_to_julian69.f
c            nodalpass.f
c
c HISTORY:
c           Darren Jackson  CIRES/ETL  April 2002
c             - Version 3 computes array of observation times, September 2002
c             - Version 4 identifies asc/des passes
c
c============================================================================
      subroutine grid_ch4(day_file,jday,iyr4,nlon,nlat,max_cnt,
     .     iparm,s_range,sum,sumsq,obs_time,icnt,gres)

      implicit none
      include    'hirs_grid.inc'
      include    'recl_parms.inc'

      integer    nlon,nlat,max_cnt,iparm,s_range(2)
      integer    icnt(nlon,nlat,node),ispot(nspot),irec,ii,jj,k
      integer    obs_time(nlon,nlat,max_cnt,node)
      integer    jday,jday69,iyr4
      integer    itmp,ijday69,prvline,prvtime,prvn,n,line

      real       sum(nlon,nlat,node),sumsq(nlon,nlat,node)
      real       alt,fscale,xlon,xlat,gres,ref
      real       parm,prvlon

      character  day_file*95
      logical    bad_data,nolimb
c
c* HIRS CH scan line data array (1 record  = 46 bytes) 
c
      integer   itime
      integer*2 ilon,ilat,iszen,ialt,iline,ipct,itc,itsfc
      byte      isp,iref,iem,iclear
      integer*2 itb(7),ideltrad(4)
c
c* Inititalize
c
      data bad_data/.false./
c
c* Read and bin scan line data
c
      irec=1
      open(1,file=day_file,access='direct',recl=ch_rec)
 60   read(1,rec=irec,err=99) itime,ilon,ilat,iline,isp,
     .     iszen,ialt,iref,itb,ipct,iem,itc,itsfc,ideltrad,
     .     iclear

      xlon=fscale(ilon,180.,100.)
      line=iline

      if(isp .gt. 28) isp = isp - 28
      if(isp .ge. s_range(1) .and. isp .le. s_range(2)) then

         call date_to_julian69(iyr4,1,jday,jday69)
         itmp = nint((float(itime) /100. / 86400.) * 1000.)
         ijday69 = jday69 * 1000 + itmp
         xlat=fscale(ilat,0.,100.)
         call binit(gres,0,0,xlon,xlat,ii,jj)
c
c* Do not process first pixel so to find node direction
c
         if(irec .gt. 1) then 

            call nodalpass(itime,xlon,line,n,prvtime,prvlon,
     .           prvline,prvn)
c
c* Convert parameters for binning
c
            if(iparm .eq. 1) then
               parm = ipct * 1. ! Cloud top pressure
            elseif(iparm .eq. 2) then
               parm = fscale(itc,100.,100.) ! Cloud top temperature
            elseif(iparm .eq. 3) then
               parm = iem * 1.  ! IR emissivity
            elseif(iparm .ge. 4 .and. iparm .le. 8) then
               parm = fscale(itb(iparm-3),100.,100.) ! Channel 4-8 brightness temperature
            elseif(iparm .eq. 10) then
               parm = fscale(itb(6),100.,100.) ! Channel 10 brightness temperature
            elseif(iparm .eq. 12) then
               parm = fscale(itb(7),100.,100.) ! Channel 12 brightness temperature
            elseif(iparm .eq. 20) then
               parm = fscale(itsfc,100.,100.) ! Surface temperature
            elseif(iparm .ge. 24 .and. iparm .le. 27) then
               parm = fscale(ideltrad(iparm-23),0.,100.) ! Channel 4-8 Rclr-Rcld
            elseif(iparm .eq. 30) then
               if((iclear .eq. 0 .and. ipct .lt. 1000) .or. 
     .              (iclear .eq. 1 .and. iem .lt. 95)) then
                  parm = 1.
               else
                  parm = 0.
               endif
            endif
c     
c* Sum resulting parameters
c
            if(parm .gt. -998.) then
               sum(ii,jj,n)=sum(ii,jj,n)+parm
               sumsq(ii,jj,n)=sumsq(ii,jj,n)+parm**2
               icnt(ii,jj,n)=icnt(ii,jj,n)+1
               if(icnt(ii,jj,n) .le. max_cnt) then
                  obs_time(ii,jj,icnt(ii,jj,n),n)=ijday69
               else
                  write(*,*) 'Number of obs. exceeds max_cnt'
                  write(*,*) 'max_cnt,icnt = ',max_cnt,icnt(ii,jj,n) 
                  write(*,*) 'year = ',iyr4,' jday = ',jday
                  write(*,*) 'longitude,latitude = ',xlon,xlat
                  write(*,*) ' '
                  obs_time(ii,jj,1,n) = -1
               endif
            endif

         endif
      endif

      irec=irec+1
      prvn=n
      prvtime=itime
      prvlon=xlon
      prvline=line
      goto 60

 99   close(1)
      write(*,*) 'Read ',irec,' records in ',day_file
      return
      end