mtclim42_vic.c

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

      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;
  
  /* b parameter, and t_fmax not varying with Tdew, so these can be
     calculated once, outside the iteration between radiation and humidity
     estimates. Requires storing t_fmax in an array. */
  for (i=0 ; i<ndays ; i++) {	
    /* b parameter from 30-day average of DTR */
    b = b0 + b1 * exp(-b2 * sm_dtr[i]);
    
    /* proportion of daily maximum transmittance */
    t_fmax[i] = 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[i] *= 0.75;
    
  }
  
  /* As a first approximation, calculate radiation assuming
     that Tdew = Tmin */
  for (i=0 ; i<ndays ; i++) {
    yday = data->yday[i]-1;
    /* start vic_change */
    /* pva = 610.7 * exp(17.38 * data->tmin[i] / (239.0 + data->tmin[i])); */
    pva = 1000. * svp(data->tmin[i]);
    /* end vic_change */
    t_tmax = ttmax0[yday] + abase * pva;
    
    /* final daily total transmittance */
    t_final = t_tmax * t_fmax[i];

    /* 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)
       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];
  }
  
  /* Estimate Tdew using the initial estimate of radiation for PET */
  /* estimate air pressure at site */
  pa = atm_pres(p->site_elev);
  for (i=0 ; i<ndays ; i++) {
    tmink = data->s_tmin[i] + 273.15;
    pet = calc_pet(data->s_srad[i],data->s_tday[i],pa,data->s_dayl[i]);
    
    /* calculate ratio (PET/effann_prcp) and correct the dewpoint */
    ratio = pet/parray[i];
    ratio2 = ratio*ratio;
    ratio3 = ratio2*ratio;
    tdewk = tmink*(-0.127 + 1.121*(1.003 - 1.444*ratio + 12.312*ratio2 
				   - 32.766*ratio3) + 0.0006*(dtr[i]));
    tdew[i] = tdewk - 273.15;
  }
  
  /* Revise estimate of radiation using new Tdew */
  for (i=0 ; i<ndays ; i++) {
    yday = data->yday[i]-1;
    /* start vic_change */
    /* pva = 610.7 * exp(17.38 * tdew[i] / (239.0 + tdew[i])); */
    pva = 1000. * svp(tdew[i]);
    /* end vic_change */
    t_tmax = ttmax0[yday] + abase * pva;
    
    /* final daily total transmittance */
    t_final = t_tmax * t_fmax[i];

    /* start vic_change */
    data->s_tskc[i] = sqrt((1.-t_fmax[i])/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)
	 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 */
    data->s_srad[i] = srad1 + srad2;

    /* start vic_change */
    /* under some conditions near the arctic circle in winter, s_srad < 0 */
    if (data->s_srad[i] < 0.0)
      data->s_srad[i] = 0.0;
    /* end vic_change */

  }
  
  /* Revise estimate of Tdew using new radiation */
  for (i=0 ; i<ndays ; i++) {
    tmink = data->s_tmin[i] + 273.15;
    pet = calc_pet(data->s_srad[i],data->s_tday[i],pa,data->s_dayl[i]);
    
    /* calculate ratio (PET/effann_prcp) and correct the dewpoint */
    ratio = pet/parray[i];
    ratio2 = ratio*ratio;
    ratio3 = ratio2*ratio;
    tdewk = tmink*(-0.127 + 1.121*(1.003 - 1.444*ratio + 12.312*ratio2 
				   - 32.766*ratio3) + 0.0006*(dtr[i]));
    tdew[i] = tdewk - 273.15;
  }
  
  /* now calculate vapor pressure from tdew */
  for (i=0 ; i<ndays ; i++) {
    /* start vic_change */
    /* pva = 610.7 * exp(17.38 * tdew[i] / (239.0 + tdew[i])); */
    pva = 1000. * svp(tdew[i]);
    /* end vic_change */
    if (ctrl->outhum) {
      /* output humidity as vapor pressure (Pa) */
      data->s_hum[i] = pva;
    }
    else {
      /* output humidity as vapor pressure deficit (Pa) */
      /* calculate saturated VP 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 */
      vpd = pvs - pva;
      if (vpd < 0.0) 
	vpd = 0.0;
      data->s_hum[i] = vpd;
    }
  } /* end for i = ndays loop */
  
  /* free local array memory */
  free(dtr);
  free(sm_dtr);
  free(parray);
  free(window);
  free(t_fmax);
  free(tdew);
  
  return (!ok);
} /* end of calc_srad_humidity_iterative() */
  
/* data_free frees the memory previously allocated by data_alloc() */
int data_free(const control_struct *ctrl, data_struct *data)
{
  int ok=1;
  if (ctrl->inyear) free(data->year);
  free(data->yday);
  free(data->tmax);
  free(data->tmin);
  free(data->prcp);
  if (ctrl->indewpt) free(data->tdew);
  free(data->s_tmax);
  free(data->s_tmin);
  free(data->s_tday);
  free(data->s_prcp);
  free(data->s_hum);
  free(data->s_srad);
  free(data->s_dayl);
  /* start vic_change */
  free(data->s_tskc);
  /* end vic_change */
  return (!ok);
}

/* calc_pet() calculates the potential evapotranspiration for aridity 
   corrections in calc_vpd(), according to Kimball et al., 1997 */
double calc_pet(double rad, double ta, double pa, double dayl)
{
  /* input parameters and units :
     double rad      (W/m2)  daylight average incident shortwave radiation
     double ta       (deg C) daylight average air temperature
     double pa       (Pa)    air pressure
     double dayl     (s)     daylength 
  */
  
  double rnet;       /* (W m-2) absorbed shortwave radiation avail. for ET */
  double lhvap;      /* (J kg-1) latent heat of vaporization of water */ 
  double gamma;      /* (Pa K-1) psychrometer parameter */
  double dt = 0.2;   /* offset for saturation vapor pressure calculation */
  double t1, t2;     /* (deg C) air temperatures */
  double pvs1, pvs2; /* (Pa)   saturated vapor pressures */
  double pet;        /* (kg m-2 day-1) potential evapotranspiration */
  double s;          /* (Pa K-1) slope of saturated vapor pressure curve */
  
  /* calculate absorbed radiation, assuming albedo = 0.2  and ground
     heat flux = 10% of absorbed radiation during daylight */
  rnet = rad * 0.72;
  
  /* calculate latent heat of vaporization as a function of ta */
  lhvap = 2.5023e6 - 2430.54 * ta;
  
  /* calculate the psychrometer parameter: gamma = (cp pa)/(lhvap epsilon)
     where:
     cp       (J/kg K)   specific heat of air
     epsilon  (unitless) ratio of molecular weights of water and air
  */
  gamma = CP * pa / (lhvap * EPS);
  
  /* estimate the slope of the saturation vapor pressure curve at ta */
  /* temperature offsets for slope estimate */
  t1 = ta+dt;
  t2 = ta-dt;
  
  /* calculate saturation vapor pressures at t1 and t2, using formula from 
     Abbott, P.F., and R.C. Tabony, 1985. The estimation of humidity parameters.
     Meteorol. Mag., 114:49-56.
  */
  /* start vic_change */
  /* pvs1 = 610.7 * exp(17.38 * t1 / (239.0 + t1)); */
  pvs1 = 1000. * svp(t1);
  /* pvs2 = 610.7 * exp(17.38 * t2 / (239.0 + t2)); */
  pvs2 = 1000. * svp(t2);
  /* end vic_change */
  
  /* calculate slope of pvs vs. T curve near ta */
  s = (pvs1-pvs2) / (t1-t2);
  
  /* calculate PET using Priestly-Taylor approximation, with coefficient
     set at 1.26. Units of result are kg/m^2/day, equivalent to mm water/day */
  pet = (1.26 * (s/(s+gamma)) * rnet * dayl)/lhvap;
  
  /* return a value in centimeters/day, because this value is used in a ratio
     to annual total precip, and precip units are centimeters */
  return (pet/10.0);
}    

/* atm_pres() calculates the atmospheric pressure as a function of elevation */
double atm_pres(double elev)
{
  /* daily atmospheric pressure (Pa) as a function of elevation (m) */
  /* From the discussion on atmospheric statics in:
     Iribane, J.V., and W.L. Godson, 1981. Atmospheric Thermodynamics, 2nd
     Edition. D. Reidel Publishing Company, Dordrecht, The Netherlands.
     (p. 168)
  */
  
  int ok=1;
  double t1,t2;
  double pa;
  
  t1 = 1.0 - (LR_STD * elev)/T_STD;
  t2 = G_STD / (LR_STD * (R / MA));
  pa = P_STD * pow(t1,t2);
  
  return(pa);
}

/* pulled_boxcar() calculates a moving average of antecedent values in an
   array, using either a ramped (w_flag=1) or a flat (w_flag=0) weighting */	
int pulled_boxcar(double *input,double *output,int n,int w,int w_flag)
{
  int ok=1;
  int tail,i,j;
  double *wt;
  double total,sum_wt;
  
  if (w > n) {
    printf("Boxcar longer than array...\n");
    printf("Resize boxcar and try again\n");
    ok=0;
  }
  
  if (ok && !(wt = (double*) malloc(w * sizeof(double)))) {
    printf("Allocation error in boxcar()\n");
    ok=0;
  }
  
  if (ok) {
    /* when w_flag != 0, use linear ramp to weight tails,
       otherwise use constant weight */
    sum_wt = 0.0;
    if (w_flag) {
      for (i=0 ; i<w ; i++) {
	wt[i] = (double)(i+1);
	sum_wt += wt[i];
      }
    }
    else {
      for (i=0 ; i<w ; i++) { 	
	wt[i] = 1.0;
	sum_wt += wt[i];
      }
    }
    
    /* fill the output array, starting with the point where a full
       boxcar can be calculated */
    for (i=w-1 ; i<n ; i++) {
      total = 0.0;
      for (j=0 ; j<w ; j++) {
	total += input[i-w+j+1] * wt[j];
      }
      output[i] = total/sum_wt;
    }
    
    /* fill the first w elements of the output array with the value from
       the first full boxcar */
    for (i=0 ; i<w-1 ; i++) {
      output[i] = output[w-1];
    }
    
    free(wt);
    
  } /* end if ok */
  
  return (!ok);
}
/* end of pulled_boxcar() */