C  Copyright(c) 1998, Space Science and Engineering Center, UW-Madison
C  Refer to "McIDAS Software Acquisition and Distribution Policies"
C  in the file  mcidas/data/license.txt

C *** $Id: kbxavhr.dlm,v 1.9 1999/05/07 21:16:16 krauss Rel $ ***

C$ Name:
C$      KBX_INI - Calibration module for AVHRR, initializes conversion type
C$
C$ Interface:
C$      INTEGER FUNCTION
C$      KBX_INI(CHARACTER*4 CIN, CHARACTER*4 COUT, INTEGER ARRAY IOPT(*))
C$
C$ Input:
C$      CIN - Input pixel type ('RAW')
C$      COUT - Output pixel type ('TEMP', 'REFL', 'RAD', 'RAW', 'BRIT')
C$      IOPT - Input calibration parameters
C$                 (array size set by NUMAREAOPTIONS)
C$               IOPT(1) - Source pixel size (1,2,4)
C$               IOPT(2) - Destination pixel size (1,2,4)
C$               IOPT(3) - unused
C$               IOPT(4) - unused
C$               IOPT(5) - unused
C$
C$ Output:
C$      To be backward compatible with the existing API, the return
C$      value does not indicate which error was encountered. You must
C$      use the trace facility to get more information.
C$
C$ Return values:
C$        0     - successful initialization of the calibration module
C$       -1     - error (undefined input or output pixel size or type)
C$
C$ Remarks:
C$      Note that for ADDE, messages cannot be sent to "stdout" pipe
C$      mixed with data.  Messages are sent to trace file if the user
C$      has set TRACE=1 on the command line.
C$
C$ Categories: 
C$      image 
C$      display 
C$      calibration 
C$      met/science 


      INTEGER FUNCTION KBXINI(CIN,COUT,IOPT)

      implicit none           !  All variables must be declared


      ! Input Parameters:

      character*4  cin        !  Input pixel type (RAW)
      character*4  cout       !  Output pixel type (TEMP, REFL, RAD, RAW, BRIT)
      integer      iopt(*)    !  Input calibration parameters


      INCLUDE      'areaparm.inc'     !  Global declarations for McIDAS areas
                                      !     (defines NUMAREAOPTIONS)



      ! Symbolic constants & shared data

      integer      itype                    !  Input pixel data type
      integer      jtype                    !  Output pixel data type
      integer      jopt(NUMAREAOPTIONS)     !  Calibration parameters

      common/AVHXX/itype,jtype,jopt


      ! Local variable




C     Move pixel sizes to common for later access by KBX-CAL
C     Input must be 2-byte RAW pixels containing 10-bit data values
C     which are left shifted 5 bits (divide by 32 to get true value)

      call movw(NUMAREAOPTIONS,iopt,jopt)
      if(iopt(1).ne.1.and.iopt(1).ne.2.and.iopt(1).ne.4) then
         go to 900
      end if
      if(iopt(2).ne.1.and.iopt(2).ne.2.and.iopt(2).ne.4) then
         go to 900
      end if


C     Set the input pixel type (currently only RAW is supported)

      itype=0
      if(cin.eq.'RAW') itype=1
      if(itype.eq.0) then
         go to 900
      end if


C     Set the output pixel type (RAD & REFL share type 2)
C     RAW is type 1, which is a "do nothing" call for kbx-cal, returning RAW
C     Consequently, any undefined output type requested will get RAW pixels
C     returned and no error is flagged.

      jtype=1
      if(cout.eq.'RAD' .or. cout.eq.'REFL') then
         jtype=2
      end if
      if(cout.eq.'TEMP') then
         jtype=3
         if(jopt(2).lt.2) go to 900
      end if
      if(cout.eq.'BRIT') jtype=4
      if(cout.eq.'MODB') jtype=4

      KBXINI=0

      return


C     Undefined input or output pixel size or type

900   continue
      KBXINI=-1

      return

      end


C============================================================================



C$ Name:
C$      KBX_CAL - Converts input DN to output of correct byte size & dimensions
C$
C$ Interface:
C$      INTEGER FUNCTION
C$      KBX_CAL(INTEGER PFX(*), INTEGER IDIR(*), INTEGER NVAL,
C$                                    INTEGER BAND, INTEGER IBUF(*) )
C$
C$ Input:
C$      PFX  - Line header/prefix associated with the IBUF data
C$      IDIR - Area directory associated with the IBUF data
C$      NVAL - Number of input pixels in the IBUF data
C$      BAND - Area band number for the IBUF data
C$
C$
C$ Input and Output
C$      IBUF - Buffer containing the pixel data numbers to be converted
C$                 (The size of the input and output pixels in IBUF is
C$                  determined by the first two IOPT values in the call
C$                  to KBX_INI.)
C$
C$
C$ Return values:
C$       0      - success
C$      -1      - failure to build a calibration lookup table
C$                 (output will always be returned in IBUF, even with an
C$                  invalid or unbuilt table, but it will be undefined,
C$                  so it is important to always check the return value)
C$
C$ Remarks:
C$      This function uses the initialized common block to perform the
C$      conversion from input pixel data (usually from an area, via a
C$      call to REDARA) to output units of the correct byte size and
C$      physical dimensions.
C$
C$      If a lookup table for the conversion hasn't been built, this routine
C$      will first build one, when possible, using information in the area's
C$      directory IDIR, the line prefix PFX, and the band number BAND.
C$
C$ Categories: 
C$      image 
C$      display 
C$      calibration 
C$      met/science 



      INTEGER FUNCTION KBXCAL(PFX,IDIR,NVAL,BAND,IBUF)


      implicit none             !  All variables must be declared


      ! Input Parameters:

      integer pfx(*)            !  Input array of calibration constants
                                !     from line header for IBUF data
      integer idir(*)           !  Area directory buffer (64 words minimum)
      integer nval              !  Number of pixels to process from input
      integer band              !  Band number (1-6)


      ! Input & Output Parameters:

      integer ibuf(*)           !  I/O array containing pixels to be modified
                                !     (will contain converted values at end)


      ! Symbolic constants & global shared data in common

      INTEGER    MAXAVHRRBUFFERS    ! number of buffers for use by the AVHRR
                                    ! calibration routines for storing
                                    ! calibration tables in single band mode.
      PARAMETER (MAXAVHRRBUFFERS=6)


      INTEGER    MAXAVHRRBANDS      ! number of bands the AVHRR/3 instrument
                                    ! on the POES satellite transmits
      PARAMETER (MAXAVHRRBANDS=6)

      INCLUDE 'areaparm.inc'         !  Global declarations for McIDAS areas
                                     !     (defines NUMAREAOPTIONS)


      integer itype
      integer jtype
      integer JOPT(NUMAREAOPTIONS)   !  Calibration parameters for conversion
      common/AVHXX/itype,jtype,jopt


      character*4 caltyp             !  Calibration type for breakpoint table
      common/BRKPNT/caltyp           !  Common block for BRKSET table type


      real bbtk                      !  bbtk = internal target temperature
      common/BBT/bbtk


      ! External functions

      integer ic
      integer lit


      ! Local variables

      integer calpos                            !  position in prefix cal block
      integer i                                 !  loop index
      integer icc1                              !  calibration constants 1-8
      integer icc2
      integer icc3
      integer icc4
      integer icc5
      integer icc6
      integer icc7
      integer icc8
      integer ides                              !  destination pixel size
      integer isou                              !  sounrce pixel size
      integer iss                               !  spacecraft id number
      integer istat                             !  status return value
      integer j                                 !  loop index
      integer kband                             !  band number
      integer lband                             !  band number
      integer kcc1(MAXAVHRRBANDS)               !  BAND=ALL integer constants
      integer kcc2(MAXAVHRRBANDS)
      integer kcc3(MAXAVHRRBANDS)
      integer kcc4(MAXAVHRRBANDS)
      integer kcc5(MAXAVHRRBANDS)
      integer kcc6(MAXAVHRRBANDS)
      integer kcc7(MAXAVHRRBANDS)
      integer kcc8(MAXAVHRRBANDS)
      integer kstbnd(MAXAVHRRBANDS)             !  BAND=ALL previous band
      integer kstcal                            !  previous conversion type
      integer kstcc1(MAXAVHRRBANDS)             !  BAND=ALL previous constants
      integer kstcc2(MAXAVHRRBANDS)
      integer kstcc3(MAXAVHRRBANDS)
      integer kstcc4(MAXAVHRRBANDS)
      integer kstcc5(MAXAVHRRBANDS)
      integer kstcc6(MAXAVHRRBANDS)
      integer kstcc7(MAXAVHRRBANDS)
      integer ktype                             !  RAW_TO conversion type
      integer length                            !  length of line prefix
      integer level                             !  band map position
      integer levpos                            !  band map address (lev block)
      integer lstbnd(MAXAVHRRBUFFERS)           !  previous band
      integer lstcal(MAXAVHRRBUFFERS)           !  previous conversion type
      integer lstcc1(MAXAVHRRBUFFERS)           !  previous cal constants 1-7
      integer lstcc2(MAXAVHRRBUFFERS)
      integer lstcc3(MAXAVHRRBUFFERS)
      integer lstcc4(MAXAVHRRBUFFERS)
      integer lstcc5(MAXAVHRRBUFFERS)
      integer lstcc6(MAXAVHRRBUFFERS)
      integer lstcc7(MAXAVHRRBUFFERS)           !  constant 8 currently unused
      integer mband(MAXAVHRRBANDS)              !  band map array
      integer nbands                            !  z-resolution (# of bands)
      integer ncal                              !  size of cal section
      integer ncod                              !  size of line prefix
      integer ndoc                              !  size of doc section
      integer nlev                              !  size of lev section
      integer*2    itemp(2048)                  !  integer temperature*100

      integer      table(1024,MAXAVHRRBUFFERS)  !  lookup table for each band
      integer*2    kable(1024,MAXAVHRRBANDS)
      equivalence  (table,kable)                !  equivalence to save space


      ! initialized variables

      integer kstsss(MAXAVHRRBANDS)
      integer lstsss(MAXAVHRRBUFFERS)

      data kstsss /   MAXAVHRRBANDS*0 /
      data lstsss / MAXAVHRRBUFFERS*0 /



C     Get the number of number of bands present and the spacecraft ID
      nbands=idir(14)
      iss=idir(3)                 ! valid spacecraft are 62(NOAA-12),
                                  !                      64(NOAA-14),
                                  !                      65(NOAA-15)

      if(band.eq.6) then          ! band 6 constants in prefix at position 3
         lband=3
      else
         lband=band
      end if


C     Get the line prefix section lengths in bytes
      ncod=0
      if(idir(36).ne.0) ncod=4
      ndoc=idir(49)
      ncal=idir(50)
      nlev=idir(51)
      length=ncod+ndoc+ncal+nlev
      length=(length+3)/4


C     Prepare the line-sensitive calibration constants for the requested band
C     This must be done for every access, even if tables are already built,
C     since we have to see if the constants have changed from the last table

      calpos=ncod+ndoc
      if(band.eq.lit('ALL')) go to 500

      if(ncal.eq.160) then                 !  AVHR line prefix
         calpos=4+calpos+(lband-1)*32      !     8 constants per band
         call movw(1,pfx(calpos/4),icc1)
         call movw(1,pfx(1+calpos/4),icc2)
         call movw(1,pfx(2+calpos/4),icc3)
         call movw(1,pfx(3+calpos/4),icc4)
         call movw(1,pfx(4+calpos/4),icc5)
         call movw(1,pfx(5+calpos/4),icc6)
         call movw(1,pfx(6+calpos/4),icc7)
         call movw(1,pfx(7+calpos/4),icc8)
         bbtk=0.01 * icc7                  ! icc7 is target temp for bands 3-5

      else if(ncal.eq.32) then             !  AVHR line prefix assumed
         calpos=4+calpos                   !     8 constants, single band
         call movw(1,pfx(calpos/4),icc1)
         call movw(1,pfx(1+calpos/4),icc2)
         call movw(1,pfx(2+calpos/4),icc3)
         call movw(1,pfx(3+calpos/4),icc4)
         call movw(1,pfx(4+calpos/4),icc5)
         call movw(1,pfx(5+calpos/4),icc6)
         call movw(1,pfx(6+calpos/4),icc7)
         call movw(1,pfx(7+calpos/4),icc8)
         bbtk=0.01 * icc7                  ! icc7 is target temp for bands 3-5

      else                                 ! Unknown line prefix format
         go to 900
      end if


C     Get the input and output pixel sizes
      isou=jopt(1)
      ides=jopt(2)

      go to (40,10,20,30) jtype      !  build or switch tables for outputs
                                     !  jtype=1 is a "do nothing"
                                     !  for RAW or unrecognized output types


C     If we switch bands, we regenerate the tables.  After that, we continue
C     to use the same table until a line prefix with a reference minor frame
C     is encountered, whereupon we build another table.  This normally occurs
C     every five lines in the image.  If we do not have any changes in the
C     spacecraft id, band number, output calibration type, or slope-offset
C     constants, and if this is not the first pass through the calibration
C     module for this band, we can omit the table regeneration, because
C     nothing has changed.  Such stability will often be the case for many
C     images.


C     RAW to RAD or REFL conversion
   10 do i=1,MAXAVHRRBUFFERS
         if(lstsss(i).eq.0) then        ! generate table on first pass
            lstsss(i)=iss
            lstbnd(i)=band
            lstcc1(i)=icc1
            lstcc2(i)=icc2
            lstcc3(i)=icc3
            lstcc4(i)=icc4
            lstcc5(i)=icc5
            lstcc6(i)=icc6
            lstcc7(i)=icc7
            lstcal(i)=jtype
            call RAWTO(0,band,iss,icc1,icc2,icc3,icc4,icc5,icc6,
     &                                           icc7,table(1,i),istat)
            if(istat.ne.0) go to 900
            go to 15
         end if
C        If same band, ss, and jtype, then we have to regenerate the table
C        only if encountering any new cc constants
         if(band.eq.lstbnd(i).and.iss.eq.lstsss(i).and.
     &                                       lstcal(i).eq.jtype) then
            if(lstcc1(i).ne.icc1.or.lstcc2(i).ne.icc2.or.
     &                   lstcc3(i).ne.icc3.or.lstcc4(i).ne.icc4) then
               lstcc1(i)=icc1
               lstcc2(i)=icc2
               lstcc3(i)=icc3
               lstcc4(i)=icc4
               lstcc5(i)=icc5
               lstcc6(i)=icc6
               lstcc7(i)=icc7
               call RAWTO(0,band,iss,icc1,icc2,icc3,icc4,icc5,icc6,
     &                                           icc7,table(1,i),istat)
               if(istat.ne.0) go to 900
            end if
            go to 15
         end if
      end do

C     Error exit -- Couldn't find old table or build a new one
      go to 900

C     Now do the table lookup
   15 continue
      call mavhtb(nval,isou,ides,ibuf,table(1,i))
      go to 100


C     RAW to TEMP conversion
   20 do i=1,MAXAVHRRBUFFERS
         if(lstsss(i).eq.0) then        ! generate table on first pass
            lstsss(i)=iss
            lstbnd(i)=band
            lstcc1(i)=icc1
            lstcc2(i)=icc2
            lstcc3(i)=icc3
            lstcc4(i)=icc4
            lstcc5(i)=icc5
            lstcc6(i)=icc6
            lstcc7(i)=icc7
            lstcal(i)=jtype
            call RAWTO(1,band,iss,icc1,icc2,icc3,icc4,icc5,icc6,
     &                                           icc7,table(1,i),istat)
            if(istat.ne.0) go to 900
            go to 25
         end if
C        If same band, ss, and jtype, then we have to regenerate the table
C        only if encountering any new cc constants
         if(band.eq.lstbnd(i).and.iss.eq.lstsss(i).and.
     &                                      lstcal(i).eq.jtype) then
            if(lstcc1(i).ne.icc1.or.lstcc2(i).ne.icc2.or.
     &                   lstcc3(i).ne.icc3.or.lstcc4(i).ne.icc4) then
               lstcc1(i)=icc1
               lstcc2(i)=icc2
               lstcc3(i)=icc3
               lstcc4(i)=icc4
               lstcc5(i)=icc5
               lstcc6(i)=icc6
               lstcc7(i)=icc7
               call RAWTO(1,band,iss,icc1,icc2,icc3,icc4,icc5,icc6,
     &                                           icc7,table(1,i),istat)
               if(istat.ne.0) go to 900
            end if
            go to 25 
         end if
      end do

C     Error exit -- Couldn't find old table or build a new one
      go to 900

C     Now do the table lookup
   25 continue
      call mavhtb(nval,isou,ides,ibuf,table(1,i))
      go to 100


C     RAW to BRIT or MODB conversion
   30 do i=1,MAXAVHRRBUFFERS
         if(lstsss(i).eq.0) then        ! generate table on first pass
            lstsss(i)=iss
            lstbnd(i)=band
            lstcc1(i)=icc1
            lstcc2(i)=icc2
            lstcc3(i)=icc3
            lstcc4(i)=icc4
            lstcc5(i)=icc5
            lstcc6(i)=icc6
            lstcc7(i)=icc7
            lstcal(i)=jtype
            call RAWTO(2,band,iss,icc1,icc2,icc3,icc4,icc5,icc6,
     &                                           icc7,table(1,i),istat)
            if(istat.ne.0) go to 900
            go to 35
         end if
C        If same band, ss, and jtype, then we have to regenerate the table
C        only if encountering any new cc constants
         if(band.eq.lstbnd(i).and.iss.eq.lstsss(i).and.
     &                                      lstcal(i).eq.jtype) then
            if(lstcc1(i).ne.icc1.or.lstcc2(i).ne.icc2.or.
     &                   lstcc3(i).ne.icc3.or.lstcc4(i).ne.icc4) then
               lstcc1(i)=icc1
               lstcc2(i)=icc2
               lstcc3(i)=icc3
               lstcc4(i)=icc4
               lstcc5(i)=icc5
               lstcc6(i)=icc6
               lstcc7(i)=icc7
               call RAWTO(2,band,iss,icc1,icc2,icc3,icc4,icc5,icc6,
     &                                           icc7,table(1,i),istat)
               if(istat.ne.0) go to 900
            end if
            go to 35
         end if
      end do

C     Error exit -- Couldn't find old table or build a new one
      go to 900

C     Now do the table lookup
   35 continue
      call mavhtb(nval,isou,ides,ibuf,table(1,i))
      go to 100


C     RAW or unknown cal, just change pixel sizes in ibuf if necessary
   40 do i=1,MAXAVHRRBUFFERS
         if(lstsss(i).eq.0) then        ! generate tables on first pass only
            lstsss(i)=iss
            do j=1,1024
               table(j,1) = (j-1)*32
            end do
         end if
         call mavhtb(nval,isou,ides,ibuf,table(1,1))
      end do


  100 continue
      KBXCAL=0
      return



C     Build calibration tables for all bands in the area at once
  500 continue
      levpos=ncod+ndoc+ncal
      level=0
      kband=0
C     If more than one band, all 5 sets of constants are generated
C     and the mband array determines which constants are applied
      if(nbands.gt.1) nbands=5
      do i=1,nbands
         level=ic(pfx,levpos)
         if(level.eq.0) go to 506
         kband=kband+1
         mband(kband)=level
         levpos=levpos+1

         if(ncal.eq.160.or.ncal.eq.32) then        !  AVHR line prefix assumed
            calpos=4+calpos+(i-1)*32               !     8 constants per band
            call movw(1,pfx(calpos/4),kcc1(kband))
            call movw(1,pfx(1+calpos/4),kcc2(kband))
            call movw(1,pfx(2+calpos/4),kcc3(kband))
            call movw(1,pfx(3+calpos/4),kcc4(kband))
            call movw(1,pfx(4+calpos/4),kcc5(kband))
            call movw(1,pfx(5+calpos/4),kcc6(kband))
            call movw(1,pfx(6+calpos/4),kcc7(kband))
            call movw(1,pfx(7+calpos/4),kcc8(kband))
            calpos=calpos+32

         else                                 !  Unrecognized prefix
         end if
      end do

  506 continue
      isou=jopt(1)
      ides=jopt(2)
      go to (510,900) itype

  510 if(jtype.eq.1) then
         go to 900
      end if
      ktype=-1
      if(jtype.eq.2) ktype=0
      if(jtype.eq.3) ktype=1
      if(jtype.eq.4) ktype=2
      if(ktype.eq.-1) then
         go to 900
      end if

      do i=1,kband
         if(kstsss(i).eq.0) then        ! generate table on first pass
            kstsss(i)=iss
            kstbnd(i)=mband(i)
            kstcc1(i)=kcc1(i)
            kstcc2(i)=kcc2(i)
            kstcc3(i)=kcc3(i)
            kstcc4(i)=kcc4(i)
            kstcc5(i)=kcc5(i)
            kstcc6(i)=kcc6(i)
            kstcc7(i)=kcc7(i)
            kstcal=ktype
            call RAWTO(ktype,mband(i),iss,kcc1(i),kcc2(i),kcc3(i),
     &                    kcc4(i),kcc5(i),kcc6(i),kcc7(i),itemp,istat)
            if(istat.ne.0) go to 900
            call movh(1024,itemp(2),2,kable(1,i),1)
            go to 515
         end if

C        If same band, ss, and ktype, then we have to regenerate
C        the table only if encountering any new cc constants
         if(mband(i).eq.kstbnd(i).and.iss.eq.kstsss(i).and.
     &                                          kstcal.eq.ktype) then
            if(kstcc1(i).ne.kcc1(i).or.kstcc2(i).ne.kcc2(i).or.
     &             kstcc3(i).ne.kcc3(i).or.kstcc4(i).ne.kcc4(i)) then
               kstcc1(i)=kcc1(i)
               kstcc2(i)=kcc2(i)
               kstcc3(i)=kcc3(i)
               kstcc4(i)=kcc4(i)
               kstcc5(i)=kcc5(i)
               kstcc6(i)=kcc6(i)
               kstcc7(i)=kcc7(i)
               call RAWTO(ktype,mband(i),iss,kcc1(i),kcc2(i),kcc3(i),
     &              kcc4(i),kcc5(i),kcc6(i),kcc7(i),itemp,istat)
               if(istat.ne.0) go to 900
               call movh(1024,itemp(2),2,kable(1,i),1)
            end if
         end if

  515    continue
      end do

C     Successful generation of tables for all bands
  600 continue
      KBXCAL=0
      return


C     Error in building calibration table(s)
  900 continue
      KBXCAL=-1
      return

      end


C============================================================================


C$ Name:
C$      KBX_OPT - Returns auxiliary parameters for setup or sets internal state
C$
C$ Interface:
C$      INTEGER FUNCTION
C$      KBX_OPT(CHARACTER*4 CFUNC, INTEGER IIN(*), INTEGER IOUT(*) )
C$
C$ Input:
C$      CFUNC   - Option or function to be executed
C$                             ('KEYS', 'BRKP', 'INFO')
C$      IIN     - Array of input parameters
C$                  (The number of input and output parameters is determined
C$                   by the option or function executed.  The calling routine
C$                   is responsible for having arrays large enough to send
C$                   and receive output for the option or function requested.)
C$
C$ Output:
C$      IOUT    - Array of output parameters
C$                  (The number of input and output parameters is determined
C$                   by the option or function executed.  The calling routine
C$                   is responsible for having arrays large enough to send
C$                   and receive output for the option or function requested.)
C$
C$ Return values:
C$       0      - success
C$      -1      - Invalid function requested (doesn't exist)
C$      -2      - Invalid band for the spacecraft or instrument (doesn't exist)
C$      -3      - Error in lookup table generation (table cannot be set up)
C$
C$ Remarks:
C$      Check the imgprobe.f code to see how the various functions are used.
C$
C$ Categories: 
C$      image 
C$      display 
C$      calibration 
C$      met/science 



      INTEGER FUNCTION KBXOPT(CFUNC,IIN,IOUT)
C     Sets internal flags or states, returns parameter or state data

      IMPLICIT NONE             !  All variables must be declared

C     Input parameters:
      CHARACTER*4 CFUNC         !  Option or function descriptor
      INTEGER IIN(*)            !  Input parameters

C     Output Parameters:
      INTEGER IOUT(*)           !  Output parameters

C     ERROR CODES -1 INVALID OPTION CFUNC
C                 -2 INVALID BAND ('INFO')
C                 -3 ERROR IN BREAKPOINT TABLE OR MISSING TABLE KEYS='BRKP'


C     Symbolic constants & global shared data in common

      INCLUDE 'areaparm.inc'         !  Global declarations for McIDAS areas
                                     !     (defines NUMAREAOPTIONS)


      character*4 caltyp             !  Calibration type for breakpoint table
      common/BRKPNT/caltyp           !  Common block for BRKSET table type


C     External functions

      integer brkset                 !  Sets breakpoint table values
      integer ischar                 !  Numeric ASCII value of character byte
      integer lit                    !  Four byte integer representing CHAR*4



C     Local variables

      character*8 cfile              !  Breakpoint table name for SU
                                     !    (stored in frame dir words 38-39)


      ! Symbolic constants & shared data

      integer itype                  !  Input type
      integer jtype                  !  Output type
      integer jopt(NUMAREAOPTIONS)   !  Pixel sizes for input and output
      common/AVHXX/itype,jtype,jopt



C     1)  "KEYS"  option:                     ! Identified valid pixel types
      if(cfunc.eq.'KEYS') then
         if(iin(4).le.2.or.iin(4).eq.6) then
            iout(1)=3
            iout(2)=lit('RAW ')
            iout(3)=lit('REFL')
            iout(4)=lit('BRIT')
         else
            iout(1)=4
            iout(2)=lit('RAW ')
            iout(3)=lit('RAD ')
            iout(4)=lit('TEMP')
            iout(5)=lit('BRIT')
         end if
         if (ischar(iin(38)).eq.1) then
            call movwc(iin(38),cfile)
            if (brkset(cfile,caltyp).ne.0) then
               KBXOPT = -3
               return
            end if
         end if
         KBXOPT = 0


C     2)  "BRKP"  option:                     ! Check for SU file existence
      else if (cfunc.eq.'BRKP') then
         call movwc(iin(1),cfile)
         if (brkset(cfile,caltyp).ne.0) then
            KBXOPT = -3
            return
         end if
         KBXOPT = 0


C     3)  "INFO"  option:            ! Return dimensions and scaling factors
      else if (cfunc.eq.'INFO') then

C     Check for valid band
         if (iin(1).lt.1.or.iin(1).gt.12) then
            KBXOPT = -2
            return
         end if

         if(iin(1).le.2.or.iin(1).eq.6) then
            iout(1)=3
            iout(2)=lit('RAW ')
            iout(3)=lit('REFL')
            iout(4)=lit('BRIT')
            iout(5)=lit('    ')
            iout(6)=lit(' %  ')
            iout(7)=lit('    ')
            iout(8)=1
            iout(9)=100
            iout(10)=1
         else
            iout(1)=4
            iout(2)=lit('RAW ')
            iout(3)=lit('RAD ')
            iout(4)=lit('TEMP')
            iout(5)=lit('BRIT')
            iout(6)=lit('    ')
            iout(7)=lit('mW**')
            iout(8)=lit('K   ')
            iout(9)=lit('    ')
            iout(10)=1
            iout(11)=1
            if(iin(1).eq.3) iout(11)=1000
            if(iin(1).eq.4.or.iin(1).eq.5) iout(11)=100
            iout(12)=10
            iout(13)=1
            KBXOPT = 0
         end if
      else
         KBXOPT = -1
      end if
      return

      end

C============================================================================


C$ Name:
C$      RAWTO - Converts AVHRR/3 raw data number to radiance, albedo, or temp
C$
C$ Interface:
C$      SUBROUTINE
C$      RAWTO(INTEGER TYPE,INTEGER BAND,INTEGER ISS,INTEGER ICC1,
C$          INTEGER ICC2, INTEGER ICC3, INTEGER ICC4,INTEGER ICC5,
C$          INTEGER ICC6, INTEGER ICC7,INTEGER TABLE,ISTAT)
C$
C$ Input:
C$      type  - calibration type desired for output
C$      band  - instrument band number
C$      iss   - spacecraft id number
C$      icc1  - first calibration constant
C$      icc2  - second calibration constant
C$      icc3  - third calibration constant
C$      icc4  - fourth calibration constant
C$      icc5  - fifth calibration constant
C$      icc6  - sixth calibration constant
C$      icc7  - seventh calibration constant
C$
C$ Output:
C$      table(*)  - calibration table that was generated
C$      istat - return status
C$                  0   - success
C$                 -1   - Invalid input parameters -- table not generated
C$
C$ Return values:
C$      none (see istat above)
C$
C$ Remarks:
C$    Conversion from RAW to RAD or REFL(type 0), TEMP (type 1),
C$    or BRIT(type2) applies to AVHRR only.
C$
C$    Calibrated values are stored with scale factors:
C$
C$       albedoes (reflectance factors) are between 0 and 10000
C$          for bands 1,2, and 3A  (i.e. per cent x 100)
C$
C$       radiances are in units of milliwatts/meter**2/steradian/cm
C$          times 100 for bands 4 and 5, and times 1000 for band 3B.
C$
C$ Categories: 
C$      image 
C$      display 
C$      calibration 
C$      met/science 


      SUBROUTINE RAWTO(TYPE,BAND,ISS,ICC1,ICC2,ICC3,ICC4,ICC5,
     &                                          ICC6,ICC7,TABLE,ISTAT)

      implicit none             !  All variables must be declared

      include 'hex80.inc'

      ! Input parameters:

      integer type
      integer band
      integer iss
      integer icc1
      integer icc2
      integer icc3
      integer icc4
      integer icc5
      integer icc6
      integer icc7


      ! Output Parameters:

      integer table(*)
      integer istat


      ! Symbolic constants & shared data

      INCLUDE 'areaparm.inc'

      integer itype
      integer jtype
      integer jopt(NUMAREAOPTIONS)
      common/AVHXX/itype,jtype,jopt

      character*4 caltyp                 !  SU file's calibration type
      common/BRKPNT/caltyp


      ! External functions

      integer brkval
      integer GRYSCL
      real AVHTMP


      ! Local variables

      integer i                        ! do loop index variable
      integer jbnd                     ! band index
      integer jsfac                    ! output scaling factor to apply
      integer jss                      ! spacecraft index
      real albval                      ! nominal reflectance values
      real albvalx                     ! NOAA-14 and earlier reflectance
      real albvaly                     ! reflectance from icc1-4 constants
      real calval                      ! unscaled calibrated pixel value
      real cc1                         ! calibration constant 1
      real cc2                         ! calibration constant 2
      real cc3                         ! calibration constant 3
      real cc4                         ! calibration constant 4
      real cc5                         ! calibration constant 5
      real cc6                         ! calibration constant 6
      real cc7                         ! calibration constant 7
      real radval                      ! temporary radiance value
      real radvaly                     ! radiance from NOAA-15 non-lin corr
      real radlin                      ! linear radiance from cc1 & cc2
      real radsp                       ! radiance of space
      real radsq                       ! square of linear radiance
      real radnl                       ! non-linear radiance
      real rawval                      ! raw data value
      real x                           ! x value of transfer function
      real y                           ! y value of transfer function
      character*80 string


      ! Initialized variables

      integer isfac(6)                 ! integer scaling for output

      real a(3,3)                      ! Nonlinear IR cal constant a(jbnd,jss)
      real b(3,3)                      ! Nonlinear IR cal constant b(jbnd,jss)
      real c(3,3)                      ! Nonlinear IR cal constant c(jbnd,jss)

      data isfac / 100,100,1000,100,100,100 /  ! these must match kbx_opt

      data a     / 0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0004524, 0.0002811 /
      data b     / 0.0,0.0,0.0,0.0,0.0,0.0,0.0, -0.0932,   -0.0659   /
      data c     / 0.0,0.0,0.0,0.0,0.0,0.0,0.0, 4.76,      3.83      /



      istat = 0

C     Calculate jss, the AVHRR spacecraft index into the data tables
      if(iss.eq.62) then
         jss=1
      else if(iss.eq.64) then
         jss=2
      else if(iss.eq.65) then
         jss=3
      else
      end if

      jbnd=band-2
      cc1=icc1/100000.0      ! slopes and intercepts are scaled by 100000
      if(band.eq.3) cc1=icc1/10000000.0     ! except band 3 slope
      cc2=icc2/100000.0
      cc3=icc3/100000.0
      cc4=icc4/100000.0
      cc5=icc5             ! raw space and target averages are unscaled
      cc6=icc6
      cc7=icc7/100.0       ! temp is scaled only by 100
      write(string,'(7F10.3)') cc1,cc2,cc3,cc4,cc5,cc6,cc7
      call m0sxtrce(string)
      if(icc1.eq.0.or.icc2.eq.0) then
         istat = -1
      end if
      if(icc1.eq.hex80.or.icc2.eq.hex80) then
C     A single band (probably 3 or 6) is missing from this line, so we
C     set the table to black output
         do i=1,1024
            table(i) = 0
         end do
      end if
      jsfac=isfac(band)              ! add band dependent scaling for McIDAS

C     Using modified brightness table against raw values
C     if caltyp is set, the SU tables should be used
      if ((caltyp.eq.'RAW') .and.(jtype.eq.4)) then
         do i=1,1024
            table(i)= brkval(real(i*32))
         end do
C        In this case, we are finished
         go to 999
      end if

C     Visible bands 1, 2, and 3A(labeled 6)
      if(band.le.2.or.band.eq.6) then

C        Type BRIT is albedo converted to 0-255
C        using formula : brit = 25.5 albedo**1/2  (0 < albedo < 100)
C        type rad is albedo times 100 (0-10000)
         do i=1,1024
            calval=i-1
            if(jss.le.2) then
c              NOAA-14 and earlier uses fixed calibration constants
c              for a single linear transformation
               albvalx=(calval*cc1+cc2)
               albvalx=amax1(0.,albvalx)
               albvalx=amin1(100.,albvalx)
               albval=(calval*cc1+cc2)
            else
C              NOAA-15 and later use a bi-linear transfer function
C              that stretches the low albedos more and provides
C              better low light and twilight resolution
               x=(cc2-cc4)/(cc3-cc1)
               y=cc1*x + cc2
               if( i.le.int(x) ) albvaly=(calval*cc1+cc2)
               if( i.gt.int(x) ) albvaly=(calval*cc3+cc4)
               albvaly=amax1(0.,albvaly)
               albvaly=amin1(100.,albvaly)
C              Neither slope can be negative or missing for NOAA-15
C              ...and for band 6 this might imply that band 3 was present
C              so we would want to output a black screen
               if(icc1.le.0.or.icc3.le.0.or.icc3.eq.HEX80) albvaly = 0
C              Theoretical values can be used as defaults in a pinch
C              and provide a check on the transmitted values of cc1-cc4
C              when using TRACE=1 output
               if(band.eq.6) then
                  if( i.le.500 ) albval=12.5 * (i-1)/499.0
                  if( i.gt.500 ) albval=12.5 + 87.5*(i-500)/500.0
               else if(band.lt.3) then
                  if( i.le.500 ) albval=25.0 * (i-1)/499.0
                  if( i.gt.500 ) albval=25.0 + 75.0*(i-500)/500.0
               end if
            end if
            albval=amax1(0.,albval)
            albval=amin1(100.,albval)
C           for operational use, we use albvaly, and force DN=0 to 0 albedo
C           so missing data gives a black screen
            if(jss.gt.2) then
               albval=albvaly
               if(i.eq.1) albval=0
            else
               albval=albvalx
               if(i.eq.1) albval=0
            end if

            if ((caltyp.eq.'REFL').and.(jtype.eq.4)) then
C              if caltyp is set, the SU tables should be used
               table(i) = brkval(albval)
            else if ((caltyp.eq.'BRIT').and.(jtype.eq.4)) then
C              if caltyp is set, the SU tables should be used
               table(i) = brkval(25.5*sqrt(albval))
            else if(type.eq.0) then
C              if caltyp is not set, we apply default scaling to REFL
               table(i)=0.5+albval*jsfac
            else if(type.eq.2) then
C              if caltyp is not set, we return BRIT
               table(i)=0.5+25.5*sqrt(albval)
            else
C              we should never get here
               istat = -1
            end if
         end do
      else

C        Infrared bands 3B, 4, and 5
         radsp=cc1*cc5+cc2                      ! radiance of space
         do i=1,1024
            rawval = i-1
            radval = cc1*rawval + cc2           ! linear radiance

            ! NOAA-12, -14 correct the "linear" temperature, not the radiance
            ! but for NOAA-15 we correct the radiance instead, including
            ! a correction for the measured radiance of space
            radlin = cc1*rawval + cc2           ! linear radiance for NOAA-15
            if(jbnd.eq.4) radlin = radlin - 4.50
            if(jbnd.eq.5) radlin = radlin - 3.61
            radsq  = radlin*radlin              ! radiance squared
            radnl = a(jbnd,jss)*radsq + b(jbnd,jss)*radlin + c(jbnd,jss)
            radvaly = radlin + radnl            ! non-linear radiance NOAA-15
            if( jss.eq.3 ) radval = radvaly

            radval = amax1(0.0,radval)          ! remove negative radiances

C           We look for requested use of SU tables (MODB: jtype=4) first
            if ((caltyp.eq.'RAD').and.(jtype.eq.4)) then
               table(i) = brkval(radval)
            else if ((caltyp.eq.'TEMP').and.(jtype.eq.4)) then
               calval = AVHTMP(jss,jbnd,radval)
               table(i) = brkval(calval)
            else if ((caltyp.eq.'BRIT/SU').and.(jtype.eq.4)) then
               calval = AVHTMP(jss+7,jbnd,radval)
               calval = GRYSCL(calval)
               table(i) = brkval(calval)
            else
C              No SU table lookup
               calval=radval*jsfac           ! Radiance (appropriately scaled)
               if(type.ge.1) then
                  calval=AVHTMP(jss+7,jbnd,radval)           ! TEMP conversion
               end if
               if(type.eq.1) then
                  calval=calval*10                           ! TEMP scaling
               end if
               if(type.eq.2) calval=GRYSCL(calval)           ! BRIT conversion
               table(i)=nint(calval)                         ! convert to int
            end if
C           Slope for IR bands must always be negative, or it's bad
C           There must be a IR target temperature (bbtk) present as well
            if(cc1.ge.0.0.or.icc1.eq.HEX80.or.
     &         cc7.eq.0.0.or.icc7.eq.HEX80) then
               table(i) = 0
            end if
         end do

C        Force 0,1023 to 0 output (this should be missing data -- unphysical)
C        so we output a black image for missing data, i.e. no band 3 or 6 data
C        returned by server
         table(1)=0
         table(1024)=0

      end if

999   continue

      return

      end


C============================================================================


C$ Name:
C$      AVHTMP - Transforms AVHRR IR radiance values to brightness temperatures
C$
C$ Interface:
C$      REAL FUNCTION
C$      AVHTMP(INTEGER ISS, INTEGER N, REAL VALUE)  (RCD)
C$
C$ Input:
C$      ISS     - (I) Input satellite code#
C$      N       - AVHRR Band No. - 2 (N=1 for band 3, i.e. IR Band No.)
C$      VALUE   - (I*2) Input radiance
C$
C$ Output:
C$      none
C$
C$ Return values:
C$      AVHTMP  - infrared brightness temperature corresponding to radiance
C$
C$ Categories: 
C$      image 
C$      display 
C$      calibration 
C$      met/science 


      REAL FUNCTION AVHTMP(iss,n,value)

      implicit none             !  All variables must be declared

      integer iss
      integer n
      real value


      ! Symbolic constants & shared data

      integer    NUMSAT
      parameter (NUMSAT=10)


      real bbtk       ! Internal calibration target (blackbody) temp.(deg.K)
      common/BBT/bbtk


      ! External functions

      real TDINT

      ! Local variables

      integer kss
      real britmp
      real cor
      real dte
      real te
      character*80 string


      ! Initialized variables

      integer nbbt
      integer ntt
      real c1
      real c2
      real veci(14)
      real vecj(3)
      real cnl14(14,3)
      real cnl24(14,3)
      real cnl15(14,3)
      real cnl25(14,3)
      real cnl35(14,3)
      real hwn(3,NUMSAT)
      real b2(3,NUMSAT)
      real b1(3,NUMSAT)


C     C1,C2=constants used to compute effective black body temp
      data C1 / 1.1910659E-5 /
      data C2 / 1.438833 /

      data nbbt /  3 /
      data ntt  / 14 /


C     VECI,VECJ,CNL are used to compute nonlinearity correction
C     VECI = 14 values of scene brightness temperature
C     VECJ =  3 values of internal target (blackbody) temperature
      data veci     /205.,215.,225.,235.,245.,255.,265.,
     &               275.,285.,295.,305.,310.,315.,320./
      data vecj     /283.16,288.16,293.16/


C     Nonlinearity corrections -- CNLIJ, I=CH.1,2,3 J=SAT.4,5
C     where CH.1=BAND4, CH.2=BAND5, CH.3=BAND3
C     CNL = correction, indices are scene temp, bb temp, ch, sat


C *** NOAA-12, CHANNEL 4 -- NON-LINEARITY CORRECTION TABLE

      DATA CNL14       /-1.6,-1.2,-1.3,-1.1,-1.2,-1.0,-0.7,
     &                  -0.4,+0.2, 0.8, 1.6, 2.0, 2.6, 3.2,
     &                  -1.8,-1.7,-1.7,-1.6,-1.4,-1.2,-1.0,
     &                  -0.8,-0.2, 0.5, 1.4, 1.9, 2.4, 2.9,
     &                  -1.3,-1.5,-1.6,-1.5,-1.6,-1.5,-1.2,
     &                  -1.1,-0.5, 0.1, 0.8, 1.3, 1.7, 2.3/

C *** NOAA-12, CHANNEL 5 -- NON-LINEARITY CORRECTION TABLE

      DATA CNL24       /-1.2,-1.2,-1.0,-0.9,-0.6,-0.5,-0.4,
     &                  -0.2, 0.1, 0.4, 0.7, 0.8, 0.8, 0.8,
     &                  -1.2,-1.2,-1.1,-0.9,-0.8,-0.5,-0.4,
     &                  -0.3,-0.1, 0.2, 0.6, 0.7, 0.8, 0.8,
     &                  -1.2,-1.2,-1.2,-1.0,-0.9,-0.6,-0.5,
     &                  -0.4,-0.2, 0.1, 0.4, 0.6, 0.7, 0.8/

C *** Note about noaa-14 correction tables.  Values given in the
C ***  correction table did not include values for differing internal
C ***  target temperatures.  Values in tables below have been duplicated
C ***  for all internal target temperatures BTTK.

C *** NOAA-14, CHANNEL 4 -- NON-LINEARITY CORRECTION TABLE
      DATA CNL15       / 5.58, 4.23, 3.04, 2.03, 1.35, 0.80, 0.39,
     &                   0.13,-0.02,-0.14, 0.05, 0.15, 0.27, 0.43,
     &                   5.58, 4.23, 3.04, 2.03, 1.35, 0.80, 0.39,
     &                   0.13,-0.02,-0.14, 0.05, 0.15, 0.27, 0.43,
     &                   5.58, 4.23, 3.04, 2.03, 1.35, 0.80, 0.39,
     &                   0.13,-0.02,-0.14, 0.05, 0.15, 0.27, 0.43/

C *** NOAA-14, CHANNEL 5 -- NON-LINEARITY CORRECTION TABLE
      DATA CNL25       / 2.22, 1.65, 1.14, 0.77, 0.47, 0.25, 0.09,
     &                  -0.00,-0.05,-0.03, 0.03, 0.09, 0.16, 0.25,
     &                   2.22, 1.65, 1.14, 0.77, 0.47, 0.25, 0.09,
     &                  -0.00,-0.05,-0.03, 0.03, 0.09, 0.16, 0.25,
     &                   2.22, 1.65, 1.14, 0.77, 0.47, 0.25, 0.09,
     &                  -0.00,-0.05,-0.03, 0.03, 0.09, 0.16, 0.25/

C *** NOAA-14, CHANNEL 3 -- NON-LINEARITY CORRECTION TABLE
      DATA CNL35       /-11.95,-7.92,-5.62,-2.58,-1.21,-0.66,-0.39,
     &                   -0.19,-0.09,-0.03, 0.01, 0.02, 0.02, 0.02,
     &                  -11.95,-7.92,-5.62,-2.58,-1.21,-0.66,-0.39,
     &                   -0.19,-0.09,-0.03, 0.01, 0.02, 0.02, 0.02,
     &                  -11.95,-7.92,-5.62,-2.58,-1.21,-0.66,-0.39,
     &                   -0.19,-0.09,-0.03, 0.01, 0.02, 0.02, 0.02/

C HALF-AREA WAVENUMBER OF THE CHANNEL SPECTRAL RESPONSE FUNCTION

      DATA HWN           /2651.41,921.33,921.33,
     &                    2673.70,951.49,951.49,
     &                    2686.70,928.59,841.54,
     &                    2653.69,915.59,915.59,
     &                    2691.95,930.75,845.84,
     &                    2673.03,911.15,911.15,
     &                    2683.61,928.71,841.78,
     &                    2653.39,921.92,837.55,
     &                    2658.25,929.70,833.92,
     &                    2695.9473,925.4075,839.8979/        ! NOAA-15

C FIRST REGRESSION COEFFICIENT

      DATA B2           /1.0023381,1.0017971,1.0017971,
     &                   1.0019172,1.0032497,1.0032491,
     &                   1.0024137,1.0014342,1.0011819,
     &                   1.0019587,1.0013928,1.0013928,
     &                   1.0025310,1.0013650,1.0011120,
     &                   1.0022960,1.0011321,1.0011321,
     &                   1.0025470,1.0013690,1.0009830,
     &                   1.0025850,1.0016330,1.0010740,
     &                   1.0029470,1.0016730,1.0008500,
     &                   1.0019890,1.0012830,1.0009770/   ! NOAA-15

C SECOND REGRESSION COEFFICIENT

      DATA B1           /-1.8188,-0.5888,-0.5888,
     &                   -1.7930,-3.0440,-3.0440,
     &                   -1.9895,-0.5730,-0.4061,
     &                   -1.8300,-0.5480,-0.5480,
     &                   -1.9310,-0.5420,-0.3950,
     &                   -1.8026,-0.5578,-0.5578,
     &                   -1.7992,-0.5032,-0.2365,
     &                   -1.9173,-0.5548,-0.3227,
     &                   -1.9571,-0.5148,-0.1482,
     &                   -1.624481,-0.338243,-0.304856/    ! NOAA-15



C     Calculate effective black body temp (te) using inverse Planck function:

C     Use the latest satellite ID if the current one isn't yet available
      kss=min(iss,NUMSAT)

      if(value.le.0.0) then     ! negative radiances give a zero temperature
         AVHTMP=0.0
         return
      end if
      dte=alog(1+c1*hwn(n,kss)**3/value)   ! calculate denominator
      if(dte.le.0.0) then
         AVHTMP=0.0
         return
      end if
      te=c2*hwn(n,kss)/dte     ! inverse Planck function

      if(iss.ge.8) then

C        Calculate the brightness temperature using the above regression tables
         britmp=b1(n,kss)+b2(n,kss)*te
         AVHTMP=britmp

C        If black body temp=0 or nonlinear correction unavailable (or
C        unneeded for NOAA-15 or later), return
         if(bbtk.eq. 0.0 .or. iss.ne.kss .or. iss.ge.10) return

C        Get and apply any nonlinearity correction
C        to thermal IR bands for satellites NOAA-12 & -14
         if(iss.eq.8.and.n.eq.2)
     &      cor=TDINT(vecj,nbbt,veci,ntt,cnl14,bbtk,britmp)
         if(iss.eq.8.and.n.eq.3)
     &      cor=TDINT(vecj,nbbt,veci,ntt,cnl24,bbtk,britmp)
         if(iss.eq.9.and.n.eq.1)
     &      cor=TDINT(vecj,nbbt,veci,ntt,cnl35,bbtk,britmp)
         if(iss.eq.9.and.n.eq.2)
     &      cor=TDINT(vecj,nbbt,veci,ntt,cnl15,bbtk,britmp)
         if(iss.eq.9.and.n.eq.3)
     &      cor=TDINT(vecj,nbbt,veci,ntt,cnl25,bbtk,britmp)
         AVHTMP=britmp+cor

      else

C        Negative radiance returned for undefined satellites (numbered 1-7).
C        The AVHR calibration module is not designed for NOAA-11 and earlier,
C        although we keep the index numbering the same for backward comparisons
C        with TIRO calibration constants.  Note that the new line prefix
C        structure and the slope-offset scaling for radiances are incompatible
C        with TIRO, so we have removed all the old unneeded CNL constants.
        
         AVHTMP=-9999.9

      end if

      return

      end


C============================================================================


C$ Name:
C$      TDINT   - Returns temperature correction interpolated from x-y grid
C$
C$ Interface:
C$      REAL FUNCTION
C$      TDINT( REAL VECJ,INTEGER NCOLS, REAL VECI, INTEGER NROWS,
C$                                       REAL GRID, REAL XIN, REAL YIN)
C$
C$ Input:
C$      VECJ    -  X values corresponsing to grid points
C$      NCOLS   -  number of columns in vecj array
C$      VECI    -  Y values corresponding to grid points
C$      NROWS   -  number of rows in veci array
C$      GRID    -  temperature correction grid
C$      XIN     -  X value at which to interpolate
C$      YIN     -  Y value at which to interpolate
C$
C$ Output:
C$      none
C$
C$ Return values:
C$      TDINT   - interpolated temperature correction
C$
C$ Categories: 
C$      calibration 
C$      met/science 


      REAL FUNCTION TDINT(vecj,ncols,veci,nrows,grid,xin,yin)

      implicit none             !  All variables must be declared

      ! Input parameters:

      integer ncols     ! dimension variables must appear before arrays
      real vecj(ncols)
      integer nrows     ! dimension variables must appear before arrays
      real veci(nrows)
      real grid(nrows,ncols)
      real xin
      real yin


      ! Local variables

      integer i
      integer i1
      integer i2
      integer j
      integer j1
      integer j2
      real dvx
      real dvy
      real dx
      real dy
      real v1
      real v2

      ! Initialized variables

      integer np
      data    np /0/


      np=np+1
      if(np.gt.30000)np=1

      if(xin.le.vecj(1))then
         j1=1
         j2=1
         go to 15
      end if

      if(xin.ge.vecj(ncols))then
         j1=ncols
         j2=ncols
         go to 15
      end if

      do 10 j=2,ncols
      if(xin.gt.vecj(j))go to 10
         j1=j-1
         j2=j
         go to 15
10    continue
15    continue

      if(yin.le.veci(1))then
         i1=1
         i2=1
         go to 25
      end if

      if(yin.ge.veci(nrows))then
         i1=nrows
         i2=nrows
         go to 25
      end if

      do 20 i=2,nrows
      if(yin.gt.veci(i))go to 20
         i1=i-1
         i2=i
         go to 25
20    continue
25    continue

      if(j1.eq.j2)then
         v1=grid(i1,j1)
         v2=grid(i2,j1)
      else
         dx=vecj(j2)-vecj(j1)
         dvx=grid(i1,j2)-grid(i1,j1)
         v1=grid(i1,j1)+(xin-vecj(j1))*dvx/dx
         dvx=grid(i2,j2)-grid(i2,j1)
         v2=grid(i2,j1)+(xin-vecj(j1))*dvx/dx
      end if

      if(i1.eq.i2)then
         TDINT=v1
      else
         dvy=v2-v1
         dy=veci(i2)-veci(i1)
         TDINT=v1+(yin-veci(i1))*dvy/dy
      end if

      return

      end


C============================================================================


C $ integer function gryscl(tempk)
C $ convert a real brightness temperature to an integer BRIT grey scale
C $ tempk = (r) input  temperature in degrees Kelvin
C $$ gryscl = converted brightness value (0-255)
C $ Remarks:
C $    Note that this function is similar to the existing function of the
C $    same name in the McIDAS library, but uses a real temperature value
C $    as input instead of a BRIT input value generated from temperature.

      INTEGER FUNCTION GRYSCL (tempk)

      implicit none
      real tempk

      ! INITIALIZED VARIABLES

      integer con1
      integer con2
      real tlim

      data con1 / 418 /
      data con2 / 660 /
      data tlim / 242.0 /

C convert a brightness temperature to a grey scale
C     tempk---temperature in degrees kelvin
      if (tempk .ge. tlim) goto 100
         GRYSCL = min0 (con1 - ifix (tempk), 255)
         goto 200
100   continue
         GRYSCL = max0 (con2 - ifix (2 * tempk), 0)
200   continue
      return

      end


C============================================================================


