c----------------------------------------------------------------------------
c     NAME:    read_scanline_c2
c
c     PURPOSE: Example program for reading HIRS all-sky c2 data.
c
c
c     LANGUAGE: FORTRAN-77
c
c     COMPILATION SEQUENCE:
c                           f77 -o read_scanline_c2.x read_scanline_c2.f
c     INPUTS:
c             Prompted inputs when executing code
c
c     OUTPUT:
c           xtime..........GMT Time in seconds
c           xlon...........longitude in degrees East from GM (0-359.99)
c           xlat...........latitude in degrees (-90S to 90N)
c           xszen..........solar zenith angle (0 to 90 -> 90 is night)
c           xalt...........satellite altitude in km
c           iline..........line number
c           ispot..........scan spot position (1-56)
c           ref............reflectance (0 - 1 valid; 2.55 and 2.54 missing)
c           tb(19).........no limb bright. temp. in Kelvin.  Array index
c                          matches channel number.
c
c
c   NOTES:  
c         1. All-sky data were produced by Little Endian computer.
c
c   CREATION HISTORY:
C                     Darren Jackson CIRES/ETL August 2004     
c----------------------------------------------------------------------------
c
      program read_scanline_c2

      implicit none

      integer    nspot,nchn2,nchn3,nang,ncoef
      parameter  (nspot=1,nchn2=19,nchn3=12,nang=6,ncoef=15)

      integer    irec,k,n,isat,iyr2,jday,i,ip,iprint
      real       fscale
      character  day_file*70,ctmp*1,satid(12)*3,i_path*50,r_path*50
      logical    limb,bad_data,mis_data,around,little
c
c* HIRS scan line input data
c
      integer   itime
      integer*2 ilon,ilat,iszen,ialt,iline
      byte      isp,isf,iref
      integer*2 itb(nchn2)
c
c* HIRS scan line output data
c
      real      xtime,xlon,xlat,xszen,xalt
      real      tb(nspot,nchn2),tblc(nspot,nchn2),ref
      integer   ispot,line
c
c* Limb correction coefficients
c
c    Night coefficients
      integer    chnord(ncoef-1,nang,nchn2),npred(nang,nchn2)
      real       rcoef(ncoef,nang,nchn2)
c    Day coefficients
      integer    chnord2(ncoef-1,nang,nchn3),npred2(nang,nchn3)
      real       rcoef2(ncoef,nang,nchn3)
c
c* Variables for byte swapping
c
      integer*2 i2
      integer*4 i4
c
c* Data arrays
c
      data satid/'N05','N06','N07','N08','N09','N10','N11',
     .           'N12','N13','N14','N15','N16'/
      data bad_data/.false./
      data mis_data/.false./
c
c-------------------------------------------------------------------
c        MAIN PROGRAM
c-------------------------------------------------------------------
c
c* Input parameters
c
      write(*,*) 'Input satellite number (5-14):'
      read(*,*) isat
      write(*,*) 'Input two digit year;'
      read(*,*) iyr2
      write(*,*) 'Input julian day:'
      read(*,*) jday
      write(*,*) 'Is operating system Little Endian? (y or n)'
      write(*,*) '(Little Endian examples: INTEL-based PC, DEC)'
      read(*,*) ctmp
      if(ctmp .eq. 'y') then
         little=.true.
      elseif(ctmp .eq. 'n') then
         little=.false.
      else
         write(*,*) 'swap answer not y or n'
         write(*,*) 'Terminate program'
         stop
      endif
      write(*,*) 'Scan line data path:'
      read(*,'(a)') i_path
      ip=index(i_path,' ')-1
      write(*,*) 'Print to screen - 1, Print to file - 2'
      read(*,*) iprint
      if(iprint .ne. 1) open(22,file='read_scanline_v2.dat')
c
c* Load HIRS limb regression coefficients
c
      if(limb) then
         call open_regcoef(r_path,isat,little,rcoef,chnord,npred,rcoef2,
     .        chnord2,npred2)
         write(*,*) 'Completed retrieval of regression coefficients'
      endif
c
c* Read scan line data
c
      write(day_file,50) i_path(1:ip),'/HIRS.',satid(isat-4),'.Y',iyr2,
     .     '.D',jday,'.V22.CS'
 50   format(a,a6,a3,a2,i2,a2,i3.3,a7)
      inquire(file=day_file,exist=around)
      if(around) then
         open(1,file=day_file,access='direct',recl=55)
         irec=1
         write(*,*) 'Finished opening ',day_file
         read(1,rec=irec,err=999) itime,ilon,ilat,iline,isp,iszen,
     &     ialt,iref,(itb(k),k=4,8),itb(10),itb(12),iem,itc,
     &     itsfc,(ideltrad(k),k=1,4),iclear(iorbit,iline,isp)
     &    ,(ipctrto(k),k=1,5),iused

c 60      read(1,rec=irec,err=99) itime,ilon,ilat,iline,isp,iszen,ialt,
c     .        isf,iref,itb

         if(.not. little) then
            call swap32 (itime,i4)
            itime = i4
            call swap16 (ilon,i2)
            ilon = i2
            call swap16 (ilat,i2)
            ilat = i2
            call swap16 (iline,i2)
            iline=i2
            call swap16 (iszen,i2)
            iszen = i2
            call swap16 (ialt,i2)
            ialt = i2
            do k=1,nchn2
               call swap16 (itb(k),i2)
               itb(k) = i2
            enddo
         endif
         xtime=itime/100.
         xlon=fscale(ilon,180.,100.)
         xlat=fscale(ilat,0.,100.)
         xszen=fscale(iszen,0.,100.)
         xalt=fscale(ialt,0.,10.)
         ispot=isp
         line=iline
         do k=1,nchn2
            tb(1,k)=fscale(itb(k),100.,100.)
         enddo
         ref=int(iref)/100.
         if(ref .lt. 0.) ref=ref+2.56
c     
c* Apply limb correction
c
         if(limb) then
            call hlimb(xszen,ispot,xalt,bad_data,rcoef,
     .           chnord,npred,rcoef2,chnord2,npred2,tb,tblc)
            if(xszen .lt. 90.) then
               do i=13,19
                  tblc(1,i)=tb(1,i)
               enddo
            endif
         endif
c     
c* write data record
c
         if(iprint .eq. 1) then
            write(*,*) irec,xtime,xlon,xlat,xszen,xalt,line,ispot
            if(limb) then
               do n=1,nchn2
                  write(*,*) 'channel=',n,tb(1,n),tblc(1,n)
               enddo
            else
               do n=1,nchn2
                  write(*,*) 'channel=',n,tb(1,n)
               enddo
            endif
            write(*,'(a12,2x,f6.2)') 'channel= 20',ref
            write(*,*) ' '
         else
            write(22,*) irec,xtime,xlon,xlat,xszen,xalt,line,ispot
            if(limb) then
               do n=1,nchn2
                  write(22,*) 'channel=',n,tb(1,n),tblc(1,n)
               enddo
            else
               do n=1,nchn2
                  write(22,*) 'channel=',n,tb(1,n)
               enddo
            endif
            write(22,'(1x,a12,4x,f6.2)') 'channel=  20',ref
            write(22,*) ' '
         endif
c
c* read next data record
c
         irec=irec+1
         goto 60

      else
         write(*,*) 'Input file ',day_file,' does not exist'
      endif

 99   close(1)
      if(iprint .ne. 1) close(22)
      end
c
c---------------------------------------------------------
c    END OF MAIN PROGRAM
c----------------------------------------------------------
c
      subroutine open_regcoef(r_path,isat,little,rcoef,chnord,npred,
     .     rcoef2,chnord2,npred2)
c
c*  Opens and reads HIRS limb regression coefficients.
c   Coefficients for all channels are computed.
c   Retrieves coefficients for night and day obserations.  
c
c   INPUT:
c          r_path.....directory path for coefficient data
c          isat.......Satellite number
c          little.....Little or Big Endian computer
c   OUTPUT:
c          rcoef.......Regression coefficients for night observations
c          chnord......Channel order for reg. coefficients (night)
c          npred.......Number of predictors (night)
c          rcoef2......Regression coefficients for day obserations
c          chnord2.....Channel order for reg. coefficients (day)
c          npred2......Number of predictors (day)
c
      implicit none

      integer   ncoef,nang,nchn2,nchn3
      parameter (ncoef=15,nang=6,nchn2=19,nchn3=12)

      real      rcoef(ncoef,nang,nchn2)
      integer   icoef(ncoef,nang,nchn2)
      integer   chnord(ncoef-1,nang,nchn2),npred(nang,nchn2)
      real      rcoef2(ncoef,nang,nchn3)
      integer   icoef2(ncoef,nang,nchn3)
      integer   chnord2(ncoef-1,nang,nchn3),npred2(nang,nchn3)
      integer   i,j,k,isat,iok,rp,i4
      logical   little,around

      character coef_file*16,chnord_file*18,npred_file*17,suffix2*8
      character infile*70,r_path*50,suffix1*4,suffix3*8,suffix4*12

      coef_file='hirs_limb_coef_n'
      chnord_file='hirs_limb_chnord_n'
      npred_file='hirs_limb_npred_n'
      suffix1='.dat'
      suffix2='.lon.dat'
      suffix3='.lir.dat'
      suffix4='.lon.lir.dat'

      rp=index(r_path,' ')-1
      iok=10

c*  Open night coefficient file
      write(infile,'(a,a16,i2.2,a8)') r_path(1:rp),coef_file,isat,
     .     suffix2
      inquire(file=infile,exist=around)
      if(around) then
         open(iok,file=infile,access='direct',recl=ncoef*nang*nchn2*4)
         read(iok,rec=1) icoef
         close(iok)
         do i=1,ncoef
            do j=1,nang
               do k=1,nchn2
                  if(little) then
                     call swap32(icoef(i,j,k),i4)
                     icoef(i,j,k)=i4
                  endif
                  rcoef(i,j,k)=icoef(i,j,k)/100000. 
               enddo
            enddo
         enddo
      else
         write(*,*) 'Input file ',infile,' does not exist'
         stop
      endif
c*  Open day coefficient file
      write(infile,'(a,a16,i2.2,a12)') r_path(1:rp),coef_file,isat,
     .     suffix4
      inquire(file=infile,exist=around)
      if(around) then
         open(iok,file=infile,access='direct',recl=ncoef*nang*nchn3*4)
         read(iok,rec=1) icoef2
         close(iok)
         do i=1,ncoef
            do j=1,nang
               do k=1,nchn3
                  if(little) then
                     call swap32(icoef2(i,j,k),i4)
                     icoef2(i,j,k)=i4
                  endif
                  rcoef2(i,j,k)=icoef2(i,j,k)/100000.
               enddo
            enddo
         enddo
      else
         write(*,*) 'Input file ',infile,' does not exist'
         stop
      endif

c*  Open night channel order file
      write(infile,'(a,a18,i2.2,a4)') r_path(1:rp),chnord_file,isat,
     .     suffix1
      inquire(file=infile,exist=around)
      if(around) then
         open(iok,file=infile,access='direct',
     .        recl=(ncoef-1)*nang*nchn2*4)
         read(iok,rec=1) chnord
         close(iok)
         if(little) then
            do i=1,ncoef-1
               do j=1,nang
                  do k=1,nchn2
                     call swap32(chnord(i,j,k),i4)
                     chnord(i,j,k)=i4
                  enddo
               enddo
            enddo
         endif
      else
         write(*,*) 'Input file ',infile,' does not exist'
         stop
      endif
c*  Open day channel order file
      write(infile,'(a,a18,i2.2,a8)') r_path(1:rp),chnord_file,isat,
     .     suffix3
      inquire(file=infile,exist=around)
      if(around) then
         open(iok,file=infile,access='direct',
     .        recl=(ncoef-1)*nang*nchn3*4)
         read(iok,rec=1) chnord2
         close(iok)
         if(little) then
            do i=1,ncoef-1
               do j=1,nang
                  do k=1,nchn3
                     call swap32(chnord2(i,j,k),i4)
                     chnord2(i,j,k)=i4
                  enddo
               enddo
            enddo
         endif
      else
         write(*,*) 'Input file ',infile,' does not exist'
         stop
      endif

c*  Open number of predictors file for night
      write(infile,'(a,a17,i2.2,a4)') r_path(1:rp),npred_file,isat,
     .     suffix1
      inquire(file=infile,exist=around)
      if(around) then
         open(iok,file=infile,access='direct',recl=nang*nchn2*4)
         read(iok,rec=1) npred
         close(iok)
         if(little) then
            do j=1,nang
               do k=1,nchn2
                  call swap32(npred(j,k),i4)
                  npred(j,k)=i4
               enddo
            enddo
         endif
      else
         write(*,*) 'Input file ',infile,' does not exist'
         stop
      endif
c*  Open number of predictors file for day
      write(infile,'(a,a17,i2.2,a8)') r_path(1:rp),npred_file,isat,
     .     suffix3
      inquire(file=infile,exist=around)
      if(around) then
         open(iok,file=infile,access='direct',recl=nang*nchn3*4)
         read(iok,rec=1) npred2
         close(iok)
         if(little) then
            do j=1,nang
               do k=1,nchn3
                  call swap32(npred2(j,k),i4)
                  npred2(j,k)=i4
               enddo
            enddo
         endif
      else
         write(*,*) 'Input file ',infile,' does not exist'
         stop
      endif

      return
      end
c
c--------------------------------------------------------------------
c
      subroutine hlimb(szen,isp,alt,bad,rcoef,chnord,npred,rcoef2,
     .     chnord2,npred2,tb,tblc)
c
c*    Computes limb correction using regression coefficient data.
c     Uses coefficients for both day and night observations.     
c
c     INPUT:
c            szen(nspot)...............Solar zenith angle (0-90 degrees)
c            isp(nspot)................HIRS scan position (1-56)
c            alt(nspot)................altitude
c            bad(nspot)................Bad values
c            rcoef(ncoef,nang,nchn2)....regression coefficients (night)
c            chnord(ncoef-1,nang,nchn2).Channel number for coefficients (night)
c            npred(nang,nchn2)..........Number of predictors (night)
c            rcoef2(ncoef,nang,nchn3)...regression coefficients (day)
c            chnord2(ncoef-1,nang,nchn3)..Channel number for coefficients (day)
c            npred2(nang,nchn3)...........Number of predictors (day)
c            tb(nspot,nchn)...............HIRS Tbs for all HIRS channels
c
c    OUTPUT:
c           tblc(nspot,nchn)....HIRS limb corrected TBs for (1-12) for day.
c                               HIRS limb corrected TBs for (1-19) for night
c
      implicit none

      integer   ncoef,nang,nchn2,nchn3,nspot
      real      xmiss
      parameter (xmiss=-99.,ncoef=15,nang=6,nchn2=19,nchn3=12,nspot=1)

      integer i,j,n,isp(nspot)
      integer chnord2(ncoef-1,nang,nchn3),npred2(nang,nchn3)
      integer chnord(ncoef-1,nang,nchn2),npred(nang,nchn2)
      integer indx1,indx2

      real    tb(nspot,nchn2),tblc(nspot,nchn2),cosx(nang),y(2),slope
      real    rcoef2(ncoef,nang,nchn3),rcoef(ncoef,nang,nchn2)
      real    alt(nspot),lza,coslzen,cosd,szen

      logical bad(nspot)

      external cosd

      data cosx/1.0, 0.9, 0.8, 0.7, 0.6, 0.5/

c*  Loop on each profile
      do n=1,nspot

         if(.not. bad(n)) then
c*  Find boundary indices
            coslzen=cosd(lza(alt(n),isp(n)))
            if(coslzen .le. 0.5) coslzen=.5001
            indx1=int(10.*(1.-coslzen))+1
            indx2=indx1+1

c*  Compute limb corrected TB for HIRS channels 1-12 (daytime obs.)
            if(szen .lt. 90.) then
               do 20 i=1,nchn3
                  tblc(n,i)=xmiss
                  y(1)=0.
                  y(2)=0.
                  do j=2,npred2(indx1,i)
                     if(tb(n,chnord2(j-1,indx1,i)) .eq. xmiss) goto 20
                     y(1)=y(1)+rcoef2(j,indx1,i)*
     .                    tb(n,chnord2(j-1,indx1,i))
                  enddo
                  y(1)=y(1)+rcoef2(1,indx1,i)
                  do j=2,npred2(indx2,i)
                     if(tb(n,chnord2(j-1,indx2,i)) .eq. xmiss) goto 20
                     y(2)=y(2)+rcoef2(j,indx2,i)*
     .                    tb(n,chnord2(j-1,indx2,i))
                  enddo
                  y(2)=y(2)+rcoef2(1,indx2,i)
                  
c*  Interpolate to find limb value
                  slope=(y(2)-y(1))/(cosx(indx2)-cosx(indx1))
                  tblc(n,i)=y(2)+slope*(coslzen-cosx(indx2))
 20            continue
c*   Limb correction using all channels (night observations)
            else
               do 30 i=1,nchn2
                  tblc(n,i)=xmiss
                  y(1)=0.
                  y(2)=0.
                  do j=2,npred(indx1,i)
                     if(tb(n,chnord(j-1,indx1,i)) .eq. xmiss) goto 30
                     y(1)=y(1)+rcoef(j,indx1,i)*
     .                    tb(n,chnord(j-1,indx1,i))
                  enddo
                  y(1)=y(1)+rcoef(1,indx1,i)
                  do j=2,npred(indx2,i)
                     if(tb(n,chnord(j-1,indx2,i)) .eq. xmiss) goto 30
                     y(2)=y(2)+rcoef(j,indx2,i)*
     .                    tb(n,chnord(j-1,indx2,i))
                  enddo
                  y(2)=y(2)+rcoef(1,indx2,i)
                  
c     *  Interpolate to find limb value
                  slope=(y(2)-y(1))/(cosx(indx2)-cosx(indx1))
                  tblc(n,i)=y(2)+slope*(coslzen-cosx(indx2))
 30            continue 
            endif
         endif
      enddo

      return
      end
c
c--------------------------------------------------------------------
c
      function lza(alt,isp)
c
c*    Computes HIRS local zenith angle
c
c     INPUT
c           alt.......Altitude in km
c           isp.......HIRS scan position (1-56)
c
c     OUTPUT
c           lza.......local zenith angle in degrees
c
      implicit none

      real    alt,scan_angle,lza,sinzen,sind,asind
      integer isp

      external sind,asind

      scan_angle=-49.5+1.8*float(isp-1)
      sinzen=abs((alt+6371.)*sind(scan_angle)/6371.)
      lza=asind(sinzen)
      return
      end
c
c----------------------------------------------------------
c* Byte swapping routines (courtesy of CIRA)
c----------------------------------------------------------
c
      subroutine swap16(in1,out1)
c
c  INPUT:   in1 - INT*2
c  OUTPUT: out1 - INT*2
c
      byte in(2),out(2)
      integer*2 in1,out1,in2,out2
      equivalence (in(1),in2), (out(1),out2)
      in2 = in1
      out(1)=in(2)
      out(2)=in(1)
      out1 = out2
      return
      end
c
c-----------------------------------------------------
c
      subroutine swap32(in1,out1)
c
c  INPUT:   in1 - INT*4
c  OUTPUT: out1 - INT*4
c
      byte in(4),out(4)
      integer*4 in1,out1,in2,out2, i
      equivalence (in(1),in2), (out(1),out2)
      in2 = in1
      do i=1,4
        out(i)=in(5-i)
      enddo
      out1 = out2
      return
      end
c
c--------------------------------------------------------
c
      function asind(y)
c
c* Computes the arcsin of ratio y (output in degrees). 
c
      real x,y
      real rtod

      data rtod/57.29577951/

      x=asin(y)
      asind=x*rtod
      return
      end 
c
c------------------------------------------------------
c 
      function sind(x)
c
c* Computes the sin of angle x (arguement in degrees). 
c
      real x,sind
      real dtor

      data dtor/0.017453292/

      sind=sin(x*dtor)
      return
      end 
c
c-----------------------------------------------------
c 
      function cosd(x)
c
c* Computes the cos of angle x (arguement in degrees). 
c
      real x,cosd
      real dtor

      data dtor/0.017453292/

      cosd=cos(x*dtor)
      return
      end 
c
c------------------------------------------------------
c
      function fscale(int2,offset,sfactor)
c
c* Makes integer*2 value from real.
c
c  INPUT
c      int2..................integer*2 value to be restored to float
c      offset................offset factor
c      sfactor...............scale factor
c  OUTPUT
c      fscale................scaled value (int*2)
c---------------------------------------------------------------------
c
      implicit none

      integer*2 int2
      real fscale,offset,sfactor

      fscale=(float(int2)/sfactor)+offset

      return
      end
 
