mtclim42_vic.c

超强的大尺度水文模拟工具

       Fig 6.14, p. 122. */
    pdif = -1.25*t_final + 1.25;
    if (pdif > 1.0) 
      pdif = 1.0;
    if (pdif < 0.0) 
      pdif = 0.0;
    
    /* estimate fraction of radiation that is direct, on an
       instantaneous basis */
    pdir = 1.0 - pdif;
    
    /* the daily total radiation is estimated as the sum of the
       following two components:
       1. The direct radiation arriving during the part of
       the day when there is direct beam on the slope.
       2. The diffuse radiation arriving over the entire daylength
       (when sun is above ideal horizon).
    */
    
    /* component 1 */
    srad1 = slope_potrad[yday] * t_final * pdir;
    
    /* component 2 */
    srad2 = flat_potrad[yday] * t_final * pdif * sky_prop; 
    
    /* save daily radiation and daylength */
    data->s_srad[i] = srad1 + srad2;
    data->s_dayl[i] = daylength[yday];
  }
  
  /* free local array memory */
  free(dtr);
  free(sm_dtr);
  
  return (!ok);
} /* end of calc_srad_humidity() */

/* without Tdew input data, an iterative estimation of shortwave radiation
   and humidity is required */
/* start vic_change */
int calc_srad_humidity_iterative(const control_struct *ctrl,
				 const parameter_struct *p, data_struct *data,
				 double *tiny_radfract)
  /* end vic_change */
{
  int ok=1;
  int i,j,ndays;
  int start_yday,end_yday,isloop;
  int ami,yday;
  double ttmax0[366];
  double flat_potrad[366];
  double slope_potrad[366];
  double daylength[366];
  double *dtr, *sm_dtr;
  double *parray, *window, *t_fmax, *tdew;
  double sum_prcp,effann_prcp;
  double tmax,tmin;
  double t1,t2;
  double pratio;
  double lat,coslat,sinlat,dt,h,dh;
  double cosslp,sinslp,cosasp,sinasp;
  double bsg1,bsg2,bsg3;
  double decl,cosdecl,sindecl,cosegeom,sinegeom,coshss,hss;
  double sc,dir_beam_topa;
  double sum_flat_potrad,sum_slope_potrad,sum_trans;
  double cosh,sinh;
  double cza,cbsa,coszeh,coszwh;
  double dir_flat_topa,am;
  double pva,t_tmax,b;
  double tmink,pet,ratio,ratio2,ratio3,tdewk;
  double pvs,vpd;
  double trans1,trans2;
  double t_final,pdif,pdir,srad1,srad2; 
  double pa;
  double sky_prop;
  double avg_horizon, slope_excess;
  double horizon_scalar, slope_scalar;
  
  /* optical airmass by degrees */
  double optam[21] = {2.90,3.05,3.21,3.39,3.69,3.82,4.07,4.37,4.72,5.12,5.60,
		      6.18,6.88,7.77,8.90,10.39,12.44,15.36,19.79,26.96,30.00};
  
  /* these are the radiation algorithm parameters, as described in
     AgForMet manuscript */
  double tbase = 0.870;
  double abase = -6.1e-5;
  double c = 1.5;
  double b0 = 0.031;
  double b1 = 0.201;
  double b2 = 0.185;
  /* start vic_change */
  long endtiny;
  long starttiny;
  long tinystep;
  long tinystepspday;
  /* end vic_change */
  
  /* number of simulation days */
  ndays = ctrl->ndays;
  
  /* local array memory allocation */
  /* allocate space for DTR and smoothed DTR arrays */
  if (!(dtr = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for DTR array\n");
    ok=0;
  }
  if (!(sm_dtr = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for smoothed DTR array\n");
    ok=0;
  }
  /* allocate space for effective annual precip array */
  if (!(parray = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for effective annual precip array\n");
    ok=0;
  }
  /* allocate space for the prcp totaling array */
  if (!(window = (double*) malloc((ndays+90)*sizeof(double)))) {
    printf("Error allocating for prcp totaling array\n");
    ok = 0;
  }
  /* allocate space for t_fmax */
  if (!(t_fmax = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for p_tt_max array\n");
      ok=0;
  }
  /* allocate space for Tdew array */
  if (!(tdew = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for Tdew array\n");
    ok=0;
  }
  
  /* calculate diurnal temperature range for transmittance calculations */
  for (i=0 ; i<ndays ; i++) {
    tmax = data->tmax[i];
    tmin = data->tmin[i];
    if (tmax < tmin) 
      tmax = tmin;
    dtr[i] = tmax-tmin;
  }
  
  /* smooth dtr array: After Bristow and Campbell, 1984 */
  if (ndays >= 30) { /* use 30-day antecedent smoothing window */
    if (pulled_boxcar(dtr, sm_dtr, ndays, 30, 0)) {
      printf("Error in boxcar smoothing, calc_srad_humidity()\n");
      ok=0;
    }
  }
  else { /* smoothing window width = ndays */
    if (pulled_boxcar(dtr, sm_dtr, ndays, ndays, 0)) {
      printf("Error in boxcar smoothing, calc_srad_humidity()\n");
      ok=0;
    }
  }
  
  /* Generate the effective annual precip, based on a 3-month
     moving-window. Requires some special case handling for the
     beginning of the record and for short records. */
  /* check if there are at least 90 days in this input file, if not,
     use a simple total scaled to effective annual precip */
  if (ndays < 90) {
    sum_prcp = 0.0;
    for (i=0 ; i<ndays ; i++) {
      sum_prcp += data->s_prcp[i];
    }
    effann_prcp = (sum_prcp/(double)ndays) * 365.25;
    /* if the effective annual precip for this period
       is less than 8 cm, set the effective annual precip to 8 cm
       to reflect an arid condition, while avoiding possible
       division-by-zero errors and very large ratios (PET/Pann) */
    if (effann_prcp < 8.0) 
      effann_prcp = 8.0;
    for (i=0 ; i<ndays ; i++) {
      parray[i] = effann_prcp;
	}
  }
  else {
    /* Check if the yeardays at beginning and the end of this input file
       match up. If so, use parts of the three months at the end
       of the input file to generate effective annual precip for
       the first 3-months. Otherwise, duplicate the first 90 days
       of the record. */
    start_yday = data->yday[0];
    end_yday = data->yday[ndays-1];
    if (start_yday != 1) {
      isloop = (end_yday == start_yday-1) ? 1 : 0;
    }
    else {
      isloop = (end_yday == 365 || end_yday == 366) ? 1 : 0;
    }
    
    /* fill the first 90 days of window */
    for (i=0 ; i<90 ; i++) {
      if (isloop) 
	window[i] = data->s_prcp[ndays-90+i];
      else 
	window[i] = data->s_prcp[i];
    }
    /* fill the rest of the window array */
    for (i=0 ; i<ndays ; i++) {
      window[i+90] = data->s_prcp[i];
    }
    
    /* for each day, calculate the effective annual precip from 
       scaled 90-day total */
    for (i=0 ; i<ndays ; i++)	{
      sum_prcp = 0.0;
      for (j=0 ; j<90 ; j++) {
	sum_prcp += window[i+j];
      }
      sum_prcp = (sum_prcp/90.0) * 365.25;
      /* if the effective annual precip for this 90-day period
	 is less than 8 cm, set the effective annual precip to 8 cm
	 to reflect an arid condition, while avoiding possible
	 division-by-zero errors and very large ratios (PET/Pann) */
      parray[i] = (sum_prcp < 8.0) ? 8.0 : sum_prcp;
    }
  } /* end if ndays >= 90 */	
  
  /*****************************************
   *                                       *
   * start of the main radiation algorithm *
   *                                       *
   *****************************************/
  
  /* before starting the iterative algorithm between humidity and 
     radiation, calculate all the variables that don't depend on 
     humidity so they only get done once. */
  
  /* STEP (1) calculate pressure ratio (site/reference) = f(elevation) */
  t1 = 1.0 - (0.0065 * p->site_elev)/288.15;
  t2 = 9.80665 / (0.0065 * (8.3143/28.9644e-3));
  pratio = pow(t1,t2);
  
  /* STEP (2) correct initial transmittance for elevation */ 
  trans1 = pow(tbase,pratio);
  
  /* STEP (3) build 366-day array of ttmax0, potential rad, and daylength */
  
  /* precalculate the transcendentals */
  lat = p->site_lat;
  /* check for (+/-) 90 degrees latitude, throws off daylength calc */
  lat *= RADPERDEG;
  if (lat > 1.5707) 
    lat = 1.5707;
  if (lat < -1.5707) 
    lat = -1.5707;
  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));
  
  /* sub-daily time and angular increment information */
  dt = SRADDT;                /* set timestep */ 
  dh = dt / SECPERRAD;        /* calculate hour-angle step */
  /* start vic_change */
  tinystepspday = 86400L/(long)SRADDT;
  /* end vic_change */

  /* begin loop through yeardays */
  for (i=0 ; i<365 ; i++) {
    /* calculate cos and sin of declination */
    decl = MINDECL * cos(((double)i + DAYSOFF) * RADPERDAY);
    cosdecl = cos(decl);
    sindecl = sin(decl);
    
    /* 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;
    
    /* calculate daylength as a function of lat and decl */
    cosegeom = coslat * cosdecl;
    sinegeom = sinlat * sindecl;
    coshss = -(sinegeom) / cosegeom;
    if (coshss < -1.0) 
      coshss = -1.0;  /* 24-hr daylight */
    if (coshss > 1.0) 
      coshss = 1.0;    /* 0-hr daylight */
    hss = acos(coshss);                /* hour angle at sunset (radians) */
    /* daylength (seconds) */
    daylength[i] = 2.0 * hss * SECPERRAD;
    
    /* start vic_change */
    if (daylength[i] > 86400)
      daylength[i] = 86400;
    starttiny = (i * 86400L)/(long)SRADDT;
    endtiny = ((i+1L) * 86400L)/(long)SRADDT - 1L;
    /* end vic_change */
    
    /* solar constant as a function of yearday (W/m^2) */
    sc = 1368.0 + 45.5*sin((2.0*PI*(double)i/365.25) + 1.7);
    /* extraterrestrial radiation perpendicular to beam, total over
       the timestep (J) */
    dir_beam_topa = sc * dt;
    
    sum_trans = 0.0;
    sum_flat_potrad = 0.0;
    sum_slope_potrad = 0.0;
    
    /* begin sub-daily hour-angle loop, from -hss to hss */
    for (h=-hss ; h<hss ; h+=dh) {
      /* 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;
      
      /* check if sun is above a flat horizon */
      if (cza > 0.0) {
	/* when sun is above the ideal (flat) horizon, do all the
	   flat-surface calculations to determine daily total
	   transmittance, and save flat-surface potential radiation
	   for later calculations of diffuse radiation */
	
	/* potential radiation for this time period, flat surface,
	   top of atmosphere */
	dir_flat_topa = dir_beam_topa * cza;
	
	/* determine optical air mass */
	am = 1.0/(cza + 0.0000001);
	if (am > 2.9) {
	  ami = (int)(acos(cza)/RADPERDEG) - 69;
	  if (ami < 0) 
	    ami = 0;
	  if (ami > 20) 
	    ami = 20;
	  am = optam[ami];
	}
	
	/* correct instantaneous transmittance for this optical
	   air mass */
	trans2 = pow(trans1,am);
	
	/* instantaneous transmittance is weighted by potential
	   radiation for flat surface at top of atmosphere to get
	   daily total transmittance */
	sum_trans += trans2 * dir_flat_topa;
	
	/* 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;
	}
	
      } /* end if sun above ideal horizon */
      else dir_flat_topa = -1;

      /* start vic_change */
      tinystep = (long) ((i * 86400L + 12L * 3600L + h * SECPERRAD)/SRADDT);
      if (tinystep < starttiny)
	tinystep = starttiny;
      if (tinystep > endtiny)
	tinystep = endtiny;
      if (dir_flat_topa > 0)
	tiny_radfract[tinystep] = dir_flat_topa;
      else
	tiny_radfract[tinystep] = 0;
      /* end vic_change */
      
    } /* end of sub-daily hour-angle loop */
    
    /* start vic_change */
    if (daylength[i] && sum_flat_potrad > 0) {
      for (tinystep = (i * 86400L)/(long)SRADDT, j = 0; j < tinystepspday;  
	   tinystep++, j++)
	tiny_radfract[tinystep] /= sum_flat_potrad;
    }
    /* end vic_change */

    /* 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 {