C Copyright(c) 2006, 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: kbxfcin.dlm,v 1.10 2024/07/31 21:05:02 tommyj Exp $ ***

      INTEGER FUNCTION KBXINI(CIN,COUT,IOPT)

      IMPLICIT NONE
      CHARACTER*4 CIN
      CHARACTER*4 COUT
      INTEGER IOPT(*)

      ! symbolic constants & shared data

      INCLUDE 'areaparm.inc'

      INTEGER ITYPE
      INTEGER JOPT(NUMAREAOPTIONS)
      INTEGER CALFLG
      INTEGER CALARR(416)
      integer           enh_str_flag    ! =1 for new BRIT stretch
      character*8       enh_str_name    ! stretch table name
      COMMON/M0FCI/ITYPE,JOPT,CALFLG,CALARR,
     &              enh_str_flag, enh_str_name

      ! external functions

      ! local variables

      CALL MOVW(NUMAREAOPTIONS,IOPT,JOPT)

      ITYPE=0
      CALFLG = 0
      if(CIN.eq.'RAW'.and. COUT.eq.'RAD ') ITYPE=1
      if(CIN.eq.'RAW'.and. COUT.eq.'TEMP') ITYPE=2
      if(CIN.eq.'RAW'.and. COUT.eq.'BRIT') ITYPE=3
      if(CIN.eq.'RAW'.and. COUT.eq.'REFL') ITYPE=4
      if(ITYPE.eq.0) GOTO 900

c --- Switch off new stretch; call kbxopt to switch on

      enh_str_flag   = 0  ! By default, new stretch is off
      enh_str_name   = ' '

      KBXINI=0
      RETURN

900   CONTINUE
      KBXINI=-1
      RETURN
      END
C
C--------------------------------------------------------------------

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

      IMPLICIT NONE
      INTEGER PFX(*)
      INTEGER IDIR(*)
      INTEGER NVAL
      INTEGER BAND
      INTEGER*2 IBUF(*)

      ! symbolic constants & shared data

      INCLUDE 'areaparm.inc'

      integer iret

      INTEGER ITYPE, BAND_OFF, BAND_VIS
      INTEGER JOPT(NUMAREAOPTIONS)
      INTEGER CALFLG
      INTEGER CALARR(416)
      integer           enh_str_flag    ! =1 for new BRIT stretch
      character*8       enh_str_name    ! stretch table name

      integer image_band, fillvalue

      real*8 warm_scale_factor, warm_add_offset
      real*8 scale_factor, add_offset
      real*8 ru_conversion_coefficient
      real*8 wavenumber
      real*8 bt_coefficient_a, bt_coefficient_b
      real*8 bt_constant_c1, bt_constant_c2
      real*8 refl_solar_irradiance, earth_sun_distance
      real*8 nom, denom, radiance

      real btemp

C--- DAVEP
C      character*20 DAVEP

      COMMON/M0FCI/ITYPE,JOPT,CALFLG,CALARR,
     &              enh_str_flag, enh_str_name

      ! external functions

      INTEGER           NINT
      INTEGER           m0graybrkset
      INTEGER           m0grayscale
      INTEGER           isenhstr

      ! local variables

      INTEGER           THIS
      INTEGER           IDES
      INTEGER           ISOU
      INTEGER           I
      INTEGER           COUNT
      INTEGER           SS
      REAL              PI
      REAL              XBRIT
      REAL              XREFL

      DATA THIS/-9999/

      PI  =  4.D0*ATAN(1.D0)

      BAND_OFF = (BAND-1)*26+1
      BAND_VIS = 0
      IF (BAND.LE.8) THEN
         BAND_VIS = 1
      ENDIF

C
C --- get SS number form AREA Directory
C
      ss = idir(3)

C --- get radiometric constants from CAL header
      IF( THIS.NE.IDIR(33) ) THEN

c --- Setup for EXPanded BRIT values

         if(enh_str_flag .eq. 1) then

            iret = isenhstr(ss, band, enh_str_name)
            if( iret .eq. 0) then
               iret = m0graybrkset(enh_str_name)
            endif

         endif

         THIS = IDIR(33)

C ------ Check to see if calibration block passed in through kbxopt
         IF(CALFLG .NE. 0) THEN
C---            CALL MOVW(416, CALARR, CALARR)
         ELSE
            CALL ARAGET(IDIR(33),IDIR(63), 1664,CALARR)  
         ENDIF

      image_band = calarr(band_off+0)

         warm_scale_factor = (1.*calarr(band_off+1))
     &       / (1.*calarr(band_off+2))
         warm_add_offset = (1.*calarr(band_off+3))
     &       / (1.*calarr(band_off+4))
         scale_factor = (1.*calarr(band_off+5))
     &       / (1.*calarr(band_off+6))
         add_offset = (1.*calarr(band_off+7))
     &       / (1.*calarr(band_off+8))
         ru_conversion_coefficient = (1.*calarr(band_off+9))
     &       / (1.*calarr(band_off+10))

         wavenumber = (1.*calarr(band_off+11))
     &       / (1.*calarr(band_off+12))
         bt_coefficient_a = (1.*calarr(band_off+13))
     &       / (1.*calarr(band_off+14))
         bt_coefficient_b = (1.*calarr(band_off+15))
     &       / (1.*calarr(band_off+16))
         bt_constant_c1 = (1.*calarr(band_off+17))
     &       / (1.*calarr(band_off+18))
         bt_constant_c2 = (1.*calarr(band_off+19))
     &       / (1.*calarr(band_off+20))

         refl_solar_irradiance = (1.*calarr(band_off+21))
     &       / (1.*calarr(band_off+22))
         earth_sun_distance = (1.*calarr(band_off+23))
     &       / (1.*calarr(band_off+24))

         fillvalue = calarr(band_off+25)

         nom = bt_constant_c2 * wavenumber

         ISOU=JOPT(1)      ! length in bytes of input data
         IDES=JOPT(2)      ! length in bytes of output data
      ENDIF

C --- See SatPy reader fci_l1c_nc.py for calibration functions

      DO I = 1, NVAL
         COUNT = IBUF(I)

C ------------ Always calculate radiance
C              Could skip this if you know you will bail later
         if (count.eq.0) then
            radiance = 0
         else
            IF (BAND.EQ.9) THEN
               IF (count.eq.fillvalue) THEN
                  radiance = 0
               ELSEIF ((2 ** 12 - 1 < count) .AND. (count <= 2 ** 13 - 1)) THEN
                  radiance = (count * warm_scale_factor) + warm_add_offset
               ELSE
                  radiance = (count * scale_factor) + add_offset
               ENDIF
            ELSE
               IF (count.ne.fillvalue) THEN
                  radiance = (count * scale_factor) + add_offset
               ELSE
                  radiance = 0
               ENDIF
            ENDIF
         endif

         IF (BAND_VIS.EQ.1) THEN
C --------- Visible
C
            IF (ITYPE.eq.1) THEN
C ------------ RAW->RAD
               IBUF(I) = nint(radiance*100.)
               IF(IBUF(I).LT.  0) IBUF(I)=0
            ELSEIF (ITYPE.eq.2) THEN
C ------------ RAW->TEMP, can't do temperature for this channel
               IBUF(I) = 0

            ELSEIF (ITYPE.eq.3) THEN
C ------------ RAW->BRIT
              XREFL = 100. * radiance * PI * earth_sun_distance ** 2
     &         / refl_solar_irradiance
              iret = m0grayscale('REFL', xrefl, xBRIT)
              IBUF(I)=NINT(XBRIT)
              IF(IBUF(I).LT.  0) IBUF(I)=0
              IF(IBUF(I).GT.255) IBUF(I)=255

            ELSEIF (ITYPE.eq.4) THEN
C ------------ RAW->REFL
              XREFL = 100. * radiance * PI * earth_sun_distance ** 2
     &         / refl_solar_irradiance
              IBUF(I)=NINT(XREFL*100.)
              IF(IBUF(I).LT.  0) IBUF(I)=0

            ENDIF

         ELSE
C --------- Infra-red
C
            IF (ITYPE.eq.1) THEN
C ------------ RAW->RAD
               IBUF(I) = nint(radiance*100.)

            ELSE IF (ITYPE.eq.2) THEN
C ------------ RAW->TEMP
               IF (radiance.EQ.0) THEN
                  IBUF(I) = 255
               ELSE

                  denom = bt_coefficient_a
     &             * log( 1. + ((bt_constant_c1 * (wavenumber ** 3))
     &             / radiance) )

                  btemp = (nom / denom)
     &             - (bt_coefficient_b / bt_coefficient_a)

                  IBUF(I) = nint(btemp * 100.)
               ENDIF

            ELSEIF (ITYPE.eq.3) THEN
C ------------ RAW->BRIT
               IF (radiance.EQ.0) THEN
                  IBUF(I) = 255
               ELSE

                  denom = bt_coefficient_a
     &             * log( 1. + ((bt_constant_c1 * (wavenumber ** 3))
     &             / radiance) )

                  btemp = (nom / denom)
     &             - (bt_coefficient_b / bt_coefficient_a)

                  iret = m0grayscale('TEMP', btemp, xBRIT)

                  IBUF(I)=NINT(XBRIT)
                  IF(IBUF(I).LT.  0) IBUF(I)=0
                  IF(IBUF(I).GT.255) IBUF(I)=255
               ENDIF

            ELSEIF (ITYPE.eq.4) THEN
C ------------ RAW->REFL, can't do reflectance for this channel
               IBUF(I) = 0

            ENDIF

         ENDIF

      ENDDO

c --- pack the array
      CALL MPIXEL(NVAL,ISOU,IDES,IBUF)

      KBXCAL=0
      RETURN

C --- READ ERROR
999   CALL M0SXTRCE('KBXFCIN: CAN NOT READ CAL HEADER')
      KBXCAL=-1
      RETURN

      END
C
C--------------------------------------------------------------------

      INTEGER FUNCTION KBXOPT(CFUNC,IIN,IOUT)

      IMPLICIT NONE
      CHARACTER*4 CFUNC
      INTEGER IIN(*)
      INTEGER IOUT(*)

      ! external functions

      INTEGER LIT
      integer kbenhstr
      character*12 cfi

      INCLUDE 'areaparm.inc'

      INTEGER ITYPE
      INTEGER JOPT(NUMAREAOPTIONS)
      INTEGER CALFLG
      INTEGER CALARR(416)
      integer           enh_str_flag    ! =1 for new BRIT stretch
      character*8       enh_str_name    ! stretch table name
      COMMON/M0FCI/ITYPE,JOPT,CALFLG,CALARR,
     &              enh_str_flag, enh_str_name
C
      KBXOPT=0

      IF(CFUNC(1:4).EQ.'KEYS') THEN
         IF( IIN(4).LE.8 ) THEN
             IOUT(1)=4
             IOUT(2)=LIT('RAW ')
             IOUT(3)=LIT('RAD ')
             IOUT(4)=LIT('REFL')
             IOUT(5)=LIT('BRIT')
         ELSE
             IOUT(1)=4
             IOUT(2)=LIT('RAW ')
             IOUT(3)=LIT('RAD ')
             IOUT(4)=LIT('TEMP')
             IOUT(5)=LIT('BRIT')
         ENDIF
      ELSEIF(CFUNC(1:4).EQ.'INFO') THEN
         IF( IIN(4).LE.8 ) THEN
             IOUT(1)=4
             IOUT(2)=LIT('RAW ')
             IOUT(3)=LIT('RAD ')
             IOUT(4)=LIT('REFL')
             IOUT(5)=LIT('BRIT')
             IOUT(6)=LIT('    ')
             IOUT(7)=LIT('WM**')
             IOUT(8)=LIT('%   ')
             IOUT(9)=LIT('    ')
             IOUT(10)=1
             IOUT(11)=100
             IOUT(12)=100
             IOUT(13)=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('WM**')
             IOUT( 8)=LIT('K   ')
             IOUT( 9)=LIT('    ')
             IOUT(10)=1
             IOUT(11)=100
             IOUT(12)=100
             IOUT(13)=1
         ENDIF

      elseif (cfunc .eq. 'STR ') then
c
c --- ORIGinal or EXPanded stretch
c
           enh_str_flag = kbenhstr(iin, enh_str_name)
           call movcw('        ', iout)
           if(enh_str_flag .lt. 0) kbxopt = -4
           call movcw(enh_str_name, iout)

      ELSEIF(CFUNC(1:4).EQ.'CALB') THEN
         CALFLG = 1
         CALL MOVW(416,IIN,CALARR)
      ELSE
         KBXOPT=-1
      ENDIF

      RETURN
      END