c==============================================================================
c
c NAME: hcloudcomp
c
c PURPOSE: Program produces a HIRS chn 8 IR clear-sky brightness temperature 5-day 
c          grid by using 5-, 15-, and 35-day composite grids.
c
c INPUTS: (via hcloudcomp.input)
c         isat  -> satellite number
c         jday  -> julian day from Jan 1, 1969
c
c INPUT DATA:
c                5-day, 15-day and 35-day *.GRD files
c
c OUTPUT DATA:
c                TBC* data files
c
c EXTERNALS:
c            hirs1c.inc
c            julian69_to_date.f
c            tb_clear.f
c            two_digit_year.f
c
c HISTORY:  Darren Jackson ETL/CIRES April 2002
c==============================================================================

      program hcloudcomp

      implicit none
      include 'hirs1c.inc'

      byte      lmask(nlon,nlat)
      integer*2 iscale,ogrid(nlon,nlat,ntim,2),tbid
      integer   i,j,k,l,n
      integer   isat,iyr2,jday,iyr4,mon,day,two_digit_year
      integer   bday,eday,biyr2,eiyr2,bjday,ejday
      real      grid30(nlon,nlat,ntim,bs),grid15(nlon,nlat,ntim,bs)
      real      grid5(nlon,nlat,ntim,bs),igrid(nlon,nlat,ntim)
      real      tbclear,tavgst,tavglt,tmaxst,tmaxlt,nclrst,nclrlt
      character misc_path*38,l_file*70,grid_file*50,dgrid_path*8

      logical   around

      data dgrid_path/'../grid/'/
      data misc_path/'/data/dlj/tovs1b/pathfinder/misc_data/'/
c
c*  Retrieve input parameters
c
      open(10,file='hcloudcomp.input')
      read(10,*) isat
      read(10,*) jday
      close(10)
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* Convert julian day to day of year
c
      bjday=jday-17
      ejday=jday+17
      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)
c
c* Retrieve 35-day grid file
c
      around=.false.
      write(grid_file,50) 'HIRS.N',isat,'.Y',biyr2,'.B',bday,
     .     '.Y.',eiyr2,'.E',eday,'.V22.GD'
      inquire(file=dgrid_path//grid_file,exist=around)

      if(around) then
         open(unit=20,file=dgrid_path//grid_file,
     .        access='direct',recl=nlon*nlat*ntim*4)
         do l=1,bs
            read(20,rec=l,err=99) igrid
            do i=1,nlon
               do j=1,nlat
                  do k=1,ntim
                     grid30(i,j,k,l)=igrid(i,j,k)
                  enddo
               enddo
            enddo
         enddo
         close(20)
 50      format(a6,i2.2,a2,i2.2,a2,i3.3,a2,i2.2,a2,i3.3,a7)
      else
         write(*,*) '35-day grid does not exist'
         write(*,*) 'No TBC grid file created for'
         write(*,*) 'satellite = ',isat,' jday = ',jday
         stop
      endif
c
c* Retrive 15-day grid files
c     
      bjday=jday-7
      ejday=jday+7
      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)
      
      write(grid_file,50) 'HIRS.N',isat,'.Y',biyr2,'.B',bday,
     .     '.Y',eiyr2,'.E',eday,'.V22.GD'

      inquire(file=dgrid_path//grid_file,exist=around)
      if(around) then
         open(unit=20,file=dgrid_path//grid_file,
     .        access='direct',recl=nlon*nlat*ntim*4)
         do l=1,bs
            read(20,rec=l,err=99) igrid
            do i=1,nlon
               do j=1,nlat
                  do k=1,ntim
                     grid15(i,j,k,l)=igrid(i,j,k)
                  enddo
               enddo
            enddo
         enddo
         close(20)
      else
         write(*,*) '15-day grid does not exist'
         write(*,*) 'No TBC grid file created for'
         write(*,*) 'satellite = ',isat,' jday = ',jday
         stop
      endif
c
c* Retrieve 5-day grid files
c
      bjday=jday-2
      ejday=jday+2
      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)

      write(grid_file,50) 'HIRS.N',isat,'.Y',biyr2,'.B',bday,
     .     '.Y',eiyr2,'.E',eday,'.V22.GD'

      inquire(file=dgrid_path//grid_file,exist=around)
      if(around) then
         open(unit=21,file=dgrid_path//grid_file,
     .        access='direct',recl=nlon*nlat*ntim*4)
         do l=1,bs
            read(21,rec=l,err=99) igrid
            do i=1,nlon
               do j=1,nlat
                  do k=1,ntim
                     grid5(i,j,k,l)=igrid(i,j,k)
                  enddo
               enddo
            enddo
         enddo
         close(21)
      else
         write(*,*) '5-day grid does not exist'
         write(*,*) 'No TBC grid file created for'
         write(*,*) 'satellite = ',isat,' jday = ',jday
         stop
      endif
c
c*  Compute 5-day clear-sky Tb
c
      do i=1,nlon
         do j=1,nlat
            if(lmask(i,j) .ge. 1) then !land or coastal region
               do k=1,ntim
                  ogrid(i,j,k,1)=iscale(-99.,100.,100.)
                  nclrst=grid5(i,j,k,5)
                  nclrlt=grid15(i,j,k,5)
                  if(nclrst .ge. 1.) then
                     tbclear=-99.
                     tavgst=grid5(i,j,k,4)/grid5(i,j,k,5)
                     tavglt=grid15(i,j,k,4)/grid15(i,j,k,5)
                     tmaxst=grid5(i,j,k,1)
                     tmaxlt=grid15(i,j,k,1)
                     call tb_clear(lmask(i,j),nclrst,nclrlt,
     .                    tavgst,tavglt,tmaxst,tmaxlt,tbclear,
     .                    tbid)
                     ogrid(i,j,k,1)=iscale(tbclear,100.,100.)
                     ogrid(i,j,k,2)=tbid
                  endif
               enddo
            else                !oceanic
               do k=1,ntim
                  ogrid(i,j,k,1)=iscale(-99.,100.,100.)
                  nclrst=grid15(i,j,k,5)
                  nclrlt=grid30(i,j,k,5)
                  if(nclrst .ge. 1.) then
                     tbclear=-99.
                     tavgst=grid15(i,j,k,4)/grid15(i,j,k,5)
                     tavglt=grid30(i,j,k,4)/grid30(i,j,k,5)
                     tmaxst=grid15(i,j,k,1)
                     tmaxlt=grid30(i,j,k,1)
                     call tb_clear(lmask(i,j),nclrst,nclrlt,
     .                    tavgst,tavglt,tmaxst,tmaxlt,tbclear,
     .                    tbid)
                     ogrid(i,j,k,1)=iscale(tbclear,100.,100.)
                     ogrid(i,j,k,2)=tbid
                  endif
               enddo
            endif
         enddo
      enddo

c
c* Output the resulting grid file
c
      write(grid_file,51) 'TBC.N',isat,'.Y',biyr2,'.B',bday,
     .     '.Y',eiyr2,'.E',eday,'.DAT'
      open(unit=20,file=dgrid_path//grid_file,
     .     access='direct',recl=nlon*nlat*ntim*2*2)
      write(20,rec=1) ogrid
      close(20)
 51   format(a5,i2.2,a2,i2.2,a2,i3.3,a2,i2.2,a2,i3.3,a4)

      stop
     
 99   write(*,*) 'Read error in one grid file'
      write(*,*) 'No TBC file created for'
      write(*,*) 'Satellite = ',isat,' jday = ',jday

      end
