c=============================================================================
c
c NAME: hcloud_lt
c
c PURPOSE: Performs second cloud test on the clear-sky data so to remove
c          regions of persistent cloud cover.  
c
c INPUTS: (via hcloud_lt.input)
c                  
c         isat     ->  Satellite number
c         jday     ->  julian day from jan 1, 1969
c 
c INPUT FILES: TBC and CS1 data files
c
c OUTPUT FILES: CS data files
c
c EXTERNALS:  (description given in README)
c             julian69_to_date.f
c             binit.f
c
c HISTORY: Darren Jackson ETL/CIRES April 2002
c
c============================================================================

      program hcloud_lt

c
      implicit none
      include 'hirs1c.inc'
      include 'recl_parms.inc'

      byte      lmask(nlon,nlat)
      integer   isat,iyr2,bday,eday,jday,ispot,line,irec,orec,n,fday30
      integer   iyr4,imn,idy,year,month,iday,rem,bjday
      integer   ii,jj,it,igl,ics,ifl,csf,npen,ipen,i,j,two_digit_year
      integer   ejday,mon,biyr2,eiyr2
      real      xlon,xlat,szen,alt,refl,tb(nchn2),tblc,del1(2)
      real      xtime,tovs_hr,perc,fscale,tblt,eres
      character grid_file*50,dgrid_path*8,tmp_file*70
      character cs_file*50,l_file*70
      integer*2 iscale,ogrid(nlon,nlat,ntim)
      logical   around,clear
c      integer*2 elev(3600,1800)
c
c* HIRS scan line data array
c
      integer   itime
      integer*2 ilon,ilat,iszen,ialt,iline
      byte      isp,iqc,isf,iref
      integer*2 itb(19)
c
c* Data values
c
      include 'datapaths.inc'
      data igl/10/,ics/11/,ifl/20/
      data del1/3.,6./
      data dgrid_path/'../grid/'/
      data eres/0.1/
c
c*  Retrieve input parameters
c
      open(10,file='hcloud_lt.input')
      read(10,*) isat
      read(10,*) jday
      close(10)

      open(41,file='hcloud_lt.err',access='append')
c
c* Retrieve land mask data
c
      l_file=misc_path//'lmask_0.5_lsc_time.dat'
      open(10,file=l_file,recl=nlon*nlat*1,access='direct')
      read(10,rec=1) lmask
      close(10)
c
c* Retrieve elevation data and swap int*2 data
c
c      l_file=misc_path//'elev_0.1d.dat'
c      open(10,file=l_file,recl=3600*1800*2,access='direct')
c      read(10,rec=1) elev
c      close(10)
c      do i=1,3600
c         do j=1,1800
c            call swap16(elev(i,j),1)
c         enddo
c      enddo
c
c* Find nearest pentad values
c
      rem = mod(jday,5)
      if(rem .eq. 0) rem = 5
      bjday = jday - rem + 1
      ejday = bjday + 4
c
c* Convert julian day to day of year
c
      call julian69_to_date(1,bjday,iyr4,mon,bday)
      biyr2=two_digit_year(iyr4)
      call julian69_to_date(1,ejday,iyr4,mon,eday)
      eiyr2=two_digit_year(iyr4)
      call julian69_to_date(1,jday,iyr4,mon,iday)
      iyr2=two_digit_year(iyr4)      
c
c* Open clear-sky composite data
c
      around=.false.
      write(grid_file,49) 'TBC.N',isat,'.Y',biyr2,'.B',bday,
     .     '.Y',eiyr2,'.E',eday,'.DAT'
 49   format(a5,i2.2,a2,i2.2,a2,i3.3,a2,i2.2,a2,i3.3,a4)
      inquire(file=dgrid_path//grid_file,exist=around)

      if(around) then

         open(unit=igl,file=dgrid_path//grid_file,
     .        access='direct',recl=nlon*nlat*ntim*2)
         read(igl,rec=1) ogrid
         close(igl)
c
c*  Open temporary clear-sky data file
c
         around=.false.
         write(tmp_file,50) 'HIRS.N',isat,'.Y',
     .        iyr2,'.D',iday,'.V22.CS1'
 50      format(a6,i2.2,a2,i2.2,a2,i3.3,a8)
         inquire(file=tmp_path//tmp_file,exist=around)
         if(around) then
            open(unit=ics,file=tmp_path//tmp_file,access='direct',
     .           recl=cs1_rec)
c
c* Open clear-sky output file
c     
            write(cs_file,51) 'HIRS.N',isat,'.Y',
     .           iyr2,'.D',iday,'.V22.CS'
 51         format(a6,i2.2,a2,i2.2,a2,i3.3,a7)
            open(unit=ifl,file=cs_path//cs_file,access='direct',
     .           recl=cs_rec)
c
c*  Loop through each FOV
c
            irec=1
            orec=1
 80         read(ics,rec=irec,err=99) xtime,xlon,xlat,szen,alt,
     .           ispot,line,csf,refl,tblc,(tb(n), n=1,nchn2)
            irec=irec+1
c
c* Find diurnal time index
c
            tovs_hr=xtime/3600.
            it=ifix(tovs_hr/tinc)+1
            if(it .ge. itmax) it=1
c
c* find lat,lon indices
c
            call binit(gres,1,1,xlon,xlat,ii,jj)
c            
c* Perform long-term cloud test
c
            clear=.true.
            tblt=fscale(ogrid(ii,jj,it),100.,100.)
            if(lmask(ii,jj) .ge. 1) then ! land or coastal region
               if(tblt-tblc .gt. del1(2)) clear=.false.
            else                ! ocean
               if(tblt-tblc .gt. del1(1)) clear=.false.
            endif
c
c* Remove observations above 1500m elevation
c
c            if(clear) then
c               call binit(eres,1,1,xlon,xlat,ii,jj)
c               if(elev(ii,jj) .gt. 1500) clear=.false.
c            endif
c
c* create output data structure
c
            if(clear) then
               itime=nint(xtime*100.)
               ilon=iscale(xlon,180.,100.)
               ilat=iscale(xlat,0.,100.)
               iszen=iscale(szen,0.,100.)
               ialt=iscale(alt,0.,10.)
               isp=ispot
               iline=line
c               iqc=oqc
               isf=csf
               if(refl .gt. 1.) then
                  if(refl .eq. missing) then
                     refl=2.55
                  else
                     refl=2.54
                  endif
               endif
               iref=nint(refl*100.)
               do n=1,nchn2
                  if(tb(n) .eq. missing) tb(n)=-99.
                  itb(n)=iscale(tb(n),100.,100.)
               enddo
               write(ifl,rec=orec) itime,ilon,ilat,iline,isp,
     .              iszen,ialt,isf,iref,itb
               orec=orec+1
            endif
            goto 80
            
 99         continue
            perc=float(orec-1)/float(irec-1)*100.
            write(*,*) isat,iyr4,jday,irec-1,orec-1,perc
            
         else
            write(41,*) 'No short term clear-sky swath file'
            write(41,*) 'Jday = ',jday
         endif
      else
         write(41,*) 'No clear-sky composite file for period'
         write(41,*) 'jday = ',jday,' bday = ',bday,' eday = ',eday
      endif


      close(41)
      end