      program BUFR_TRANHIRS2
c$$$  main program documentation block
c                .      .    .                                       .
c main program: BUFR_TRANHIRS2
C   PRGMMR: KEYSER           ORG: NP22        DATE: 2004-01-23
c
c abstract: read raw hirs-2 1b format file, decode, write selected
c   observations to output file.
c
c program history log:
c 1998-04-02  russ treadon, np2 - original author
c 1998-04-06  russ treadon, np2 - 4 digit year for obs date/time
c 1998-04-07  russ treadon, np2 - remove calls to w3fs12o
c 1999-09-23  X.Su,np2          - seperate decoding msu from hirs-2
c 2000-09-06  woollen,np2       - added second output in BUFR
c 2001-01-31  keyser   -- added error handling when no output is
c                         created so that subsequent tranjb's are
c                         skipped
c 2002-02-11  woollen  -- modifications and corrections to output BUFR
c                         dataset: "SAID" (0-01-007) corrected to
c                         proper WMO Code Table value (was 14 for
c                         NOAA-14, etc.), "SIID" (0-02-019) repl.
c                         "SIDU" (0-02-021) which didn't seem to
c                         be correct, "HMSL" (0-07-002) corrected
c                         to proper units of meters (was being stored
c                         in km), "LSQL" (0-08-012) corrected to proper
c                         WMO Code Table value (0-land/1-sea) (was
c                         backwards), "TMBR" (0-12-163) corrected to
c                         proper units of K (was being stored as
c                         K + 273.15), channel 20 "TMBR" set to missing
c                         for HIRS-2 and HIRS-3 types
c 2004-01-23  keyser   -- based on new namelist switch "compress", now
c                         has option to write compressed BUFR messages
c                         using WRITCP instead of WRITSB (removes the
c                         need for the downstream program BUFR_COMPRESS);
c                         corrected error which prevented the last
c                         message from being written to the BUFR file
c
c usage:
c
c   input files:
c     unit 05     - namelist input
c     unit 11     - binary hirs-2 1b file
c     unit lundx  - BUFR mnemonic table
c
c   output files:
c     unit 06     - printout
c     unit 51     - binary (ieee) hirs-2 brightness temperature file
c     unit lubfr  - BUFR file containing hirs-2 data
c
c   subprograms called:
c     unique:    - calday   ctime    lbyte    nsat1b   chars
c                  ichars   bufr1b
c     library:
c       w3lib    - w3tagb  w3tage  errexit w3fs26
c       bufrlib  - openbf  closbf  openmb  writsb  writcp  ufbseq
c                  mesgbc
c
c
c   Exit states
c     0 = no errors detected
c     2 = unknown surface type
c     3 = problem reading header record of 1b hirs-2 file
c     6 = unknown satellite id
c     7 = unknown satellite instrument
c   253 = no records written, all subsequent processing disabled
c
c attributes:
c   language: fortran 90
c   machine:  ibm-sp
c
c$$$
c-----------------------------------------------------------------------
c
c     Set parameters
c
c     Declare variables
      integer stdin, stdout
      character*80 infile,outfile,coefile
      character*8  compress
      namelist /input/ infile,outfile,coefile,compress

      common/switches/compress
c
c     Set i/o unit numbers
      data stdin, stdout / 5, 6 /
c
c************************************************************************
c     Start bufr_tranhirs2 here.
c
c     Call w3 routines to signal start of code
      call w3tagb('BUFR_TRANHIRS2',2004,0023,0068,'NP22')

      print *
      print *, 'WELCOME TO BUFR_TRANHIRS2 - Version 01/23/2004'
      print *
c
c

c      Get user input from namelist.  Echo to stdout
      read(stdin,input)
      write(stdout,*)'Run bufr_tranhirs2 with following user input'
      write(stdout,input)
      write(stdout,*)' '
c
c     Process 1b satellite data.
      write(stdout,*)' ' 
      write(stdout,*)'Call hirs2'
      call hirs2(infile,outfile)
c
c
      call w3tage('BUFR_TRANHIRS2')
c
c
c     End of program
      stop
      end


      real(8) function bright(r,k)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    bright      does what function bright does
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: calculates brightness temperature from radiance and 
c   hirs-2 channel number.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    bright(r,k)
c   input argument list:                                               
c     r        - real argument containing radiance.
c     k        - integer argument containing hirs-2 channel number.
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
      integer,parameter::real_32=selected_real_kind(6,37)
      integer,parameter::real_64=selected_real_kind(15,307)
c
c
c     Declare variables
      integer k
      real(real_64) r,fk1(20),fk2(20),fnu(20),tc(2,20)
      real(real_64) expn,tt
c
c     Declare common block
      common /plnctn/fnu,fk1,fk2,tc
c
c
c******************************************************************
c     Start bright here.
c
c     Compute apparent brightness temperature, tt, using
c     inverse Planck function.  This is equation 3.3.1-2
c     in NOAA POD User's Guide.  You must call PFCOEF 
c     prior to calling this routine to load the hirs-2
c     channel dependent coefficients fk1 and fk2.
c
      expn = fk1(k)/r + 1.
      tt = fk2(k)/log(expn)
c
c     Evaluate equation 4.5.1-2 of NOAA POD User's Guide.  
c     This converts the apparent brightness temperature,
c     tt, to the "true" brightness temperature.
c
c     You must call PFCOEF prior to calling this routine to
c     load the hirs-2 channel dependent coefficients tc(1,k)
c     and tc(2,k)
c
      bright = (tt-tc(1,k))/tc(2,k)
c
c
c     End of function
      return
      end
      subroutine calday(nyrday,kyear,nmon,kalday,cmon)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    calday      finds month, day from year, day-of-year
c   prgmmr: katz,b           org: np2        date: 97-11-06             
c                                                                       
c abstract: calculates month and day-of-month given year and day-of-
c   year.
c                                                                       
c program history log:                                                  
c   97-11-06  katz                                                 
c                                                                       
c usage:    call calday(nyrday,kyear,nmon,kalday,cmon)
c   input argument list:                                               
c     nyrday   - integer argument holding year and day-of-year 
c                in yyddd format.
c
c   output argument list:
c     kyear    - integer argument holding 2-digit year
c     nmon     - integer argument holding month-of-year
c     kalday   - integer argument holding day-of-month
c     cmon     - character*4 argument holding abbreviated month
c
c   subprograms called:
c     w3lib    - w3fs26
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 24.03.92
c ... scalar arguments ..
      integer kalday,kyear,nmon,nyrday
      character*4 cmon
c ...
c ... local scalars ..
      integer jday,month
c ...
c ... local arrays ..
      character*4 cmons(12)
c ...
c ... external subroutines ..
      external w3fs26
c ...
c ... intrinsic functions ..
      intrinsic mod
c ...
c ... data statements ..
c ... nyrday is in the form 'yyddd'
c ... kyear,nmon,kalday are binary integer year,month,day. cmon is
c ... the month expressed in left-justified three-letter code.
      data cmons/'jan ','feb ','mar ','apr ','may ','jun ','jul ',
     +    'aug ','sep ','oct ','nov ','dec '/
c ...
      JULIAN(IYR,IDYR) = -31739 + 1461 * (IYR + 4799) / 4
     &     -3 * ((IYR + 4899) / 100) / 4 + IDYR
c
      kyear = nyrday/1000
      jday  = mod(nyrday,1000)
      iyr = kyear
      idy = jday
      IF (IYR.GE.0.AND.IYR.LE.99) THEN
         IF (IYR.LT.21) THEN
            KYR = IYR + 2000
         ELSE
            KYR = IYR + 1900
         END IF
      ELSE
         KYR = IYR
      ENDIF
      jdn = julian(kyr,idy)
      call w3fs26(jdn,iyear,jmo,jda,idaywk,idayyr)
      month = jmo
      kalday = jda
      cmon = cmons(month)
      nmon = month
      return
      end
      subroutine chars(iword,len,cword)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    chars       converts integer to character          
c   prgmmr: katz,b           org: np2        date: 97-11-06             
c                                                                       
c abstract: turns integer into character string of specified length,
c   starting at low-order byte of integer.
c                                                                       
c program history log:                                                  
c   97-11-06  katz                                                 
c                                                                       
c usage:    call chars(iword,len,cword)
c   input argument list:                                               
c     iword    - integer argument.
c     len      - integer argument holding number of low-order bytes 
c                of first argument to convert into character.
c
c   output argument list:
c     cword    - character argument.
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
      character*1 cword(len)
      integer iword
      do ic = len , 1 , -1
        ibeg = (len - ic) * 8
        ichr = ibits(iword,ibeg,8)
        cword(ic) = char(ichr)
      enddo
      return
      end
      subroutine ctime(hour,chour)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    ctime       takes real hours and gives char hh:mm:ss    
c   prgmmr: katz,b           org: np2        date: 97-11-06             
c                                                                       
c abstract: takes floating-point hours as input and returns hh:mm:ss
c   in character form.
c                                                                       
c program history log:                                                  
c   97-11-06  katz                                                 
c                                                                       
c usage:    call ctime(hour,chour)
c   input argument list:                                               
c     hour     - real argument holding hour
c
c   output argument list:
c     chour    - character*8 argument holding 'hh:mm:ss'
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... convert real*4 hour to hh:mm:ss
c ... scalar arguments ..
      real hour
      character*8 chour
c ...
c ... local scalars ..
      integer ih,im,is
      character co
c ...
c ... intrinsic functions ..
      intrinsic float,mod,nint
c ...
c ... data statements ..
      data co/':'/
c ...
 
      ih = hour
      is = nint(3600* (hour-float(ih)))
      im = is/60
c ... trap special case, when im = 60
      if (im.eq.60) then
        ih = ih + 1
        im = 0
      end if
      is = mod(is,60)
      write (chour,fmt='(i2.2,a1,i2.2,a1,i2.2)') ih,co,im,co,is
      return
      end
      subroutine dattim(idt,ndt)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    dattim      performs time-stamp conversion           
c   prgmmr: katz,b           org: np2        date: 97-11-06             
c                                                                       
c abstract: converts year, day-of-year, and second-of-day into year,
c   month-of-year, day-of-month, hour-of-day, minute-of-hour, and
c   second-of-minute.
c                                                                       
c program history log:                                                  
c   97-11-06  katz                                                 
c                                                                       
c usage:    call dattim(idt,ndt)
c   input argument list:                                               
c     idt      - integer array argument containing three members:
c                year, day-of-year, and second-of-day.
c
c   output argument list:
c     kdt      - integer array argument containing six members:
c                year, month-of-year, day-of-month, hour-of-day, 
c                minute-of-hour, and second-of-minute.
c
c   subprograms called:
c     w3lib    - w3fs26
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 30.07.79
 
c ... array arguments ..
      integer idt(3),ndt(6)
c ...
c ... local scalars ..
      integer iday,ihr,imin,imon,isec,iyr,jday,jsec
c ...
c ... external subroutines ..
      external w3fs26
c ...
c ... intrinsic functions ..
      intrinsic mod
c ...
c
       JULIAN(IYR,IDYR) = -31739 + 1461 * (IYR + 4799) / 4
     &                    -3 * ((IYR + 4899) / 100) / 4 + IDYR
c
      iyr  = idt(3)
      jday = idt(2)
c
c     If year is two digits, convert to 4 digits
      if (iyr.ge.0.and.iyr.le.99) then
         if (iyr.lt.21) then
            kyr = iyr + 2000
         else
            kyr = iyr + 1900
         endif
      else
         kyr = iyr
      endif
c
c     Compute julian day number as number days after 4713 bc
      idy = jday
      jdn = julian(kyr,idy)
c
c
      call w3fs26(jdn,iyear,jmo,jda,idaywk,idayyr)
c
      imon = jmo
      iday = jda
      jsec = idt(1)
      ihr = jsec/3600
      imin = mod(jsec,3600)/60
      isec = jsec - 3600*ihr - 60*imin
      ndt(1) = iyr
      ndt(2) = imon
      ndt(3) = iday
      ndt(4) = ihr
      ndt(5) = imin
      ndt(6) = isec
      return
      end
cfpp$ expand(lansea)
      subroutine genval(il,ie,iheader,rbuf,ndata,rdata,nhd,nval)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    genval      creates values for output to tovs1b files 
c   prgmmr: katz,b           org: np2        date: 97-11-07             
c                                                                       
c abstract: generates values for output to tovs1b file.  calls timsat
c   to determine time within hirs-2 scans.
c                                                                       
c program history log:                                                  
c   97-11-07  katz                                                 
c   99-09-23  modified by X.Su for hirs-2 data
c                                                                       
c usage:    call genval(il,ie,iheader,inbuf,nsub,values)
c   input argument list:                                               
c     il       - integer argument containing line number.
c     ie       - integer argument containing spot number.
c     iheader   - integer array argument holding file header information.
c     inbuf    - integer array argument holding individual observation.
c     nsub     - integer argument.
c
c   output argument list:                                               
c     values   - real array argument.
c
c   subprograms called:
c     unique:  - timsat                          
c     library:
c       oneblib- calday   lansea   hrtopo
c       w3lib  - w3fs26   w3log    w3tage
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c
c     Set parameters
      integer,parameter::real_32=selected_real_kind(6,37)
      integer,parameter::real_64=selected_real_kind(15,307)
c
c
c     Declare variables.
      integer stdout
      integer iheader(nhd)
      real(real_64) rbuf(nval)
      real(real_32) rdata(ndata)
      real(real_64) rlat,rlon,xlat,xlon,rmiss
      character*4 cmon
c
c     Set unit for standard out.  Set missing data value.
      data stdout / 6 /
      data imiss,rmiss / 99999, 9999.99d0 /
c
c***********************************************************
c     Start genval here
c
c
c     Intialize variables.  Initialize output array values
c     to missing data value
      rdata = rmiss
c
c
c     Generate date/time steps.
      call calday(iheader(1),iys,ims,ids,cmon)
      ihrs = iheader(2)/10000
      imns = mod(iheader(2)/100,100)
      iscs = mod(iheader(2),100)
      call timsat(iys,ims,ids,ihrs,imns,iscs,il,ie,
     *     jys,jms,jds,jhrs,jmns,jscs)
c
c
c     Setup satellite id.  noaa11=11, noaa12=12, noaa14=14
      rdata(1) = iheader(5) + 1.d-6
c
c     Set instrument type.
         instr = 1
         nchan = 20
      rdata(2) = instr + 1.d-6

c
c     
c     Load date/time info into output array values.
c     If year is not Y2K compliant, make it so.
cy2k
      if (jys.lt.1900) then
c
c        Assume that if the year is 76-99, we mean 1976-1999
         if (jys.gt.75) then
            joff = 1900
c
c        Otherwise we assume 00-75 means 2000-2075
         else
            joff = 2000
         endif
         jys = joff + jys
      endif
cy2k
      rdata(3) = jys + 1.d-6
      rdata(4) = jms + 1.d-6
      rdata(5) = jds + 1.d-6
      rdata(6) = 3600.d0*jhrs + 60.d0*jmns + jscs + 1.d-6
c
c
c
c     Set line and scan position information.
      rdata(8) = ie + 1.d-6
      if ( (ie.lt.1) .or. (ie.gt.56) ) then
        call w3tage('BUFR_TRANHIRS2')
        call errexit(2)
      endif
c
c
c
c     Set observation (lat,lon).  If not within geographic
c     limits, we have a problem with this obs.  Write message
c     to standard out and return to calling routine.
      xlat = rbuf(1)
      xlon = rbuf(2)
      if ( (xlat.le. 90.) .and. (xlat.gt. -90.) .and.
     x     (xlon.le.180.) .and. (xlon.ge.-180.) ) then  
         rdata(9)  = xlat          
         rdata(10) = xlon           
      else
         write(stdout,*) 'funny lat/lon data:',xlat,xlon,
     x        ' skip this data point'
         return
      endif
c
c
c     Set local zenith angle (lza).  Currently (980102) this is 
c     the local zenith angle at the furthest scan position.  At
c     a later date it may be the actual scan position for the
c     given scan position.
c
c     Array element 42 in rbuf contains the lza.  Load into 
c     output array element 3.  Units are degree.
c
      rdata(11) = rbuf(42)
c
c     Future code to compute "true" local zenith angle
cx      dg2rad = acos(-1.)/180.
cx      erad   = 6371.2
cx      hgt    = rbuf(41)
cx      rlza   = rbuf(42)
cx      x      = sin(dg2rad*rlza)
cx      ratio  = 1. + hgt/erad
cx      y      = x/ratio
cx      if (y.gt. 1.) write(6,*)'WARNING y>+1 ',hgt,rlza,x,y
cx      if (y.lt.-1.) write(6,*)'WARNING y<-1 ',hgt,rlza,x,y
cx      y      = max(-1.,y)
cx      y      = min(y,1.)
cx      zref   = -1.*abs(asin(y)) / dg2rad
c
c     Set step size (degrees) for hirs-2
cx      anginc = 1.80
c
c     Variable ie contains the scan step position
c
cx      truzla     = zref + (ie-1)*anginc
cx      rvalues(3) = truzla
c
c
c
c     Set solar zenith angle if within physical bounds.
      if ( (rbuf(3).le.90.) .and. 
     x     (rbuf(3).ge.0.) ) then
         rdata(12) = rbuf(3)
      endif
c
c
c     Set surface type information based on resolution option.
c     use low resolution data to set surface type.
c     
         ils = lansea(xlat,xlon,ll)
c
c        Sea point
         if (ils.eq.2) then     
            rdata(7)  = 0 + 1.d-6
            rdata(13) = 0.0
c
c        Land point
         else if (ils.eq.1) then
            rdata(7)  = 1 + 1.d-6
            rdata(13) = 1.0*ll
c
c        Unknown surface type.
         else
            write(stdout,*) '** unknown surf type at:',
     x           xlat,xlon,ll,ils
            call w3tage('BUFR_TRANHIRS2')
            call errexit(2)
         endif
c
c
c     Set satellite height (km).  Input array rbuf contains 
c     this information in array element 41.  Load into output
c     array element 6.  At a later date, this may be replaced
c     with the temperature of the satellite sensor.
c
      rdata(14) = rbuf(41)
c
c
c     Set equivalent black body brightness temperatures
         do ichan = 1,nchan-1  
            rdata(14+ichan) = rbuf(3+ichan)
         end do    
         rdata(14+20) = rbuf(23)
c
c
c     If temperatures exceed 400 in absolute values, set the
c     temperature to the missing value.  
c
      do ichan = 1,nchan
         if (abs(rdata(14+ichan)).gt.400.) then
            rdata(14+ichan) = rmiss
         endif 
      end do
c
c
c     End of routine.
      return
      end
cfpp$ expand(bright,rbiref,lchar,mbyt16,xfloat,lbit,lbyt16,lbyte,lngtst)
cfpp$ expand(rsolzn)
      subroutine hirs2(infile,outfile)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram    hirinw      performs hirs-2 processing in tovinw
c   prgmmr: Keyser         org: np22       date: 2004-01-23
c                                                                       
c abstract: performs hirs-2 processing, including error checking and
c   calibration.
c                                                                       
c program history log:                                                  
c   97-11-07  katz                                                 
c   99-09-23  X.Su seperate decoding msu from hirs-2
c 2004-01-23  keyser   -- corrected error which prevented the last
c                         message from being written to the BUFR file
c                                                                       
c usage:    call hirinw(iu,nu,idsat)
c   input argument list:                                               
c     iu       - integer argument containing input hirs-2 unit number.
c     nu       - integer argument containing output hirs-2 unit number.
c     idsat    - integer argument containing satellite id.
c
c   input files:
c     unit iu  - binary file of hirs-2 data
c
c   output files:
c     unit 06  - printout
c     unit 35  - binary file of correction coefficients
c     unit nu  - binary file of calibrated and checked hirs-2 data
c
c   subprograms called:
c     unique:    arcsin   sine (statement functions)
c     oneblib  - bright   xfloat   lbit     lbyt16   lbyte    lchar
c                lngtst   mbyt16   nsat1b   calday   pfcoef
c                ichars
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c
c
c     Include machine dependent parameters.
      include 'rfac.inc'
      
c     Declare/set parameters.
c      nbi  = total number of bytes (4253) in HIRS-2 data record 
c      leno = number of bytes (3844) in HIRS-2 data record with info
c      nwo  = number of 4-byte words in leno bytes (3844/4=961)
c      nbx  = closest multiple of 4 greater than nbi (4253/4=1063.25,
c              so nbx = 1064*4 = 4256 = 4253+3 = nbi+3)
c      nwi  = number of 4-byte words in nbx bytes (4256/4=1064)
c
c      NOTE:  Byte counts for HIRS-2 data record come from Table
c             4.1.2.1-1 of NOAA POD User Guide (Aug '97 version).  For
c             more info point your web browser to 
c             -->  http://www2.ncdc.noaa.gov/POD/podug/index.hmtl
c
      integer,parameter::real_32=selected_real_kind(6,37)
      integer,parameter::real_64=selected_real_kind(15,307)
      real(real_64) eps
      integer nbi,nbx,nwi,nwo,leno
      parameter (nbi=4253,nbx=nbi+3,nwi=nbx/4,nwo=961,leno=nwo*4)
c
      integer mrow,mcol,lenb,nsub,lens,nid,nnd,nhd
      parameter (mrow=1100,mcol=56)
      parameter (lenb=112,nsub=4,lens=lenb/nsub)
      parameter (nc=19,np=30,nid=3,nnd=6)
      parameter (nval=42)
      parameter (eps=1.d-3,mch=20)
c
c     Set parameter for structure of output data file.
      parameter (nreal=14,nchan=20)
      parameter (ndata=nreal+nchan)
      parameter (nhd=30)
c
c
c     Declare variables.
      logical badhgt,irok,hirsok,locok
c
      character*1 kbuf(nbx)
      character*4 cmon,indat(nwo),jbuf(nwi),cins(0:1),cflag
      character*8 cunk,time,csat(0:14),subset
      character*12 date
      character*80 infile,outfile
c
      integer idsat,iu,nu,linmax,stdout
      integer init,ioff,ispot,jline,k,locw,misf,nbadl,nbads,nbadx
      integer i,initms,initmt,irec,iyyddd,j,jb,jbyte,jch,jddd,jdispl
      integer jhour,jmin,jnkodr,joff,jpos,jsat,jsec,jstart,jstop,jtype
      integer jyd,jyear,jyyddd,ksat,lerrs,linco,line,lkalib
      integer lnode,lset,ltflag,mfn,mtdif,mtime,mtlast,nbad,nch,nd,nm
      integer nopos,nprs,nrec,nscans,nspac,ny
      integer newpos(20),idt(nid),ndt(nnd)
      integer ihdrhir(nhd)
      integer ichan(nchan)
c
      real(real_64) height,zangle,hgtsav,hgtnom,sinzen,r,rtime
      real(real_64) ang,counts,rads,sctime,spacnt,spavg
      real(real_64) test,viscon,x,cosz,expn,tb,f1,f2,arg,ta
      real(real_64) rtlast,rtdif,rlat,rlon,rzen,rvis,rszn,rbdr
      real(real_64) rtc4,rmisf
      real(real_64) fk1(20),fk2(20),fnu(20),tc(2,20)
      real(real_64) c0(20),c1(20),c2(20),cvis(14),r0(20),r1(20),r2(20)
      real(real_64) reason(20,2),sathgt(2),ylats(56),ylons(56),zens(56)
      real(real_64) values(nval),rbuf(np)
      real(real_64) rbufout(nval),rtem(np),tbright(mch),rtsol(3)
      real(real_32) rdata(ndata)
c
      double precision two22,two30,two44
c
c     Declare equivalences
      equivalence (kbuf(1),jbuf(1))
c
c     Declare common block
      common /plnctn/fnu,fk1,fk2,tc
c    
c     Declare statement functions.
      real(real_64) arcsin,sine
c    
c     Data statements.
      data (cvis(i),i=1,4)/4*0.d0/
      data cvis(5)/2.84547d0/
      data cvis(6)/7.6975d0/
      data cvis(7)/2.84547d0/
      data cvis(8)/4.696d0/
      data cvis(9)/4.56359d0/
      data (cvis(i),i=10,14)/5*4.6d0/
c
      data sathgt/860.d0,820.d0/,spavg/-3617.d0/
c
      data newpos/1,17,2,3,13,4,18,11,19,7,8,20,10,14,6,5,15,12,16,9/
c
      data two22/4194304.d+0/
      data two30/1.0737418d+9/
      data two44/1.7592186d+13/
c
      data c0/20*0.d0/,c1/20*0.d0/,c2/20*0.d0/
      data r0/20*0.d0/,r1/20*0.d0/,r2/20*0.d0/
c
      data (reason(i,1),i=1,12)/7*15.d0,6.d0,3.5d0,1.d0,.5d0,.25d0/
      data (reason(i,1),i=13,20)/.003d0,3*.002d0,.001d0,.0004d0,
     x     .0001d0,0.d0/
      data (reason(i,2),i=1,12)/5*120.d0,130.d0,145.d0,245.d0,
     x     100.d0,175.d0,2*50.d0/
      data (reason(i,2),i=13,20)/7*10.d0,0.d0/
c    
c     More data statements.
      data cins/'unkn','hirs'/
      data csat/5*' ','tiros-n','noaa-6','noaa-7','noaa-8','noaa-9',
     +     'noaa-10','noaa-11','noaa-12','noaa-13','noaa-14'/
      data cunk/'unknown'/
      data rmisf/9999.99d0/
      data nite / 90.00 /
c
c     Set flag for first write to output file.  Set
c     parameter for header record.
      data ichan / 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
     x     12, 13, 14, 15, 16, 17, 18, 19, 20 /
c
c     Data statements for standard out and hirs-2 1b file.
      data stdout /  6 /
      data lundat, lundx, lunout, lubfr / 11, 12, 51, 52 / 
c    
c    
c     Define statement functions.
      sine(x) = sin(.0174533d0*x)
      arcsin(x) = 57.29578d0*atan2(x,sqrt(1.-x*x))
c
c
c*********************************************************************
c     Start hirs-2 here
c 
      write(stdout,*)' '
      write(stdout,*)' begin hirs-2  ingest'
c
c     Initialize variables
      csat(0) = cunk
c
c
c     Open the BUFR  output file.
      call openbf(lubfr,'OUT',lundx)
      subset = 'NC021021'
 
c     open and Write header record to binary output file
      open(lunout,file=outfile,form='unformatted',status='unknown')
      write(lunout) nreal,nchan,(ichan(i),i=1,nchan)
cx      write(99,*) nreal,nchan
cx      write(99,*) (ichan(i),i=1,nchan)
c
c
c     Open unit to raw 1b hirs-2 data file.  Read header record.
c     See if valid satellite.  If not, exit routine.
c
      open(lundat,file=infile,recl=nbi/rfac,
     &      access='direct',status='old',err=1950)
      nri = 1
      read (lundat,rec=nri,err=1900) (kbuf(i),i=1,nbi)
      jsat = mova2i(kbuf(1))
      if (jsat.eq.0) then
        write(stdout,*)' unable to extract satellite id from header.'
        goto 2100
      end if
c
c
c     Load header record into work array.
c
c     A 1b HIRS-2 data record nbi=4253 bytes long.  The closest multiple
c     of four greated than 4253 is nbx=4256.  4256 bytes correspond to
c     nwi=1064 four-byte words.  
c
c     Only the first 3844 bytes of the HIRS-2 1b data record contain data.
c     leno=3844 bytes corresponds to nwo=961 4-byte words.
c
c     We read the 1b HIRS-2 data into elements 1 through nbi of character*1
c     array KBUF.  This array is equivalenced to character*4 array JBUF 
c     (length nwi).  Load data records into working array INDAT.  Since 
c     bytes 3845 through 4253 are not used, only load the first nwo=961
c     character*4 words (or 3844 bytes) of JBUF into INDAT.
c
      do i = 1,nwo
         indat(i) = jbuf(i)
      end do
c
c
c     Unpack date/time/satellite id info from data record.
c
      jyear = lbyte(17,7,indat)
c
cy2k
c     If year is two digit, convert to 4 digit.
      iyr = jyear
      if (iyr.ge.0.and.iyr.le.99) then
         if (iyr.lt.21) then
            kyr = iyr + 2000
         else
            kyr = iyr + 1900
         endif
      else
         kyr = iyr
      endif
      jyear = kyr

      idsat = nsat1b(jsat,jyear)
      if (idsat.eq.0) then
        write(stdout,*)' unrecognized satellite ident ',idsat
        call w3tage('BUFR_TRANHIRS2')
        call errexit(3)
        stop ' *** unrecognized satellite ident *** '
      else
        write(stdout,*)' data source is satellite number ',idsat
      end if
c
      jtype = lchar(2,indat)
      if (jtype.eq.81) goto 110
      if (jtype.eq.82) goto 110
      write(stdout,*)' *** unrecognized data type, jtype=',jtype
      call w3tage('BUFR_TRANHIRS2')
      call errexit(3)
      stop ' *** unrecognized data type *** '
c
 110  continue
      inst = (jtype-81)/16 + 1
      if (inst.ne.1) inst = 0
      if (inst.eq.0 .or. idsat.eq.1) goto 2100
      lines = lbyte(65,16,indat)
c
      write(stdout,1100) jtype,cins(inst),jsat,csat(idsat),lines
 1100 format (' data type = ',i2,' => ',a4,', s/c id = ',i1,' => ',a8,
     +       '... no.lines = ',i4)
c
      jstart = .001*float(lbyte(38,27,indat)) + .5
      jddd   = lbyte(24,9,indat)
      jyd    = 1000*jyear + jddd
      call calday(jyd,ny,nm,nd,cmon)

      jstop  = .001*float(lbyte(102,27,indat)) + .5
      jhour  = jstart/3600
      jmin   = (jstart-3600*jhour)/60
      jsec   = jstart - 3600*jhour - 60*jmin
      write(stdout,120) jhour,jmin,jsec,nd,cmon,ny,lines
  120 format (' hirs-2 starting time ',3i3,2x,i3,1x,a4,i4,i6,
     +       ' scan lines')

      jhour  = jstop/3600
      jmin   = (jstop-3600*jhour)/60
      jsec   = jstop - 3600*jhour - 60*jmin
      write(stdout,130) jhour,jmin,jsec
  130 format (7x,' ending time ',3i3,/)


c     Load necessary information for output file header.


        idt(1)=.001*float(lbyte(38,27,indat))+0.5/64.0
        idt(2)=mod(jyd,1000)
        idt(3)=jyd/1000


        call dattim(idt,ndt)

        ihdrhir     =0
        ihdrhir(1) =jyd
        ihdrhir(2) =ndt(4)*10000+ndt(5)*100+ndt(6)
        ihdrhir(4) =lines
        ihdrhir(5) =idsat

c***********************************************************************
c   Initialize common blocks which contain vital coefficients used in  *
c   later routines.  The call to pfcoef reads in data from the user    *
c   spcified coefficient data file                                     *
c***********************************************************************

        write (stdout,*) '  '
        call pfcoef(idsat)
        write(stdout,*) 'HIR2: call pfcoef for idsat=',idsat


c
c
c     Prepatory initializations prior to reading in satellite data.
c
c       mtlast:  large value ensures program will not make a missing
c                data call on the first s.e.b.
c       lersum:  logical sum of all errors
c       nopos :  number of undefined encoder positions
c       nbad  :  number of bad altitudes
c       spacnt:  space count accumulator for channel 20
c       nspac :  number of space looks in channel 20
c        
      ksat   = 2 - mod(idsat,2)
      hgtnom = sathgt(ksat)
      hgtsav = 0.
      viscon = cvis(idsat)
      rtlast = -6400000.
      nopos  = 0
      spacnt = 0.
      nspac  = 0
      nprblm = 0
      nbadl = 0
      nbadh = 0
      nbadk = 0
      ncldy = 0
      nrec  = 0
      nrep  = 0
c
c
c
c     Main loop over number of scans
c
      nlo = 0
 1200 continue
c
c        Increment record counter.  Read data record.  Load 
c        into local work array.  Increment line counters.
         nri = nri + 1
         read(lundat,rec=nri,err=1600) (kbuf(i),i=1,nbi)
         do i = 1,nwo
            indat(i) = jbuf(i)
         end do
         nlo  = nlo + 1
         line = nlo

c
c
c        Unpack satellite height (km) and local zenith angle (128ths
c        of a degree).  Convert angle to whole degrees.
c
c        According to Table 4.1.2.1-1 of the NOAA POD Users Guide
c        bytes 737-738 contain the height (or altitude)
c        bytes 739-340 contain the local zenith angle (lza)
c        
         height = lbyt16(369,indat)
         zangle = lbyt16(370,indat)/128.d0
c
c        Check for bad satellite height.  If bad value found, set
c        to last good or norminal value following logic below.
         badhgt = .false.
         if (height.lt.760.d0 .or. height.gt.920.d0) then
            badhgt = .true.
            if (nbadh.lt.20) 
     x           write(stdout,*)' bad satellite height ',height,
     x           ' in record ',nri
            nbadh = nbadh + 1
c
c           Replace bad height with last good value.  If there
c           is no good last value, replace bad value with nominal
c           value.
            if (hgtsav.ne.0.) then
               height = hgtsav
            else
               height = hgtnom
            end if
         endif
         if (.not.badhgt) hgtsav = height

c      Extract time information
c
            msec = lbyte(38,27,indat)
            hour = float(msec)/3600000.
            jy   = lbyte(17,7,indat)
c
cy2k
c           If year is two digit, convert to 4 digit.
            if (jy.ge.0.and.jy.le.99) then
               if (jy.lt.21) then
                  jy = jy + 2000
               else
                  jy = jy + 1900
               endif
            endif

            jd   = lbyte(24,9,indat)
            iyd  = 1000*jy + jd
            call calday(iyd,iy,im,id,cmon)
            call ctime(hour,time)

c        Periodically echo date/time stuff to standard out
         if (mod(nlo,250).eq.1) then
            write (date,1400) id,cmon,iy
 1400       format (i2.2,1x,a4,i4)
            write(stdout,*)' line ',nlo,' ... ',date,'  ',
     x           time,'  ',height,zangle
         end if
c
c        Test for 'fatal' and 'synch' error conditions.
         lset = 2*lbit(65,indat) + lbit(82,indat)
c
c        Extract line/time info.  (sctime = time of day
c        of start of scan in seconds).  Save date/time
c        info for first record.
         linco  = lbyt16(1,indat)
c

         jyyddd =iyd 
         mtime  =msec 
         sctime = 1.d-3*float(msec)
c
c
c        Loop over 56 scan steps.  For each, pull out the 
c        earth location (lat,lon) data.  Units for lat,lon
c        in 1b file is 128ths of a degree.  Thus, divide by
c        128 to get whole degrees.
c
c        According to Table 4.1.2.1-1 the earth location 
c        (lat,lon) data for the scan steps are stored in
c        two-byte fields (128ths of a degree) starting at
c        byte 741.
c
c        From the loop below we make the following correspondence
c
c         j   jb   jb+1   ylats    ylons
c         1   371  372   741-742  743-744
c         2   373  374   745-746  747-748
c         .   ...
c        56   481  482   962-963  964-965
c
c        
         do j = 1,56
            jb = 371 + 2* (j-1)
            ylats(j) = float(mbyt16(jb,indat))/128.d0
            ylons(j) = float(mbyt16(jb+1,indat))/128.d0
            ang = -49.5d0 + 1.8d0*float(j-1)
            sinzen = abs((height+6371.d0)*sine(ang)/6371.d0)
            zens(j) = arcsin(sinzen)
         end do
c
c
c        Set flag LNODE for ascending or descening satellite pass.
cx         lnode = 1
cx         if (lngtst(ylons(1),ylons(28),ylons(56)).eq.0) lnode = -1
c
c
c        Loop over 20 channels.  Pull out calibration coefficients.
c
c        The following information is from Section 4.1.2.1 of the 
c        NOAA POD Users Guide.  Table 4.1.2.1-1 indicates that the
c        calibration coefficients are in bytes 17 through 736.

c        The calibration coefficients are stored in three groups 
c        of 240 bytes each.  The first group contains manually 
c        calibrated coefficients for each of the 20 channels.  This
c        routine does not use this information.  
c
c        The next two groups contain auto-calibrated and normalization
c        coefficients, respectively.  For each group there are three
c        coefficients for each channel.
c        and 
c
c        From the loops below we make the following connections
c
c                      byte locations
c        nch   jb  jch    c0      c1      c3       r2      r1      r0
c         1   249   1  497-500 501-504 505-508  257-260 261-264 265-268
c         2   255  17  509-512 513-516 517-520  269-272 273-276 277-280
c         .   ...  ..
c        20   363   9  725-728 729-732 733-736  485-488 489-492 493-496
c
c
         do nch = 1,20
            jb = 249 + 6* (nch-1)
c
c           The counts are not in order of channel number. The
c           order is given in array NEWPOS (comment by Lauritsen,
c           3 Sep '79).  Set JCH to be the channel number.
            jch = newpos(nch)
c
            if (jch.eq.20) goto 150
c
c           Get auto-calibrated coefficients.  Divide each by the 
c           appropirate scale factor.
            c0(jch) = xfloat(jb,indat)/two22
            c1(jch) = xfloat(jb+2,indat)/two30
            c2(jch) = xfloat(jb+4,indat)/two44
c
c           Unpack the normalization coefficients used in the
c           auto-calibration.  Note that the normalization coefficients
c           are given in the reverse order of the auto-calibration
c           coefficients.  Apply the scaling factor accordingly.
            jb = jb - 120
            r2(jch) = xfloat(jb,indat)/two44
            r1(jch) = xfloat(jb+2,indat)/two30
            r0(jch) = xfloat(jb+4,indat)/two22
c
c           For channels 3 through 20, skip the following tests.
            if (jch.gt.2)    goto 150
c
c           Skip the following tests for noaa satellites 5,9,11
            if (idsat.eq.5)  goto 150
            if (idsat.eq.9)  goto 150
cold        if (idsat.eq.11) goto 150   noaa-11 has same problem as 14!
c
c
c           Specific tests for channel 2.
            if (jch.eq.2) then
               if (idsat.ne.12) goto 150
               test = abs(r0(2))
               if (test.ge.200.d0) goto 150
               test = test + 512.d0
               r0(2) = sign(test,r0(2))
               goto 150
            end if
            test = abs(r0(1))
c
c           Test for channels 1,2 for satellites other than noaa-12
            if (idsat.ne.12) then
               if (test.ge.200.d0) goto 150
               test = test + 512.d0
               r0(1) = sign(test,r0(1))
               goto 150
            end if
c
            if (test.lt.200.d0) then
               test = test + 2048.d0
            else
               test = test + 1536.d0
            end if
            r0(1) = sign(test,r0(1))
c
 150        continue
c
c        End of loop unpacking calibration and normalization
c        coefficients.
         end do
c
c
c        Loop over the 56 scan steps.  
         do jpos = 1,56
c
c           Initialize variables.
            lerrs   = lset
            mfn     = jpos - 1
            tbright = rmisf
            rtem    = rmisf
c
c           Unpack the encoder position.  Possible values for
c           the encoder position are
c                jnkodr   meaning
c                1-56     earth view in scan position jnkodr 
c                68       space view
c                105      cold black body view
c                156      warm black body view
c           
c
            jbyte  = 965 + 44* (jpos-1)
            jnkodr = lchar(jbyte,indat)
c
c           Encoder position should equal scan step position.  This
c           condition does not apply for the first four scan steps.
c           If encoder position >= 68, we either have bad data or 
c           are not looking at the earth.
c
            if (jnkodr.ge.68) goto 180
            if (mfn.ge.4 .or. jnkodr.eq.jpos) goto 190
            nprblm = nprblm + 1
c
 160        lerrs  = ior(lerrs,32)
c
c           Encoder position 68 is a space view.
            if (jnkodr.eq.68) goto 170
c
c           We have a cold or warm black body encoder position
c           (or possibly bad data).  
            if (nspac.eq.0) goto 220
            spavg  = spacnt/float(nspac)
            nspac  = 0
            spacnt = 0.
            goto 220
c
c           Process data for space view.  Accumulate channel 
c           20 counts from space.
 170        jb = 22*mfn + 496
c
c           Channel 20 is in 12th position.
            spacnt = spacnt + float(mbyt16(jb,indat))
            nspac = nspac + 1
            goto 220
c
c           No data.  Goto next scan step.
 180        if (jnkodr.eq.68 .or. jnkodr.eq.105 .or.
     x           jnkodr.eq.156) goto 160
            nopos = nopos + 1
            goto 250
c
c           If we reach this point, we have a normal earth view.
c
 190        continue
c
c           The data for each of the 56 encoder positions takes 44
c           bytes or half-words.  This information is preceded by 
c           484 half-words.  
c
            jdispl = 22*mfn + 484
c
c           Loop over channels.  Convert linearized counts to radiance.
            do nch = 1,20
               jb  = jdispl + nch
               jch = newpos(nch)
               if (jch.eq.20) goto 210
               kounts = mbyt16(jb,indat)
               counts = c0(jch) + kounts* (c1(jch)+kounts*c2(jch))
               rads   = r0(jch) + counts* (r1(jch)+counts*r2(jch))
c
               if ( (rads.ge.reason(jch,1) .and.
     x              rads.le.reason(jch,2)) ) goto 200
c
               nbadk = nbadk + 1
               if (nbadk.gt.20) goto 210
               write(stdout,*)' bad counts and/or radiance in chan ',
     x              jch,counts,rads
               goto 210
c
c              Convert radiance to brightness temperature.
 200           continue
               tbright(jch) = bright(rads,jch)
c
 210           continue
            end do
c
c
c           Now that we have data unpacked from 1b data file,
c           save the data we want for further processing.
 220        continue
c
            iloc  = 65536*linco + mfn + 1
            idate = jyyddd
            rtime = 64.d0 * (sctime + 0.1d0*float(mfn))
c
c           Test for missing data preceding this scan position.
c           The time unit used in these tests is 64ths of a second.
            ltflag = 0
            rtdif  = rtime - rtlast
            if (jpos.eq.1 .and. rtdif.gt.95) ltflag = 128
            if (jpos.gt.1 .and. rtdif.gt.9)  ltflag = 128
            lerrs  = ior(lerrs,ltflag)
            rtlast = rtime
c
            rlat = ylats(jpos)
            rlon = ylons(jpos)
            rzen = zens(jpos)
c
            jb = 22*mfn + 496
c
            rvis = 0.0
            if (jnkodr.le.56 .and. abs(spavg).gt.0.) then
               rvis = viscon*
     x              (float(mbyt16(jb,indat))-spavg)
            endif
c
c
cnew
c     Add hirprw4 and soundw2 processing here.
c
c
c           Extract parameters.  Load data values into
c           array item.  Check for bad (missing) data
c           values.  Set bad data flag accordingly.
            irbad = 0
            do i = 1,nc
               rtem(i) = tbright(i)
               if ( (abs(rtem(i)).lt.eps) .or.
     x              (abs(rtem(i)-rmisf).lt.eps) ) 
     x              irbad = irbad + 1
            end do
            irok     = irbad .eq. 0
c         
c           Load data values into output array.
            rtem(20) = rvis
            rlat     = rlat
            rlon     = rlon
            rzen     = abs(rzen)
            iflag    = lerrs
c
c           Are the (lat,lon) for this location ok?  
c           If the location is bad, skip to next spot.
            locok    = (abs(rlat).le.90.00) 
     x           .and. (abs(rlon).le.180.00)
            if (.not.locok) then
              write(stdout,*)'goto 161.  rlat,rlon,locok=',
     x         rlat,rlon,locok
               nbadl = nbadl + 1
               goto 161
            endif
c
c           Set date/time parameters.
c                idt(1) = seconds of day,
c                idt(2) = day of year,
c                idt(3) = year
c
            idt(1) = rtime/64.d0
            idt(2) = mod(idate,1000)
            idt(3) = idate/1000
c     
c           Check the ingest quality flag for other fatal
c           (bad data) conditions.
c
            hirsok = .true.
            if (iflag.eq.0) goto 211
            call chars(iflag,4,cflag)
            lflag  = 0
            lflag  = lflag + lbit(25,cflag)
            lflag  = lflag + lbit(27,cflag)
            lflag  = lflag + lbit(32,cflag)
            hirsok = lflag .eq. 0
 211        continue
c
c           Process data and store in staging buffer.
c           If bad data detected, reset channels to missing.
            if (hirsok .and. irok) goto 231
            do i = 1,nc
               rtem(i) = rmisf
            end do
 231        continue
            rszn = rmisf
            rbdr = rmisf
            rtc4 = rmisf
c
c           Compute solar zenith angle.
            rszn = rsolzn(rlat,rlon,idt,cosz)

c           If is it night time skip over corrections below.
c           If data is not ok, skip over corrections.
            if (abs(rszn-nite).lt.eps) goto 241
            if (.not.hirsok) goto 241
c
c           Compute bidirectional reflectance for visible channel
c           radiance.  Use this to correction 2511 channel for
c           reflected sunlight (after limb correction).
            call dattim(idt,ndt)
            rbdr = rbiref(rtem(20),ndt(2),cosz)
c
            if (.not.irok) goto 241
c
c           LG 31-mar-95 note the next statement was added to ensure
c               consistent behaviour - kfluor was not set previously
            kfluor = 0
            if (kfluor.eq.0) goto 241
c
c           For satellites 1-10 and 12, apply correction for 2360
c           channel (tirosn - noaa10) for fluorescence (before 
c           limb correction).
            if (idsat.lt.11 .or. idsat.eq.12) 
     x           rtem(17) = rfluor(rtem(17),cosz)
c
 241        continue
c
c           Load data into output buffer.
            rbuf    = rmisf
            rbuf(1) = rlat
            rbuf(2) = rlon
            rbuf(3) = rszn
            do k = 1,nc
               rbuf(k+3) = rtem(k)
               if (abs(rtem(k)-rmisf).lt.eps) rbuf(k+3) = rmisf
            end do
            rbuf(23) = rbdr
            rbuf(25) = rzen
c
c           The following check is from the HIRS-2 limb correction
c           routine.  I was executed regardless of the limb
c           correction on/off switch.
            if (idsat.lt.11 .or. idsat.eq.12) then
               rtem1 = 0.5d0*(rtem(2)+rtem(17))
               if (abs(rtem(1)-rtem1).gt.5.) then
                  rtem(1) = rtem1
                  rbuf(4) = rtem(1)
               endif
            endif
c
c           For daytime data, compute albedo.
            if (abs(rszn-nite).gt.eps) then
               rtsol(1) = rtem(8)
               rtsol(2) = rtem(18)
               rtsol(3) = rtem(19)
               rtc4 = rrefls(rtsol)
               rbuf(24) = rtc4
            endif
c
c           We have a good spot.  Load output array.
            rbufout = rmisf
            do ii = 1,30
               rbufout(ii) = rbuf(ii)
            end do
            rbufout(41) = height
            rbufout(42) = zangle
c
c           Prepare hirs-2 data for output.
            ie = jpos
            il = line
            nrec = nrec + 1
            call genval(il,ie,ihdrhir,rbufout,ndata,rdata,nhd,nval)
 
c           Write hirs-2 for current spot to binary output file.
            write(lunout) (rdata(i),i=1,ndata)

c           Write hirs-2 for current spot to BUFR output file.
            call bufr1b(lubfr,subset,nreal,nchan,rdata,nrep)

cx            write(99,*)'nrec =',nrec
cx            write(99,1234) (rdata(i),i=1,ndata)
cx            write(99,*)' '
cx1234       format(5(g12.6,1x))
c
           if (nrec.eq.1) then
              write(stdout,*)'header information below'
              write(stdout,*)'nreal,nchan=',nreal,nchan
              write(stdout,*)'channel numbers below'
              write(stdout,*) (ichan(i),i=1,nchan)
              write(stdout,*)' data for record nrec=',nrec
              write(stdout,1010) (rdata(i),i=1,ndata)
 1010         format(8(g12.6,1x))
           endif
c
c
 161        continue
c
c
c     End of hirprw and soundw2 insert
c
c
c
 250        continue
c
c        End of loop over scan steps
         end do
c
c
c     End of loop reading in raw 1b data.
      if (nlo.lt.lines) goto 1200
cx      if (nlo.lt.3) goto 1200
c
c
c     Done reading from 1b files.  Close unit.
 1600 continue
      write(stdout,*)'hirs-2 ingest stats'
      write(stdout,*)' no. lines read       = ',nlo,lines
      write(stdout,*)' no. bad heights      = ',nbadh
      write(stdout,*)' no. bad counts       = ',nbadk
      write(stdout,*)' no. bad locations    = ',nbadl
      write(stdout,*)' no. encoder problem  = ',nprblm
      write(stdout,*)' no. bad encoder pos  = ',nopos
      write(stdout,*)' no. space views      = ',nspac
      write(stdout,*)' no. cloudy columns   = ',ncldy
      write(stdout,*)' no. records written  = ',nrec
      write(stdout,*)' no. BUFR rpts written= ',nrep
      goto 2000
c
c
c     Error reading 1b file.
 1900 write(stdout,*)' *** error reading hdr record of file ',
     x     infile
      close(lundat)
      close(lunout)
      call closbf(lubfr)
      call w3tage('BUFR_TRANHIRS2')
      call errexit(3)
c
 1950 write(stdout,*)' *** error opening file=',infile
      close(lundat)
      close(lunout)
      call closbf(lubfr)
c
c
c
c     End of HIRS-2 1b ingest.  Write some stats to standard out.
 2000 continue
      write(stdout,*)' hirs-2 1b ingest completed'
      write(stdout,*)' '
c
c
 2100 continue
      close(lundat)
      close(lunout)
      call closbf(lubfr)

      if(nrec.eq.0) then
        write(stdout,1003)
 1003   format(/' NO RECORDS WRITTEN -- DISABLING ALL SUBSEQUENT ',
     1         'PROCESSING.'/)
        call w3tage('BUFR_TRANHIRS2')
        call errexit(253)
      endif

      call mesgbc(lubfr,msgt,icomp)
      if(icomp.eq.1) then
         print'(/"OUTPUT BUFR FILE MESSAGES   C O M P R E S S E D"/'//
     .    '"FIRST MESSAGE TYPE FOUND IS",I5/)', msgt
      elseif(icomp.eq.0) then
         print'(/"OUTPUT BUFR FILE MESSAGES   '//
     .    'U N C O M P R E S S E D"/"FIRST MESSAGE TYPE FOUND IS",I5/)',
     .    msgt
      elseif(icomp.eq.-1) then
         print'(//"ERROR READING OUTPUT BUFR FILE - MESSAGE '//
     .        'COMPRESSION UNKNOWN"/)'
      elseif(icomp.eq.-3) then
         print'(/"OUTPUT BUFR FILE DOES NOT EXIST"/)'
      elseif(icomp.eq.-2) then
         print'(/"OUTPUT BUFR FILE HAS NO DATA MESSAGES"/"FIRST '//
     .    'MESSAGE TYPE FOUND IS",I5/)', msgt
      endif
      close(lubfr)

c
c
c     End of routine
      return
      end
cfpp$ expand(ichars,lbit)
      subroutine hrtopo(flat,flon,iele,isfc)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    hrtopo      high-resolution elevation, land/sea status
c   prgmmr: katz,b           org: np2        date: 97-11-06             
c                                                                       
c abstract: using high-resolution topography data set, takes latitude
c   and latitude and returns topography and land/sea condition.
c                                                                       
c program history log:                                                  
c   97-11-06  katz                                                 
c                                                                       
c usage:    call hrtopo(flat,flon,iele,isfc,*)
c   input argument list:                                               
c     flat     - real argument containing latitude.
c     flon     - real argument containing longitude.
c
c   output argument list:
c     iele     - integer argument containing surface elevation.
c     isfc     - integer argument containing land/sea condition.
c
c   input files:
c     unit 42  - binary high-resolution topography file
c
c   subprograms called:
c     oneblib  - lbit   
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    

c ... version of 03.05.95
c ... input = lat (0-89,+n), lon (0-180,+e,-w)
c ... output = elevation (meters) + sfc type (0=water, 1=land)

      include 'rfac.inc'
      
c     .. scalar arguments ..
      real flat,flon
      integer iele,isfc
c     ..
c     .. local scalars ..
      real fele
      integer i,ibit,ilat,ilon,irec,iu,iwrd,jlat,jlon,lenb,lenh,lenm,
     +        lenr,lrec,mb,mlat,mwrd,nrec,numb,nwrd,iok
      character*12 bfile,hfile
c     ..
c     .. local arrays ..
      character*4 mbit(70),mbuf(68,1080)
      integer numbs(1080)
      character*1 lbuf(4320)
c     ..
c     .. external functions ..
      integer lbit
      external lbit
c     ..
c     .. intrinsic functions ..
      intrinsic float,iabs,nint
c     ..
c     .. data statements ..
      data bfile/'mapbitls.dat'/,lenb/70/
      data hfile/'maptopog.dat'/,lenh/1080/
      data iu/42/,nrec/1080/,nwrd/2160/,lrec/-1/
c     ..
      save lrec,numbs,mbuf,mwrd,nlat

      if (lrec.lt.0) then
        lenr = lenb*4
        open (iu,file=bfile,access='direct',recl=lenr/rfac,
     &        status='old')
        lenm = lenb - 2
        numb = 0
        do 120 irec = 1,nrec
          read (iu,rec=irec,iostat=iok) mbit
	  if (iok.ne.0) go to 999
          numbs(irec) = numb
          mb = ichars(mbit(lenb),4)
          numb = numb + mb
          do 110 i = 1,lenm
  110     mbuf(i,irec) = mbit(i)
  120   continue
        close (iu)
        lenr = lenh*4
        open (iu,file=hfile,access='direct',recl=lenr/rfac,
     &        status='old')
        mwrd = nwrd*2
        lrec = 0
        mlat = nrec/2
      end if
c ...
      iele = 0
      isfc = 0
      ilat = nint(flat*6.)
      if (flat.lt.0.) ilat = iabs(ilat) + mlat
      ilon = -nint(flon*6.)
      if (ilon.lt.0) ilon = ilon + nwrd
      jlat = ilat + 1
      jlon = ilon + 1
      if (jlon.gt.nwrd) jlon = jlon - nwrd
      do 130 i = 1,lenm
  130 mbit(i) = mbuf(i,jlat)
      isfc = lbit(jlon,mbit)
      if (isfc.eq.0) return
      numb = 1
      if (ilon.eq.0) go to 150
      do 140 ibit = 1,ilon
  140 numb = numb + lbit(ibit,mbit)
  150 numb = numb + numbs(jlat)
      irec = (numb-1)/mwrd
      iwrd = numb - mwrd*irec
      irec = irec + 1
      if (irec.eq.lrec) go to 160
      read (iu,rec=irec,iostat=iok) lbuf
      if (iok.ne.0) go to 999
      lrec = irec
  160 iele = mova2i(lbuf(iwrd)) + 128
c ... convert from 100's of feet to meters
      fele = 30.48*float(iele)
      iele = nint(fele)
 999  return 
      end
      integer function ichars(cword,len)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    ichars      converts characters into integer
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: turns character string of specified length into integer.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    ichars(cword,len)
c   input argument list:                                               
c     cword    - character*1 array argument.
c     len      - integer argument holding length of cword.
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
      character*1 cword(len)
      lchars = 0
      do ic = len , 1 , -1 
        ibeg = (len - ic) * 8
        icmove = mova2i(cword(ic))
        call mvbits(icmove,0,8,lchars,ibeg)
      enddo
      ichars = lchars
      return
      end
cfpp$ expand(ichars,lbit)
      integer function lansea(rlat,rlon,level)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    lansea      calculates topography, land/sea status    
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: calculates topography, land/sea status from latitude and 
c   longitude.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    lansea(rlat,rlon,level)
c   input argument list:                                               
c     rlat     - integer argument containing scaled latitude.
c     rlon     - integer argument containing scaled longitude.
c
c   output argument list:
c     level    - integer argument containing scaled topography.
c
c   input files:
c     unit 41  - binary low-resolution topography file
c
c   subprograms called:
c     library:
c       oneblib- lbit     ichars
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 01.01.90

      include 'rfac.inc'
      integer,parameter::real_64=selected_real_kind(15,307)
      
c ... scalar arguments ..
      integer ilat,ilon,level
c ...
c ... local scalars ..
      real(real_64) rlat,rlon
      real slat,slon,xlon
      integer iopn,iu,lan,last,lat,lenr,lon
      character*12 name
c ...
c ... local arrays ..
      character*4 iflag(12),kelev(192)
      character*2 ielev(360)
c ...
c ... external functions ..
      integer lbit,ichars
      external lbit,ichars
c ...
c ... intrinsic functions ..
      intrinsic float,max0
c ...
c ... equivalences ..
      equivalence (iflag(1),kelev(1)), (ielev(1),kelev(13))
c ...
c ... data statements ..
      data name/'lowtopog.dat'/,iu/41/,lenr/768/,last/0/,iopn/0/
c ...
      save iopn,kelev,last
 
      if (iopn.eq.0) then
        open (iu,recl=lenr/rfac,
     &        file=name,access='direct',status='old')
        iopn = 1
      end if
      lan = 0
      level = 0
      slat = rlat
      slon = rlon
      lat = slat + 1.
      if (slat.lt.0.) lat = slat
      lat = max0(lat,-87)
      lat = 91 - lat
      if (lat.eq.last) go to 100
      read (iu,rec=lat) kelev
      last = lat
  100 xlon = slon
      if (xlon.lt.0.) xlon = xlon + 360.
      lon = xlon
      if (lon.eq.360) lon = 0
      lon = lon + 1
      lan = lbit(lon,iflag)
      if (lan.ne.0) then
        ltemp = ichars(ielev(lon),2)
        if (btest(ltemp,15)) then
          ltemp = ior(ltemp,-65536)
        endif
        level = ltemp
      endif
      lansea = 2 - lan
      return
      end
cfpp$ expand(ichars)
      integer function lbit(j,array)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    lbit        extracts jth bit from character array
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: extracts jth bit from array of character*4.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    lbit(j,array)    
c   input argument list:                                               
c     j        - integer argument.
c     array    - character*4 array argument.
c
c   subprograms called:
c     library:
c       oneblib- ichars
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 06.12.84
 
c ... scalar arguments ..
      integer j
c ...
c ... array arguments ..
      character*4 array(*)
c ...
c ... local scalars ..
      integer ibit,jout,jw,jword,nbit
c ...
c ... external functions ..
      integer ichars
      external ichars
c ...
c ... intrinsic functions ..
      intrinsic btest
c ...
      jw = (j-1)/32
      nbit = j - jw*32
      jword = ichars(array(jw+1),4)
      ibit = 32 - nbit
      jout = 0
      if (btest(jword,ibit)) jout = 1
      lbit = jout
      return
      end
      integer function lchar(j,ldata)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    lchar       extracts jth char from character array
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: extracts jth character from array of character*1 and
c   returns the integer equivalent.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    lchar(j,ldata)    
c   input argument list:                                               
c     j        - integer argument.
c     ldata    - character*1 array argument.
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 13.05.92
 
c ... scalar arguments ..
      integer j
c ...
c ... array arguments ..
      character*1 ldata(*)
c ...
c ... intrinsic functions ..
      intrinsic ichar,btest,ior
c ...
      itemp = mova2i(ldata(j))
      if (btest(itemp,7)) then
        itemp = ior(itemp,-256)
      endif
      lchar = itemp
      return
      end
      integer function lngtst(a,x,b)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    lngtst      compares three real numbers
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: compares rightmost byte of three scaled floating
c   point numbers. returns "1" if first argument is greater than
c   third argument and second argument is not between them. also 
c   returns "1" if first argument is less than or equal to third
c   argument and second argument is between them. otherwise "0"
c   is returned.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    lngtst(a,x,b)    
c   input argument list:                                               
c     a        - real argument.
c     x        - real argument.
c     b        - real argument.
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 04.12.87
c ... scalar arguments ..
      real a,b,x
c ...
c ... local scalars ..
      real fac
      integer(8)  ja,jb,jx
      integer(8)  mask
c ...
c ... intrinsic functions ..
      intrinsic iand
c ...
c ... data statements ..
      data fac/182.0444/
      data mask/x'000000000000ffff'/
c ...
 
      ja = fac*a
      jx = fac*x
      jb = fac*b
      ja = iand(ja,mask)
      jx = iand(jx,mask)
      jb = iand(jb,mask)
      lngtst = 0
      if (ja.gt.jb) go to 50
      if (jx.ge.ja .and. jx.le.jb) lngtst = 1
      return
   50 if (jx.ge.ja .or. jx.le.jb) lngtst = 1
      return
      end
cfpp$ expand(ichars)
      integer function mbyt16(nb,iarray)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    mbyt16      converts character*2 element to integer
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: extracts (nb)th element from array of character*2 and
c   converts it to integer.  entry point mbyt16 propagates sign bit 
c   in result; entry point lbyt16 does not.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    mbyt16(nb,iarray)    
c   input argument list:                                               
c     nb       - integer argument.
c     iarray   - character*2 array argument.
c
c   subprograms called:
c     oneblib  - ichars
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 13.05.92
c ... scalar arguments ..
      integer nb
c ...
c ... array arguments ..
      character*2 iarray(*)
c ...
c ... local scalars ..
      integer jj,mask
c ...
c ... external functions ..
      integer ichars
      external ichars
c ...
c ... intrinsic functions ..
      intrinsic btest,ior,iand
c ...
c ... entry points ..
      integer lbyt16
c ...
c ... data statements ..
      data mask/65535/,negmsk/-65536/
c ...
c ...
      jj = ichars(iarray(nb),2)
      if(btest(jj,15)) then
        jj = ior(jj,negmsk)
      endif
      mbyt16 = jj
      return
c ...
      entry lbyt16(nb,iarray)
 
      jj = ichars(iarray(nb),2)
      jj = iand(jj,mask)
      lbyt16 = jj
      return
      end
cfpp$ expand(ichars)
      integer function mbyte(j,length,jarray)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    mbyte       converts bit string to integer
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: extracts bit string from array of character*4 and
c   converts it to integer.  entry point mbyte propagates sign bit 
c   in result; entry point lbyte does not.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    mbyte(j,length,jarray)
c   input argument list:                                               
c     j        - integer argument containing starting bit.
c     length   - integer argument containing number of bits.
c                maximum value: 32.
c     jarray   - character*4 array argument.
c
c   output files:
c     unit 06  - printout
c
c   subprograms called:
c     oneblib  - ichars
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 03.08.95
c ... scalar arguments ..
      integer j,length
c ...
c ... array arguments ..
      character*4 jarray(*)
c ...
c ... local scalars ..
c ... local scalars ..
      integer inleft,jbit,kompl,mflag,n,nlj,nrj,
     +        nsh,nword
      integer(8) item,jleft,jrite,mask
c ...
c ... local arrays ..
      integer(8)  kounts(33)
c ...
c ... external functions ..
      integer ichars
      external ichars
c ...
c ... intrinsic functions ..
      intrinsic iand,ior,mod
c ...
c ... entry points ..
      integer lbyte
c ...
c ... data statements ..
      data kounts/1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,
     +     16384,32768,65536,131072,262144,524288,1048576,2097152,
     +     4194304,8388608,16777216,33554432,67108864,134217728,
     +     268435456,536870912,1073741824,2147483648_8,0/
c ...
c ...
c ... entry mbyte
      mflag = 1
  110 nword = (j-1)/32 + 1
      if (length.lt.1 .or. length.gt.32) write (*,fmt=120) length
  120 format (' improper byte length in mbyte or lbyte',i10)
      nlj = mod(j-1,32)
      nrj = 32 - length - nlj
      if (nrj.lt.0) go to 150
      item = ichars(jarray(nword),4)
      kompl = 33 - length - nlj
      mask = -kounts(kompl)
      item = iand(item,mask)
      item = item/kounts(nrj+1)
      mask = kounts(length+1) - 1
      item = iand(item,mask)
  130 if (mflag.eq.0) go to 140
c ... means logical byte
      mbyte = item
      jbit = iand(kounts(length),item)
      if (jbit.eq.0) return
c ... need sign extension
      mask = -mask - 1
      item = ior(item,mask)
      mbyte = item
      return
c ...
      entry lbyte(j,length,jarray)
 
      mflag = 0
      go to 110
c ...
  140 lbyte = item
      return
c ...
c ... byte spans two words
  150 inleft = length + nrj
      mask = kounts(inleft+1) - 1
      jleft = ichars(jarray(nword),4)
      jleft = iand(jleft,mask)
      nsh = 1 - nrj
      jleft = jleft*kounts(nsh)
      n = 1 - nrj
      jrite = ichars(jarray(nword+1),4)
      kompl = 33 + nrj
      mask = -kounts(kompl)
      jrite = iand(jrite,mask)
      jrite = jrite/kounts(nrj+33)
      mask = kounts(n) - 1
      jrite = iand(jrite,mask)
      item = ior(jrite,jleft)
      mask = kounts(length+1) - 1
      go to 130
      end
      integer function nsat1b(inscid,jy)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    nsat1b      converts internal id to satellite no.
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: converts tovs1b internal spacecraft-id to noaa number.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    nsat1b(inscid,jy)    
c   input argument list:                                               
c     inscid   - integer argument containing internal spacecraft id.
c     jy       - integer argument containing two-digit year.
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c
c     Declare variables.
      integer inscid,jy
      integer nsat
c
c
c     The integer argument nsat1b is the index location in array
c     csat (defined) in the calling program for the spacecraft 
c     name which corresponds to the spacecraft id.  This name-id
c     matching is based on Table 2.0.4-3 of the NOAA POD User's
c     Guide (Aug 97 revision)
c
c     The below declarations need to be in the routine calling
c     function nsat1b in order to get the proper id-name matching.
c
c     character*8  csat(0:14)
c     data csat/5*' ','tiros-n','noaa-6','noaa-7','noaa-8','noaa-9',
c     +     'noaa-10','noaa-11','noaa-12','noaa-13','noaa-14'/
c     data cunk/'unknown'/
c     csat(0) = cunk
c
c
c*************************************************************8
c     Start nsat1b here.
c
c     Initialize output to zero (unknown spacecraft)
      nsat = 0
c
c
c     Given spacecraft id, set index for array csat of 
c     corresponding spacecraft name
      if (inscid.eq.1) then
        if (jy.gt.1987) then
c         NOAA-11
          nsat = 11
        else
c         TRIOS-N
          nsat = 5
        end if
c
      else if (inscid.eq.2) then
         if (jy.lt.1988) then
c           NOAA-6
            nsat = 6
         elseif (jy.eq.1993) then
c           NOAA-13
            nsat = 13
         else
c           UNKNOWN
            nsat = 0
         endif
c
      else if (inscid.eq.3) then
c        NOAA-14
         nsat=14
c
      else if (inscid.eq.4) then
c        NOAA-7 
         nsat = 7
c
      else if (inscid.eq.5) then
c        NOAA-12 
         nsat = 12
c
      else if (inscid.eq.6) then
c        NOAA-8
         nsat = 8
c
      else if (inscid.eq.7) then
c        NOAA-9
         nsat = 9
c
      else if (inscid.eq.8) then
c        NOAA-10
         nsat = 10
c
      else
c        UNKNOWN
         nsat = 0
      end if
c
c     Set index in array csat corresponding to
c     the given satellite id.
      nsat1b = nsat
c
c     End of function
      return
      end
      integer function nthday(ny,nm,nd)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    nthday      calculates day-of-year from year, month, day
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: calculates day-of-year from year, month, day.   
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    nthday(ny,nm,nd)
c   input argument list:                                               
c     ny       - integer argument containing two-digit year.
c     nm       - integer argument containing month of year.
c     nd       - integer argument containing day-of-month.
c
c   subprograms called:
c     library:
c       w3lib  - w3fs13  
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 29.09.79
c ... to get the nth day of the year, given calendar date
c ... scalar arguments ..
      integer nd,nm,ny
c ...
c ... local scalars ..
      integer nth
c ...
c ... external subroutines ..
      external w3fs13
c ...
      call w3fs13(ny,nm,nd,nth)
      nthday = nth
      return
      end
      subroutine pfcoef(noaa)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    pfcoef      reads in correction coefficients
c   prgmmr: katz,b           org: np2        date: 97-11-07             
c                                                                       
c abstract: reads in hirs-2 planck function, band correction,
c   l-w flux, and skin-temp coefficients into common blocks
c   for use by other routines.
c                                                                       
c program history log:                                                  
c   97-11-07  katz                                                 
c                                                                       
c usage:    call pfcoef(noaa,*)
c   input argument list:                                               
c     noaa     - integer argument containing satellite id.
c
c   input files:
c     unit 35  - binary file of correction coefficients.
c
c   output files:
c     unit 06  - printout
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c
c
c     Include machine dependent parameters.
      integer,parameter::real_32=selected_real_kind(6,37)
      integer,parameter::real_64=selected_real_kind(15,307)
c 
c     Declare/set parameters.
      integer lenc,nb,nt,lenp,lenf,lens,lent,iuc,lenr
      parameter (lenc=1200,nb=20,nt=2,lenp=nb* (nt+3))
      parameter (lenf=13,lens=5,lent=8)
      parameter (iuc=35,lenr=lenc*8)
c 
c     Declare varaibles.
      character*80 cfile
      integer noaa,i,ioff,irec,loff,moff,noff
      real(real_64) cflx(lenf),csun(lens),dh(nb),gh(nb),pbuf(lenp)
      real(real_64) tscc(lent),tsce(lent)
      real(real_64) cbuf(lenc)
c
c     Declare common blocks.
      common /cfdb/ctype
      common /flxcof/cflx
      common /gamhir/gh,dh
      common /plnctn/pbuf
      common /solcof/csun
      common /tsurfc/tscc
      common /tsurfe/tsce
c
c     Data statements
      data cfile/'hirsrtcf_ibm.dat'/
c
c
c*******************************************************************
c     Start pfcoef here.
c
c
c     Based on user specified coefficient database type, set
c     offset.  
         ioff = 100
c
c
c     Open unit to coefficient database.  In case of error, print
c     message and stop program.
      open (iuc,file=cfile,recl=lenr,access='direct',
     x     status='old',err=100)
 
c
c     Given satellite id, set record number for data to pull.
c     Read data into a buffer.  Close the unit.
      irec = (noaa-10)*nb
      read (iuc,rec=irec) cbuf
      close (iuc)
c
c
c     Extract information from buffer.  Load into appropriate 
c     arrays for later use.
c
c     Planck function coefficients.
      do i = 1,lenp
        pbuf(i) = cbuf(i)
      end do
c
c     These cofficients are NOT used (as of 2 Dec 97).
      do i = 1,nb
        gh(i) = cbuf(i+100)
        dh(i) = cbuf(i+120)
      end do
c
c     Coefficients used by function irefls
      do i = 1,lens
        csun(i) = cbuf(i+140)
      end do
c
c     Coefficients used by function lwflux
      loff = ioff + 45
      do i = 1,lenf
        cflx(i) = cbuf(i+loff)
      end do
c
c     These coefficients are NOT used (as of 2 Dec 97).
      moff = ioff + 80
      noff = ioff + 90
      do i = 1,lent
        tscc(i) = cbuf(i+moff)
        tsce(i) = cbuf(i+noff)
      end do
c
c
c     Successful end of routine.
      return
c
c
c     ERROR handling:  error opening data file
  100 write(*,'(''error opening file '',a12,'' in pfcoef'')') cfile
      call w3tage('BUFR_TRANHIRS2')
      call errexit(5)
      stop
c
c     ERROR handling:  invalid coefficient type specified by user
  200 write(*,'(''invalid coeff-type '',a4,'' passed to pfcoef'')')
     +     ctype
      call w3tage('BUFR_TRANHIRS2')
      call errexit(5)
      stop
c
c     End of routine
      end
      real(8) function planck(t,k)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    planck      does what function planck does
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: calculates planck function from temperature and 
c   hirs-2 channel number.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    planck(t,k)
c   input argument list:                                               
c     t        - real argument containing temperature.
c     k        - integer argument containing hirs-2 channel number.
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 05.12.79
      integer,parameter::real_32=selected_real_kind(6,37)
      integer,parameter::real_64=selected_real_kind(15,307)
 
c ... scalar arguments ..
      real(real_64) t
      integer k
c ...
c ... arrays in common ..
      real(real_64) fk1(20),fk2(20),fnu(20),tc(2,20)
c ...
c ... local scalars ..
      real(real_64) expn,tt
c ...
c ... intrinsic functions ..
      intrinsic exp
c ...
c ... common blocks ..
      common /plnctn/fnu,fk1,fk2,tc
c ...
      tt = tc(1,k) + tc(2,k)*t
      expn = exp(fk2(k)/tt) - 1.
      planck = fk1(k)/expn
      return
      end
      real(8) function rbiref(visr,mon,cosz)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    ibiref      does what function ibiref does
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: calculates ibiref. 
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    ibiref(ivisr,mon,cosz)
c   input argument list:                                               
c     ivisr    - integer argument.
c     mon      - integer argument.
c     cosz     - real argument.
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 30.01.80
      integer,parameter::real_32=selected_real_kind(6,37)
      integer,parameter::real_64=selected_real_kind(15,307)
c ... scalar arguments ..
      real(real_64) cosz,visr
      integer ivisr,mon
c ...
c ... local scalars ..
      real(real_64) b,solcon,v
c ...
c ... local arrays ..
      real(real_64) alfa(12)
c ...
c ... intrinsic functions ..
      intrinsic float,max0
c ...
c ... data statements ..
      data alfa/1.034d0,1.028d0,1.018d0,1.002d0,.985d0,
     x     .972d0,.965d0,.969d0,.983d0,.997d0,
     +     1.016d0,1.018d0/
      data solcon/22.407d0/
c ...
 
      v = 1.d-3*visr
      b = v/ (cosz*solcon*alfa(mon))
cx      ibiref = b*100. + .5
cx      ibiref = max0(ibiref,0)
      rbiref = b
      rbiref = max(rbiref,0.d0)
      return
      end
cfpp$ expand(bright,planck)
      real(8) function rfluor(rtb,cosz)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    ifluor      does what function ifluor does
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: calculates scaled quantity ifluor.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    ifluor(itb,cosz)
c   input argument list:                                               
c     itb      - integer argument containing scaled brightness 
c                temperature.
c     cosz     - real argument containing cosine of zenith angle.
c
c   subprograms called:
c     unique:    fs (statement function)
c     oneblib  - bright   planck
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 30.07.79
      integer,parameter::real_32=selected_real_kind(6,37)
      integer,parameter::real_64=selected_real_kind(15,307)
c ... scalar arguments ..
      real(real_64) cosz,itb
c ...
c ... local scalars ..
      real(real_64) con,rc,tc
      integer i,ic
c ...
c ... external functions ..
      real bright,planck
      external bright,planck
c ...
c ... intrinsic functions ..
      intrinsic max,float
c ...
c ... statement functions ..
      real fs
c ...
c ... data statements ..
      data ic/17/,con/.070d0/
c ...
c ...
      rc = planck(rtb,ic) - con*cosz
      rc = max(rc,.001d0)
      tc = bright(rc,ic)
      rfluor = tc
      return
      end
cfpp$ expand(planck)
      real(8) function rrefls(rtsol)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    irefls      does what function irefls does
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: calculates scaled quantity irefls.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    irefls(itsol)
c   input argument list:                                               
c     itsol    - integer argument.
c
c   subprograms called:
c     unique:    fs2 (statement function)
c     oneblib  - planck
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 18.12.91
      integer,parameter::real_32=selected_real_kind(6,37)
      integer,parameter::real_64=selected_real_kind(15,307)
c ... array arguments ..
cx      integer itsol(*)
      real(real_64) rtsol(*)
c ...
c ... arrays in common ..
      real(real_64) c(5)
c ...
c ... local scalars ..
      real(real_64) a,arg,f,t,w3,w4,tt,expn
      real(real_64) fk1(20),fk2(20),fnu(20),tc(2,20)
      integer i,i2,i3
c ...
c ... external functions ..
      real(real_64) planck
      external planck
c ...
c ... intrinsic functions ..
      intrinsic log,float
c ...
c ... common blocks ..
      common /solcof/c
      common /plnctn/fnu,fk1,fk2,tc
c ...
c ... statement functions ..
cx      real fs2
c ...
c ... data statements ..
      data i2/2/,i3/3/
c ...
      a = c(5)
cx      w4 = planck(rtsol(i2),18)
cx      w3 = planck(rtsol(i3),19)

      k = 18
      t = rtsol(2)
      tt = tc(1,k) + tc(2,k)*t
      expn = exp(fk2(k)/tt) - 1.
      w4 = fk1(k)/expn

      k = 19
      t = rtsol(3)
      tt = tc(1,k) + tc(2,k)*t
      expn = exp(fk2(k)/tt) - 1.
      w4 = fk1(k)/expn

      rrefls = rtsol(i2)
      arg = a*w4 - w3
      if (arg.le.0.) return
      f = log(arg)
      t = c(1) + f* (c(2)+f* (c(3)+f*c(4)))
      rrefls = t
      return
      end
      real(8) function rsolzn(rlat,rlon,idt,cosz)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    isolzn      calculates solar zenith angle 
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: calculates scaled solar zenith angle.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    isolzn(rlat,rlon,idt,cosz)
c   input argument list:                                               
c     rlat     - integer argument containing scaled latitude.
c     rlon     - integer argument containing scaled longitude.
c     idt      - integer array of three elements, containing 
c                (1) seconds of day, (2) day of year, and (3) year.
c     cosz     - real argument containing cosine of zenith angle.
c
c   subprograms called:
c     unique:    scale (statement function)
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 07.12.79
      integer,parameter::real_32=selected_real_kind(6,37)
      integer,parameter::real_64=selected_real_kind(15,307)
c
c ... scalar arguments ..
      real(real_64) cosz
      real(real_64) rlat,rlon
c ...
c ... array arguments ..
      integer idt(3)
c ...
c ... local scalars ..
      real(real_64) cosa,cosd,cosx,dpd,dtr,dyr,g,gha,
     x     sind,sinx,sinz,soza,ss,
     +     xlon,z,z2,za
      integer i,iday,init,isec,iyr,jlon,nite
c ...
c ... local arrays ..
      real(real_64) dcon(5),tcon(4)
c ...
c ... intrinsic functions ..
      intrinsic atan2,cos,float,mod,sin,sqrt
c ...
c ... statement functions ..
      real(real_64) scale
c ...
c ... data statements ..
      data dtr/.0174533d0/,init/1/,nite/90.00/
      data dcon/279.9348d0,1.914827d0,-.079525d0,.019938d0,-.00162d0/
      data tcon/-.03d0,-.12d0,.165d0,-.0008d0/
      save init,sind,cosd,g
c ...
c ... statement function definitions ..
 
      scale(i) = .01d0*float(i)
c ...
      isec = idt(1)
      if (init.eq.0) go to 110
      iday = idt(2)
      iyr = idt(3)
      dyr = 365.d0
      if (mod(iyr,4).eq.0) dyr = 366.d0
      dpd = 360.d0/dyr
      z = dpd*float(iday-1)
      za = z*dtr
      z2 = za*2.
      ss = dcon(1) + z + dcon(2)*sin(za) + dcon(3)*cos(za) +
     +     dcon(4)*sin(z2) + dcon(5)*cos(z2)
      sind = .39785d0*sin(ss*dtr)
      cosd = sqrt(1.-sind*sind)
      za = dtr*dpd* (float(iday)-dyr+286.d0)
      z2 = za*2.d0
      g = tcon(1)*sin(za) + tcon(2)*cos(za) + tcon(3)*sin(z2) +
     +    tcon(4)*cos(z2)
      init = 0
  110 sinx = sin(dtr*rlat)
      cosx = sqrt(1.-sinx*sinx)
      jlon = -100*rlon + 0.5
      if (jlon) 120,130,130
  120 jlon = 36000 + jlon
  130 xlon = dtr*scale(jlon)
      gha = dtr*15.d0* (g+12.d0+float(isec)/3600.d0)
      cosa = cos(xlon-gha)
      cosz = sind*sinx + cosd*cosx*cosa
      sinz = sqrt(1.-cosz*cosz)
      soza = atan2(sinz,cosz)/dtr
      if (cosz.le.0.) go to 140
      rsolzn = soza
      return
  140 rsolzn = nite
      return
      end
      SUBROUTINE TIMSAT(JYEAR, JMON, JDAY, JHOUR, JMIN, JSEC, 
     *                  IL   , IE  , 
     *                  JFYEAR,JFMON,JFDAY,JFHOUR,JFMIN,JFSEC)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK                                    
C                .      .    .                                       .  
C SUBPROGRAM:    TIMSAT      GENERATES TIME STAMP FOR TOVS1B OBS
C   PRGMMR: KATZ,B           ORG: NP2        DATE: 97-11-10             
C                                                                       
C ABSTRACT: GIVEN STARTING TIME OF ORBIT AND LINE AND SPOT NUMBER WITHIN
C   ORBIT, GENERATES TIME STAMP FOR INDIVIDUAL SPOTS IN HIRS-2 SCANS.
C                                                                       
C PROGRAM HISTORY LOG:                                                  
C   97-11-10  KATZ                                                 
C   99-09-23  X.Su modified just for hirs-2 data
C                                                                       
C USAGE:    CALL TIMSAT(JYEAR, JMON, JDAY, JHOUR, JMIN, JSEC,
C    *                  IL   , IE  , 
C    *                  JFYEAR,JFMON,JFDAY,JFHOUR,JFMIN,JFSEC)
C   INPUT ARGUMENT LIST:                                               
C     JYEAR    - INTEGER ARGUMENT CONTAINING ORBIT YEAR.
C     JMON     - INTEGER ARGUMENT CONTAINING ORBIT MONTH OF YEAR.
C     JDAY     - INTEGER ARGUMENT CONTAINING ORBIT DAY OF MONTH.
C     JHOUR    - INTEGER ARGUMENT CONTAINING ORBIT HOUR OF DAY.
C     JMIN     - INTEGER ARGUMENT CONTAINING ORBIT MINUTE OF HOUR.
C     JSEC     - INTEGER ARGUMENT CONTAINING ORBIT SECOND OF MINUTE.
C     IL       - INTEGER ARGUMENT CONTAINING ORBIT LINE NUMBER.
C     IE       - INTEGER ARGUMENT CONTAINING ORBIT SPOT NUMBER.
C
C   OUTPUT ARGUMENT LIST:                                               
C     JYEAR    - INTEGER ARGUMENT CONTAINING SPOT YEAR.
C     JMON     - INTEGER ARGUMENT CONTAINING SPOT MONTH OF YEAR.
C     JDAY     - INTEGER ARGUMENT CONTAINING SPOT DAY OF MONTH.
C     JHOUR    - INTEGER ARGUMENT CONTAINING SPOT HOUR OF DAY.
C     JMIN     - INTEGER ARGUMENT CONTAINING SPOT MINUTE OF HOUR.
C     JSEC     - INTEGER ARGUMENT CONTAINING SPOT SECOND OF MINUTE.
C
C   SUBPROGRAMS CALLED:
C     LIBRARY:
C       ONEBLIB- CALDAY   NTHDAY
C
C REMARKS:                                                              
C                                                                       
C ATTRIBUTES:                                                           
C   LANGUAGE: FORTRAN 77                                                
C   MACHINE:  CRAY                                                     
C                                                                       
C$$$  
c     Set parameters
      integer,parameter::real_32=selected_real_kind(6,37)
      integer,parameter::real_64=selected_real_kind(15,307)
c
c
      CHARACTER*4 CMON
      real(real_64) els
 
         ELS=(IL-1)*6.4d0 + (IE-1)*0.1d0 + JSEC
 
      
      JFSEC=MOD(ELS,60.d0) + 0.5d0
      IF(JFSEC.EQ.60) THEN
        ielm=int(els/60.d0) + JMIN + 1
        JFSEC=0
      ELSE
        IELM=INT(ELS/60.d0) + JMIN
      ENDIF
 
      JFMIN=MOD(IELM,60)

      IELH=IELM/60 + JHOUR
      JFHOUR=MOD(IELH,24)

      IELD=IELH/24

      IF(IELD.GE.1) THEN
        IF(JMON.EQ.12 .AND. (JDAY+IELD).GT.31) THEN
          IF(JYEAR.EQ.99) JYEAR=-1
          IYYDDD=(JYEAR+1)*1000 + JDAY + IELD - 31
        ELSE
          NDAY= NTHDAY(JYEAR,JMON,JDAY)
          IYYDDD=JYEAR*1000 + NDAY + IELD
        ENDIF
        CALL CALDAY(IYYDDD,JFYEAR,JFMON,JFDAY,CMON)
C       WRITE(*,*) 'CALDAY:YYDDD= ',IYYDDD,JFYEAR,JFMON,JFDAY
      ELSE
        JFDAY=JDAY
        JFMON=JMON
        JFYEAR=JYEAR
      ENDIF

      RETURN
      END
cfpp$ expand(ichars)
      real function xfloat(jb,iarray)
c$$$  subprogram documentation block                                    
c                .      .    .                                       .  
c subprogram:    xfloat      converts character*2 array to real
c   prgmmr: katz,b           org: np2        date: 97-11-05             
c                                                                       
c abstract: takes two consecutive elements of character*2 array and
c   forms a floating point number from them.
c                                                                       
c program history log:                                                  
c   97-11-05  katz                                                 
c                                                                       
c usage:    xfloat(jb,iarray)
c   input argument list:                                               
c     jb       - integer argument containing array location.
c     iarray   - character*2 array argument.
c
c   subprograms called:
c     oneblib  - ichars
c
c remarks:                                                              
c                                                                       
c attributes:                                                           
c   language: fortran 77                                                
c   machine:  cray                                                     
c                                                                       
c$$$                                                                    
c ... version of 05.10.79
c ... to make a floating number from two consecutive 16-bit
c ... integers, the first of which is the jb-th word of iarray.
 
c ... scalar arguments ..
      integer jb
      integer(8) mask
      data mask/x'ffffffff00000000'/
c ...
c ... array arguments ..
cx      character*2 iarray(2)
      character*2 iarray(2000)
      character*4 conv
c ...
c ... local scalars ..
      real xf
      integer j
      integer(8) jj
c ...
c ... local arrays ..
      integer in(2)
c ...
c ... external functions ..
      integer ichars
      external ichars
c ...
c ... intrinsic functions ..
      intrinsic btest,ior
c ...
      j = jb
      conv(1:2) = iarray(j)
      conv(3:4) = iarray(j+1)
      jj = ichars(conv,4)
      if (btest(jj,31_8)) then
        jj = ior(jj,mask) 
      endif
      xf = jj
      xfloat = xf
      return
      end
