mtclim42_vic.c

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

  return (!ok);
}
/* end of calc_tair() */

/* calc_prcp() calculates daily total precipitation */
int calc_prcp(const control_struct *ctrl, const parameter_struct *p, 
	      data_struct *data)
{
  int ok=1;
  int i,ndays;
  double ratio;
  
  ndays = ctrl->ndays;

  /* start vic_change */
  if ( p->site_isoh < 1e-10 && p->base_isoh < 1e-10 ) {
    /* If base_isoh and site_isoh are both small or 0, set the ratio to 1.
       This handles case in which annual precip is 0, resulting in base_isoh
       and site_isoh being 0 and their ratio being undefined. */
    ratio = 1.;
  }
  else if (p->base_isoh == 0) {
    vicerror("Error in calc_prcp(): base_isoh == 0 and site_isoh/base_isoh is undefined.");
  }
  else {
    ratio = p->site_isoh / p->base_isoh;
  }
  /* end vic_change */

  if (ok) {
    for (i=0 ; i<ndays ; i++) {
      data->s_prcp[i] = data->prcp[i] * ratio;
    }
  }
  
  return (!ok);
}
/* end of calc_prcp() */

/* when dewpoint temperature observations are available, radiation and
   humidity can be estimated directly */
/* start vic_change */
int calc_srad_humidity(const control_struct *ctrl, const parameter_struct *p, 
		       data_struct *data, double *tiny_radfract)
/* end vic_change */
{
  int ok=1;
  int i,j,ndays;
  double pva,pvs,vpd;
  int ami,yday;
  double ttmax0[366];
  double flat_potrad[366];
  double slope_potrad[366];
  double daylength[366];
  double *dtr, *sm_dtr;
  double tmax,tmin;
  double t1,t2;
  double pratio;
  double lat,coslat,sinlat,dt,dh,h;
  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 trans1,trans2;
  double t_tmax,b,t_fmax;
  double t_final,pdif,pdir,srad1,srad2; 
  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;
  
  /* calculate humidity from Tdew observations */
  for (i=0 ; i<ndays ; i++) {
    /* convert dewpoint to vapor pressure */
    /* start vic_change */
    /* pva = 610.7 * exp(17.38 * data->tdew[i] / (239.0 + data->tdew[i])); */
    pva = 1000. * svp(data->tdew[i]);
    /* end vic_change */
    if (ctrl->outhum) {
      /* output humidity as vapor pressure */
      data->s_hum[i] = pva;
    }
    else {
      /* output humidity as vapor pressure deficit */
      /* calculate saturation vapor pressure at tday */
      /* start vic_change */
      /* pvs = 610.7 * exp(17.38 * data->s_tday[i] / (239.0 + data->s_tday[i])); 
       */
      pvs = 1000. * svp(data->s_tday[i]);
      /* end vic_change */
      /* calculate vpd */
      vpd = pvs-pva;
      if (vpd < 0.0) 
	vpd = 0.0;
      data->s_hum[i] = vpd;
    }
  }
  
  /* estimate radiation using Tdew observations */
  /* 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;
  }
  
  /* 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 using a 30-day antecedent smoothing window */
  if (ndays >= 30) {
    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;
    }
  }
  
  /*****************************************
   *                                       *
   * start of the main radiation algorithm *
   *                                       *
   *****************************************/
  
  /* 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 {
      ttmax0[i] = 0.0;
      flat_potrad[i] = 0.0;
      slope_potrad[i] = 0.0;
    }
  } /* end of i=365 days loop */
  
  /* force yearday 366 = yearday 365 */
  ttmax0[365] = ttmax0[364];
  flat_potrad[365] = flat_potrad[364];
  slope_potrad[365] = slope_potrad[364];
  daylength[365] = daylength[364];

  /* start vic_change */
  for (tinystep = 365L*86400L/(long)SRADDT, j = 0 ; j < tinystepspday;
       tinystep++, j++)
    tiny_radfract[tinystep] = tiny_radfract[tinystep-tinystepspday];
  /* end vic_change */

  /* STEP (4)  calculate the sky proportion for diffuse radiation */
  /* uses the product of spherical cap defined by average horizon angle
     and the great-circle truncation of a hemisphere. the interaction of
     these two terms is not exact, but pretty close. this factor does not
     vary by yearday. */
  avg_horizon = (p->site_ehoriz + p->site_whoriz)/2.0;
  horizon_scalar = 1.0 - sin(avg_horizon * RADPERDEG);
  if (p->site_slp > avg_horizon) 
    slope_excess = p->site_slp - avg_horizon;
  else 
    slope_excess = 0.0;
  if (2.0*avg_horizon > 180.0) 
    slope_scalar = 0.0;
  else {
    slope_scalar = 1.0 - (slope_excess/(180.0 - 2.0*avg_horizon));
    if (slope_scalar < 0.0) slope_scalar = 0.0;
  }
  sky_prop = horizon_scalar * slope_scalar;
  
  /* STEP (5)  final calculation of daily total radiation */
  for (i=0 ; i<ndays ; i++) {
    /* correct this day's maximum transmittance for vapor pressure */
    yday = data->yday[i]-1;
    /* start vic_change */
    /* pva = 610.7 * exp(17.38 * data->tdew[i] / (239.0 + data->tdew[i])); */
    pva = 1000. * svp(data->tdew[i]);
    /* end vic_change */
    t_tmax = ttmax0[yday] + abase * pva;
    
    /* b parameter from 30-day average of DTR */
    b = b0 + b1 * exp(-b2 * sm_dtr[i]);
    
    /* proportion of daily maximum transmittance */
    t_fmax = 1.0 - 0.9 * exp(-b * pow(dtr[i],c));
    
    /* correct for precipitation if this is a rain day */
    if (data->prcp[i]) 
      t_fmax *= 0.75;
    
    /* final daily total transmittance */
    t_final = t_tmax * t_fmax;

    /* start vic_change */
    data->s_tskc[i] = sqrt((1.-t_fmax)/0.65);
    /* end vic_change */

    /* estimate fraction of radiation that is diffuse, on an
       instantaneous basis, from relationship with daily total
       transmittance in Jones (Plants and Microclimate, 1992)
       Fig 2.8, p. 25, and Gates (Biophysical Ecology, 1980)