/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Procedures in this file are for the retrieval of surface temperature.
 They are: 
   call_surftemp_clear
   surftemp_clear

 Compilation:  See the file README.C for instructions on compiling this file.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ 

#include <assert.h>
#include <errno.h>
#include <stdio.h>
#include <math.h>
#include "globals.h"

long call_surftemp_clear(short argc, void *argv[]);
float surftemp_clear (float tb4, float tb5, float scanang, short hemis, 
                      short satellite, char surftype);


/******************** CALL CLEAR SKY SURFACE TEMPERATURE ********************/

long call_surftemp_clear(short argc, void *argv[])
/*--------------------------------------------------------------------------
  This function loops through the input arrays and calls 'surftemp_clear'
  (below) for each pixel.  It is called from an IDL program with the 
  following statement:
  
  ERROR = call_external(objfln,'call_surftemp_clear',ncols,nrows,tb4,tb5,$
                        scanang,hemis,satellite,surftype,cmask,badpix,tsurf)

  Input (as pointers to the specified type):
    ncols, nrows   : Dimensions of 'tsurf' and the arrays pointed to by 'array' 
                     [short integer].  While these arrays are 2D, they are
                     treated here as 1D.  If scalar values are passed in,
                     both of these values must be 1.  If a 1D array is passed
                     in, then 'nrows' must be 1.
    tb4,tb5        : AVHRR channels 4 and 5 brightness temperatures (K) [float]
    scanang        : Sensor scan angle (degrees) [float]
    hemis          : Hemisphere (0=southern,1=northern) [short integer] 
    satellite      : satellite number (7,9,11,12,14,15(K),16(L))  [short integer]
    surftype       : Surface type as LAND, OCEAN, ICE, SNOW (see 'globals.pro' for 
                     values).  LAND means snow-free land (vegetation cover), 
                     OCEAN means open water, ICE means sea ice (with or without 
                     snow), and SNOW means snow-covered land.
                     However, corrections are currently done for only land
                     and ice surfaces.  Any other surface type uses the ice
                     procedure. [char]
    cmask          : Cloud mask array, cloudy=CLOUDY, clear=CLEAR ('globals.h') 
                     [char]
    badpix         : Valid pixels array, where invalid=TRUE ('globals.h') [char]

  Output:
    tsurf          : Surface temperature (K) [float]
                     If badpix=TRUE, MISSING values are returned.  If the 
                     pixel is CLOUDY then the values that this variable
                     has on input is unchanged.

  Additional Variables (from globals.h):
    TRUE, FALSE    : true and false flags to be used for returning error 
                    condition

  Called by: surftemp_clear (an IDL procedure)

  Calls: surftemp_clear (the C function below)
----------------------------------------------------------------------------*/
{

/*  parameters passed in from IDL  */

  short  ncols, nrows, hemis, satellite;
  char   *cmask, *surftype, *badpix;
  float  *tb4, *tb5, *scanang, *tsurf;

/*  local parameters  */

  long i;
  long  ERROR;                 /* Error flag, returned to IDL */

/*  Check that correct number of parameters was passed  */

  if (argc != 11) {
    fprintf(stdout, "argc = %d\r\n",argc);
    fprintf(stdout, "ERROR >>> CALL_SURFTEMP_CLEAR: INCORRECT NUMBER OF ARGUMENTS\r\n");
    ERROR = TRUE;
    return ERROR;
  }

/*  Cast passed parameters to local vars  */

  ncols     = *(short *) argv[0];
  nrows     = *(short *) argv[1];
  tb4       = (float *) argv[2];
  tb5       = (float *) argv[3];
  scanang   = (float *) argv[4];
  hemis     = *(short *) argv[5];
  satellite = *(short *) argv[6];
  surftype  = (char *) argv[7];
  cmask     = (char *) argv[8];
  badpix    = (char *) argv[9];
  tsurf     = (float *) argv[10];

/*  Loop through each pixel in the arrays.  If it is a "bad" pixel, skip
    it and assign the fluxes a MISSING value.  */

  ERROR = FALSE;

  for (i=0; i<nrows*ncols; i++) {
    if (badpix[i] != (char) TRUE) {
      if (cmask[i] == (char) CLEAR) 
        tsurf[i] = surftemp_clear(tb4[i],tb5[i],scanang[i],hemis,satellite,
                                  surftype[i]);
    }
    else   
      tsurf[i] = MISSING;
  }

return ERROR;
}


/********************** CLEAR SKY SURFACE TEMPERATURE ***********************/

float surftemp_clear (float tb4, float tb5, float scanang, short hemis, 
                      short satellite, char surftype)
/*---------------------------------------------------------------------------- 
  Surface temperature retrieval for clear sky conditions.  
 
  This procedure operates on scalar values only.
  
  Input: 
     tb4       - Channel 4 brightness temperature (K) 
     tb5       - Channel 5 brightness temperature (K) 
     scanang   - scan angle (degrees) 
     hemis     - hemisphere (0=southern,1=northern) 
     satellite - satellite number (7,9,11,12,14,15(K),16(L)) 
     surftype  - A value representing the surface type (global variables
                 LAND, OCEAN, ICE, SNOW* see 'globals.h' for values).
                 LAND means snow-free land (vegetation cover), OCEAN means 
                 open water, ICE means sea ice (with or without snow), and
                 SNOW means snow-covered land.
                 However, corrections are currently done for only land
                 and ice surface.  Any other surface type uses the ice
                 procedure.
  
  Returned as value of function: 
     tsurf     - Clear sky surface temperature (K)
 
  Calls: None
 
  Called by: call_surftemp_clear  
 
  Algorithm Summary:
 
  The clear sky temperature of ice and snow (and currently open ocean) 
  is estimated as a function of the split-window brightness temperature 
  difference, the channel 4 temperature, and the sensor scan angle.  
  The algorithm uses coefficients for three temperature ranges and is
  based on simulations using LOWTRAN with Arctic Ocean and Antarctic 
  radiosonde data.  The form of the correction equation is:
 
       tsurf = a + b*Tb4 + c*(Tb4-Tb5) + d*(Tb4-Tb5)*(sec(scanang)-1)
 
  where there is a unique set of coefficients a, b, c, and d for the 
  temperature ranges: < 240K, 240-260K, and > 260K.
 
  For snow-free land surfaces a somewhat different equation is used:
 
       tsurf = a + b*Tb4 + c*Tb5 + e*emiss4 + f*emiss5
 
  where 'emiss4' and 'emiss5' are the emissivities of the surface in 
  channels 4 and 5.  The coefficients in this equation are also derived
  from a regression analysis of modelled radiances, but radiosonde profiles 
  from only Arctic land areas are used in the analysis.  The two emissivities
  are set to constant values below.  If other emissivities are desired,
  this file must be edited.
---------------------------------------------------------------------------*/
{
  short sat;
  short temprange;
  float tsurf;

/* SNOW-ICE COEFFICIENTS: 
   a,b,c,d (i,j,k)  where i=0,1 hemisphere: 0=Antarctic,1=Arctic 
                          j=0..2 temp range: 0 = <240, 1 = 240-260, 2 = >260
                          k=0..6 NOAA satellite: 7,9,11,12,14,15,16  */

  float a[2][3][7] =
  {{{-1.216191, -1.762823, -1.466105, -0.800189, -1.180768, -1.258905, -0.292169},
    {-6.400722, -8.083512, -7.100426, -4.823715, -5.647473, -6.381229, -3.650344},
    {-7.000356, -7.985411, -7.398464, -6.114496, -6.602028, -7.036723, -5.768286}},

   {{-3.824689, -5.482065, -4.655325, -2.798275, -3.846001, -4.034918, -1.453829},
    {-4.605042, -6.541145, -5.393341, -3.475960, -4.397758, -4.728985, -2.304536},
    {-4.415813, -5.254914, -4.769341, -4.121089, -4.498451, -4.633262, -3.676576}}};

  float b[2][3][7] =
  {{{1.004336, 1.007446, 1.005666, 1.002277, 1.003951, 1.004476, 1.000107},
    {1.025619, 1.032868, 1.028625, 1.019081, 1.022697, 1.025568, 1.013983},
    {1.027369, 1.031757, 1.029138, 1.023614, 1.025938, 1.027533, 1.021676}},

   {{1.014527, 1.021786, 1.018096, 1.010394, 1.014939, 1.015444, 1.004601},
    {1.017617, 1.025859, 1.020959, 1.013116, 1.017179, 1.018188, 1.008018},
    {1.016484, 1.020431, 1.018134, 1.015023, 1.017059, 1.017383, 1.012527}}};

  float c[2][3][7] =
  {{{1.365568, 0.477679, 1.092875, 1.729552, 1.616824, 1.430976, 1.506015},
    {0.981037, 0.600573, 0.857094, 1.138661, 0.919952, 0.973153, 1.181797},
    {1.079770, 0.921390, 1.035730, 1.174924, 0.963384, 1.085787, 1.215382}},

   {{2.228758, 1.995829, 2.196793, 2.100035, 1.904878, 2.222340, 1.946308},
    {1.795312, 1.647280, 1.763992, 1.681569, 1.411066, 1.757905, 1.637119},
    {1.666478, 1.635749, 1.664891, 1.669003, 1.402800, 1.659278, 1.690164}}};

  float d[2][3][7] =
  {{{-0.650609, -0.080109, -0.477556, -0.757764, -0.729174, -0.661893, -0.566436},
    { 0.562568,  1.158426,  0.766613,  0.383119,  0.552105,  0.607435, 0.194899},
    { 0.889363,  1.433513,  1.073912,  0.676141,  0.793529,  0.920915, 0.477659}},

   {{-1.294083, -1.183647, -1.268939, -1.027160, -0.879747, -1.228140, -0.907209},
    {-0.080299,  0.278680,  0.041162, -0.018824,  0.213224, -0.001810, -0.131031},
    { 0.684021,  1.147774,  0.847497,  0.547265,  0.690478,  0.733954, 0.347890}}};

/* SNOW-FREE LAND COEFFICIENTS: 
   al,bl,cl,dl,el (j,k) where
         j=0..2 temp range: 0 = <240, 1 = 240-260, 2 = >260
         k=0..6 NOAA satellite: 7,9,11,12,14,15,16  */

  float al[3][7] = 
   {{26.0309, 23.0055, 24.5757, 29.1836, 27.8690, 26.0307, 31.7447}, 
    {32.1194, 29.3755, 30.9222, 34.8680, 33.5189, 32.0778, 37.6640},
    {44.4224, 41.5469, 43.0879, 46.9049, 45.9367, 44.1398, 49.0214}};

  float bl[3][7] = 
   {{4.0147, 4.4368, 4.2369, 3.4836, 3.2520, 3.9893, 3.2120},
    {3.5683, 3.6499, 3.5992, 3.5896, 3.2796, 3.5843, 3.5496},
    {3.6507, 3.7915, 3.7034, 3.6529, 3.3803, 3.6687, 3.5950}};

  float cl[3][7] = 
   {{-2.9919, -3.4103, -3.2127, -2.4606, -2.2260, -2.9655, -2.1935},
    {-2.5444, -2.6167, -2.5714, -2.5732, -2.2587, -2.5600, -2.5422},
    {-2.6387, -2.7710, -2.6874, -2.6470, -2.3725, -2.6554, -2.5950}};

  float dl[3][7] = 
   {{-165.0710, -181.5454, -173.8222, -144.4215, -133.0235, -164.3357, -133.4750},
    {-164.3970, -167.8258, -165.6568, -166.2492, -150.3413, -165.3660, -164.6871},
    {-181.1707, -188.3021, -183.7980, -181.5388, -166.7875, -182.2059, -178.4917}};

  float el[3][7] = 
   {{133.5685, 152.2116, 143.4666, 109.7186, 98.9686, 132.5978, 97.1853}, 
    {126.3626, 130.4036, 127.9483, 127.2383, 111.6787, 127.3099, 124.9140},
    {133.4351, 141.5502, 136.5114, 132.7192, 118.6037, 134.4543, 128.7903}};

  float emiss4 = 0.985;
  float emiss5 = 0.975;
/*--------------------------------------------------------------------------*/

  if ((satellite != 7) && (satellite != 9) && (satellite != 11) &&
      (satellite != 12) && (satellite != 14) && (satellite != 15)  && 
      (satellite != 16)) {
    fprintf (stderr, 
      "ERROR >>> SURFTEMP_CLEAR: UNSUPPORTED SATELLITE NUMBER, exiting\n");
    exit (-1);
  }

/* Convert satellite number to an index and use the channel 4 temperature
   to get an index for the correct temperature range.  */

  if (satellite == 7)
    sat = 0;
  else if (satellite == 9)
    sat = 1;
  else if (satellite == 11)
    sat = 2;
  else if (satellite == 12)
    sat = 3;
  else if (satellite == 14)
    sat = 4;
  else if (satellite == 15)
    sat = 5;
  else if (satellite == 16)
    sat = 6;

  if (tb4 < 240.)
    temprange = 0;
  else if ((tb4 >= 240.) && (tb4 < 260.))
    temprange = 1;
  else if (tb4 >= 260)
    temprange = 2;

  /* Operate on the appropriate surface type. */
  
   if (surftype == (char)LAND) {
    
       /* SNOW-FREE LAND */

       tsurf = al[temprange][sat] + bl[temprange][sat]*tb4 + 
         cl[temprange][sat]*tb5 + dl[temprange][sat]*emiss4 + 
         el[temprange][sat]*emiss5;
    } 
    else {
       
       /* OCEAN, SNOW/ICE ALGORITHM */
       
       tsurf = a[hemis][temprange][sat] + b[hemis][temprange][sat] * tb4 +
         c[hemis][temprange][sat] * (tb4 - tb5)
         + d[hemis][temprange][sat] * (tb4 - tb5) * 
         ((1.0 / cos (scanang * DEGRAD)) - 1.0);
    }

   return tsurf;
}