/*
 * NAME: 
 *
 *    diurnal_gentemann
 *
 * FUNCTION:
 *
 *    C implementation of Gentemann (2003) diurnal warming parameterization
 *
 * TYPE/LANGUAGE: 
 *
 *    C module
 *
 * DESCRIPTION:  
 *
 *    Contains routines to estimate diurnal warming based on the parameterization
 *    of Gentemann et al. (2003) on a global grid, given an input of wind speed on
 *    the grid and date/time information.  Gentemann (2003) parameterization is 
 *    based on theoretical calculation of top-of-atmosphere insolation.
 *
 *    Note that only the main routine is considered external for use outside this module
 *    and others are internal (defined as static)
 *
 * ROUTINES TO BE USED EXTERNALLY:
 *
 *        diurnal_gentemann - estimates diurnal warming for global grid for given windspeed and date/time
 *
 * ROUTINES USED INTERNALLY:
 *
 *        leapYear          - return 1 if year is leap year
 *        eqTimeCorr        - return correction for difference between mean & true local solar time
 *
 * HISTORY:
 * 
 *   Original version by Andy Harris on July 5th, 2012
 *
 */

/*  Standard C library includes  */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/*  Include header file for this module  */

#include "diurnal_gentemann.h"

/*  Prototypes of routines for internal use only  */

static float eqTimeCorr( const int dayOfYear, const float totDays );
static int leapYear( const int year );

/*  Define constants  */

#define	SPRING_EQUINOX	    80		/*  Day of year for equinox  */
#define	TRUE		     1
#define	FALSE		     0
#define	SUCCESS		     0
#define	FAIL		    -1
#define	TROPIC_LAT	    23.5	/*  Latitude of tropics  */
#define	S0		  1367.0	/*  Solar constant  */
#define	Q_THRESH	   132.0	/*  When insolation is less than this, dSST is set to zero  */
#define WIND_THRESHOLD	    50.0	/*  Consider winds above this to be bad  */
#define	BADVAL		  -999.0	/*  If wind is deemed "bad" then set warming to this value  */

#define	DEBUG

/*  Define macros  */

#define	MIN(a,b)	(((a)> (b))? (b):(a) )
#define	MAX(a,b)	(((a)> (b))? (a):(b) )
#define	ABS(a)		(((a)>(-a))? (a):(-a))

/*
 * NAME: 
 *
 *    diurnal_gentemann
 *
 * FUNCTION:
 *
 *    C implementation of Gentemann (2003) diurnal warming parameterization
 *
 * TYPE/LANGUAGE: 
 *
 *    C module
 *
 * DESCRIPTION:  
 *
 *    Contains routines to estimate diurnal warming based on the parameterization
 *    of Gentemann et al. (2003) on a global grid, given an input of wind speed on
 *    the grid and date/time information.  Gentemann (2003) parameterization is 
 *    based on theoretical calculation of top-of-atmosphere insolation.  The
 *    reference for the model is:
 *
 *    Gentemann, C.L, C.J. Donlon, A. Stuart-Menteth, F.J. Wentz, 
 *    "Diurnal signals in satellite sea surface temperature measurements", 
 *    GRL, 30(3), 1140-1143, 2003
 *
 *    Note that only the main routine is considered external for use outside this module
 *    and others are internal (defined as static)
 *
 * ROUTINES TO BE USED EXTERNALLY:
 *
 *        diurnal_gentemann - estimates diurnal warming for global grid for given windspeed and date/time
 *
 * ROUTINES USED INTERNALLY:
 *
 *        leapYear          - return 1 if year is leap year
 *        eqTimeCorr        - return correction for difference between mean & true local solar time
 *
 * FILES NEEDED: 
 *    
 * SUBROUTINE/FUNCTION CALLS: 
 *
 *    leapYear
 *    eqTimeCorr
 *
 * CALLING SEQUENCE:  
 *
 *    success = diurnal_gentemann( year, month, day, hour, mins,
 *					wind, nlats, nlons, dsst )
 *
 * INPUTS:  
 *    
 *   year            - year
 *   month           - month of year
 *   day             - day of month
 *   hour            - hour of day (UTC)
 *   mins            - minutes of hour
 *   wind            - wind speed array (N.B. speed implies magnitude, i.e. values should be >= 0)
 *   nlats           - number of latitudes in array
 *   nlons           - number of longitudes in array
 *
 * OUTPUTS:  
 * 
 *   dsst            - output delta-SSTs predicted by model
 *
 * RETURNS:
 *
 *    Returns a status
 *
 * SYSTEM CALLS:
 *
 *    sin
 *    cos
 *    tan
 *    acos
 *
 * HISTORY:
 * 
 *   Original version by Andy Harris on July 5th, 2012
 *
 */

int diurnal_gentemann( const int year, const int month, const int day, const int hour, const int mins, 
			const float *wind, const int nlats, const int nlons, float *dsst )
{
    int		ilon = 0, ilat = 0, i = 0, indx = 0, dayOfYear = 0, nDays[12], iResult = 0;
    float	dQ = 0., pi, d2r, phi = 0., delta = 0., H = 0., totDays = 0.;
    float	Qsol = 0., w = 0., lat = 0., lon = 0., *ft_wt = NULL;
    float	t = 0., w_t = 0., omega = 0.;
    float	tanPhi_tanDelta = 0., latIncr = 0., latOff = 0., lonIncr = 0., lonOff = 0.;

/*  5-term Fourier series describing time-dependence of delta-T  */

    float	fT[] = {6.814, -6.837, -8.427, 1.447, 4.274, -0.407, -0.851, 0.457, -0.555, -0.101, 0.375};

    float	aET = 0.;

/*  Precision-limited parameter definitions  */

    pi    = 2.0*asin(1.0);		/*  Pi  */
    d2r   = pi/180.0;			/*  degrees to radians conversion factor  */
    omega = 2.0*pi/24.0;		/*  Radians per hour  */

    nDays[0]  = 0;
    nDays[1]  = 31;
    nDays[2]  = 59;
    nDays[3]  = 90;
    nDays[4]  = 120;
    nDays[5]  = 151;
    nDays[6]  = 181;
    nDays[7]  = 212;
    nDays[8]  = 243;
    nDays[9]  = 273;
    nDays[10] = 304;
    nDays[11] = 334;

/*  Adjust for leap year  */

    totDays = 365.0;

    if(TRUE == (iResult = leapYear(year)))
	{
	for(i=2; i<12; ++i)
	    {
	    nDays[i] += 1;
	    }

	totDays += 1.0;
	}

    dayOfYear = nDays[month - 1] + day;

/*  For now, anticipate 360x181 fields, i.e. 1 degree, but can adjust for (e.g.) 1/2 degree  */

    lonIncr = 360.0/((float)(nlons));
    latIncr = lonIncr;

    lonOff = lonIncr/2.0;	/*  E.g. first cell will be -179.5, last will be +179.5  */
    latOff = -90.0;		/*  Assumes it's (e.g.) 181 cells in latitude  */

/*  Zero the output array.  Diurnal warming is only calculated if Qsol is > threshold & deemed zero elsewhere  */

    for(i=0; i<nlats*nlons; ++i);
	{
	*(dsst + i) = 0.0;
	}

/*  Allocate memory for Gentemann FT-derived "shape"  */

    if((ft_wt = (float *)malloc(nlons*sizeof(float))) == NULL)
	{
	fprintf(stderr, "diurnal_gentemann: Failed to allocate memory for ft_wt\n");
	return(FAIL);
	}

/*  Latitude of sub-solar "point" (or "track") for the day in question  */

    delta = TROPIC_LAT*d2r*sin(2.0*pi*(dayOfYear - SPRING_EQUINOX)/totDays);

/*  Equation of time correction for the day in question  */

    aET = eqTimeCorr(dayOfYear, totDays);	/*  Units are hours  */

/*  
 *  Loop to pre-calculate warming "shape"
 */

    for(ilon=0; ilon<nlons; ++ilon)
	{
	lon = ilon*lonIncr + lonOff;

/*  Convert from 0 -> 360 to -180 -> +180  */

	if(lon > 180.0)
	    {
	    lon -= 360.0;
	    }

/*  True local solar time for this longitude  */

	t = hour + mins/60.0 + lon/15.0 - aET;

	if(t > 24.0)
	    {
	    t -= 24.0;
	    }

	if(t < 0.0)
	    {
	    t += 24.0;
	    }

/*  Value of diurnal warming shape at this local time  */

	w_t = omega*t;		/*  Angular time (radians)  */

	*(ft_wt + ilon) = (fT[0]
			+ fT[1]*cos(w_t)     + fT[2]*sin(w_t)
			+ fT[3]*cos(2.0*w_t) + fT[4]*sin(2.0*w_t)
			+ fT[5]*cos(3.0*w_t) + fT[6]*sin(3.0*w_t)
			+ fT[7]*cos(4.0*w_t) + fT[8]*sin(4.0*w_t)
			+ fT[9]*cos(5.0*w_t) + fT[10]*sin(5.0*w_t))
			*0.001;
	}


/*
 *  Loop through latitudes first.  This is more efficient in terms of
 *  calculating theoretical insolation and memory mapping.  Also, since
 *  the use of theoretical insolation is a major approximation, it is
 *  deemed to make very little difference if the part of the globe is
 *  from the day before (or after).  This was also done in the Gentemann
 *  et al. (2003) paper.
 */

    for(ilat=0; ilat<nlats; ++ilat)
	{
	lat = ilat*latIncr + latOff;

	/*  Calculate length of 1/2 day  */

	phi = d2r*lat;							/*  Latitude in radians  */

	tanPhi_tanDelta = tan(phi)*tan(delta);

	H = pi;		/*  If Sun doesn't set, this is the correct value for the half-day  */

	if((ABS(tanPhi_tanDelta)) < 1.0)
	    {
	    H = acos(-1.0*tanPhi_tanDelta);	/*  This is the "usual" case, and always true if -66.5 < lat < +66.5  */
	    }

	Qsol = 0.0;

	if(tanPhi_tanDelta > -1.0)		/*  Will be true unless Sun does not rise at this location  */
	    {
	    Qsol = (S0/pi)*(sin(phi)*sin(delta)*H + cos(phi)*cos(delta)*sin(H));
	    }

	if(Qsol > Q_THRESH)			/*  If it's at or below this threshold then warming will be zero  */
	    {

	    dQ = Qsol - Q_THRESH;
/*
 *  Loop through longitudes for this latitude
 */

	    for(ilon=0; ilon<nlons; ++ilon)
		{

/*  Since value of shape is precalculated for each longitude, just need to convolve with dQ & wind for this location  */

		indx = ilat*nlons + ilon;

		w = *(wind + indx);

		if((w <= WIND_THRESHOLD) && (w >= 0.))	/*  Catch erroneous winds  */
		    {		    
		    *(dsst + indx) = *(ft_wt + ilon)*(dQ - 9.632e-4*dQ*dQ)*exp(-0.53*w);
		    }
		else
		    {
		    *(dsst + indx) = BADVAL;
		    }
		}
	    }
	}

    free(ft_wt);

    return(SUCCESS);
}
    
/*
 * NAME: 
 *
 *    eqTimeCorr
 *
 * FUNCTION:
 *
 *    C implementation of correction for equation of time
 *
 * TYPE/LANGUAGE: 
 *
 *    C module
 *
 * DESCRIPTION:  
 *
 *    Contains routine to calculate equation of time correction, i.e.
 *    the difference between mean & true local solar time.
 *
 * ROUTINES TO BE USED EXTERNALLY:
 *
 *        eqTimeCorr - correction to mean local solar time 
 *		       N.B. intended to be local to diurnal_gentemann
 *
 * ROUTINES USED INTERNALLY:
 *
 * FILES NEEDED: 
 *    
 * SUBROUTINE/FUNCTION CALLS: 
 *
 * CALLING SEQUENCE:  
 *
 *    result = eqTimeCorr( dayOfYear, totDays )
 *
 * INPUTS:  
 *    
 *   dayOfYear       - day of year
 *   totDays         - total number of days in year
 *
 * OUTPUTS:  
 * 
 * RETURNS:
 *
 *    Returns time correction in hours
 *
 * SYSTEM CALLS:
 *
 *    sin
 *    cos
 *    tan
 *    acos
 *
 * HISTORY:
 * 
 *   Original version by Andy Harris on July 5th, 2012
 *
 */

static float eqTimeCorr(const int dayOfYear, const float totDays)
{
  float	pi = 0., r2d = 0., aD = 0., aET = 0.;
  float e[] = {0.000075, 0.001868, 0.032077, 0.014615, 0.040849};

/*  Initialize precision-limited parameter definitions  */

    pi    = 2.0*asin(1.0);		/*  Pi  */
    r2d   = 180.0/pi;			/*  radians to degrees conversion factor  */

    aD = 2.0*pi*dayOfYear/totDays;	/*  Angular position of earth around sun  */

/*  Equation of time correction...  */

    aET = e[0] + e[1]*cos(aD) - e[2]*sin(aD) - e[3]*cos(2.0*aD) - e[4]*sin(2.0*aD);

    aET = aET*r2d/15.;		/*  Convert from angle (radians) to time (hours) N.B. 15 degrees per hour  */

    return(aET);

}

/*
 * NAME: 
 *
 *    leapYear
 *
 * FUNCTION:
 *
 *    C implementation for leap year determination
 *
 * TYPE/LANGUAGE: 
 *
 *    C module
 *
 * DESCRIPTION:  
 *
 *    Contains routine to calculate whether or not the input year
 *    is a leap year
 *
 * ROUTINES TO BE USED EXTERNALLY:
 *
 *        leapYear - leap year determination 
 *		     N.B. intended to be local to diurnal_gentemann
 *
 * ROUTINES USED INTERNALLY:
 *
 * FILES NEEDED: 
 *    
 * SUBROUTINE/FUNCTION CALLS: 
 *
 * CALLING SEQUENCE:  
 *
 *    result = eqTimeCorr( year )
 *
 * INPUTS:  
 *    
 *   year       - year
 *
 * OUTPUTS:  
 * 
 * RETURNS:
 *
 *    Returns TRUE if leap year
 *
 * SYSTEM CALLS:
 *
 * HISTORY:
 * 
 *   Original version by Andy Harris on July 5th, 2012
 *
 */

static int leapYear( const int year )
{
    if( (0 == year%4 && 0 != year%100) || (0 == year%400) )
	{
	return(TRUE);	/*  If it's a leap year  */
	}
    else
	{
	return(FALSE);	/*  If it's not a leap year  */
	}
}