Add a slope-aspect adjustment to solar coming out of metsim

[bartnijssen]

Slope and aspect were calculated in the original version of MTCLIM and were included in the old solar code of MTCLIM that was part of VIC 4.2. The best place to get started may be to look in the VIC 4.2.d repo, especially the following files:

• https://github.com/UW-Hydro/VIC/blob/support/VIC.4.2.d/src/mtclim_vic.c from which the code snippets below are taken
• https://github.com/UW-Hydro/VIC/blob/support/VIC.4.2.d/src/mtclim_constants_vic.h
• https://github.com/UW-Hydro/VIC/blob/support/VIC.4.2.d/src/mtclim_parameters_vic.h

Note that the indentation in the following is not that meaningful as the lines of codes are executed in different code blocks (different time intervals)

    coslat = cos(lat);
    sinlat = sin(lat);
    cosslp = cos(p->site_slp * RADPERDEG);
    sinslp = sin(p->site_slp * RADPERDEG);
    cosasp = cos(p->site_asp * RADPERDEG);
    sinasp = sin(p->site_asp * RADPERDEG);
    /* cosine of zenith angle for east and west horizons */
    coszeh = cos(1.570796 - (p->site_ehoriz * RADPERDEG));
    coszwh = cos(1.570796 - (p->site_whoriz * RADPERDEG));
...
/* do some precalculations for beam-slope geometry (bsg) */
    bsg1 = -sinslp * sinasp * cosdecl;
    bsg2 = (-cosasp * sinslp * sinlat + cosslp * coslat) * cosdecl;
    bsg3 = (cosasp * sinslp * coslat + cosslp * sinlat) * sindecl;
...
/* potential radiation for this time period, flat surface, top of atmosphere */
    dir_flat_topa = dir_beam_topa * cza;
...
 /* precalculate cos and sin of hour angle */
    cosh = cos(h);
    sinh = sin(h);
      
/* calculate cosine of solar zenith angle */
    cza = cosegeom * cosh + sinegeom;
...      
/* calculate cosine of beam-slope angle */
    cbsa = sinh * bsg1 + cosh * bsg2 + bsg3;

/* keep track of total potential radiation on a flat surface for ideal horizons */
    sum_flat_potrad += dir_flat_topa;
...	
/* keep track of whether this time step contributes to component 1 (direct on slope) */
    if ((h<0.0 && cza>coszeh && cbsa>0.0) ||
        (h>=0.0 && cza>coszwh && cbsa>0.0)) {
	  
 /* sun between east and west horizons, and direct on slope. this period contributes to component 1 */
    sum_slope_potrad += dir_beam_topa * cbsa;
...
/* calculate maximum daily total transmittance and daylight average flux density for a flat surface and the slope */
    if (daylength[i]) {
        ttmax0[i] = sum_trans / sum_flat_potrad;
        flat_potrad[i] = sum_flat_potrad / daylength[i];
        slope_potrad[i] = sum_slope_potrad / daylength[i];
    }
    else {
        ttmax0[i] = 0.0;
        flat_potrad[i] = 0.0;
        slope_potrad[i] = 0.0;
    }