mtclim42_vic.c

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

/*
  Note:  The following code is largely taken from mtclim4.2, a weather
         preprocessor developed by the NTSG group in the School of Forestry at
	 the University of Montana.  The code has been left intact as much as
	 possible, but small adaptations have been made to allow integration
	 with VIC.  Adaptations have been made so the functions, etc. are local
	 to this routine, which should allow easier updating.
	 The original comments have largely been left intact
	 typedefs and function definitions have been moved to mtclim42_vic.h
	 functions are called from mtclim42_wrapper.c
  Changes:  The only changes to the original MTCLIM code are 
            - cloudiness is returned for each day as well through the variable
	    tskc (calculated using equation 2.29 in Bras, R. L., "Hydrology, an
	    introduction to hydrologic science", Addison-Wesley, Reading,
	    Massachusetts, 1990).
	    - For each day of the year the fraction of the total daily potential
	    radiation that occurs during each radiation time step is also
	    returned through the variable tiny_radfract
	    - SRADDT has been changed to be 300 sec
            - 04-Jun-04 If data->s_srad becomes negative, we set it to 0.0.     TJB
	    - 15-Jun-04 If p->base_isoh and p->site_isoh are both small or 0, set
	      ratio = 1.  This handles the case in which annual precip for the grid
	      cell is 0, resulting in both p->base_isoh and p->site_isoh being 0,
	      and their ratio being undefined.					TJB
	    Changes are preceded by the comment * start vic_change *  and
	    followed by the comment * end vic_change *
  Author: Most of the code was written by Peter E. Thornton at the Univeristy of 
          Montana (see below).
	  Adapted by: Bart Nijssen
          Adaptation started on Sat Aug 21 using the mtclim4.2 code from
	  5/7/1999 
*/

/*
mtclim42.c

MTCLIM
VERSION 4.2  

Peter Thornton
NTSG, School of Forestry
University of Montana
5/7/99

***************************************
** Questions or comments? Contact... **
** Peter E. Thornton                 **
** NTSG, School of Forestry          **
** University of Montana             **
** Missoula, MT 59812                **
** email: [e-mail]         **
** phone: [phone]               **
***************************************

CHANGES FROM VERSION 4.1 TO VERSION 4.2
1) PUT THE SLOPE AND ASPECT CORRECTIONS TO RADIATION BACK IN. THEY
HAD BEEN REMOVED DURING THE DEVELOPMENT OF THE NEW RADIATION CODE.
This includes the estimation for diffuse radiation on sloping surfaces
during periods when the sun is above a flat horizon, but not providing
direct illumination to the slope.  Site east and west horizon corrections
have also been reintroduced.

CHANGES FROM VERSION 4.0 TO VERSION 4.1
1) ADDITIONAL REVISION OF RADIATION CODE, FOLLOWING SUBMISSION OF AG FOR MET
MANUSCRIPT DESCRIBING ANALYSIS OF SAMSON DATA.  The current code
follows exactly the algorithm detailed in
Thornton, P.E. and S.W. Running, 1999. An improved algorithm for estimating 
incident daily solar radiation from measurements of temperature, humidity,
and precipitation. Ag For Met 93:211-228.

changes from version 3.1 to version 4.0:
1) radiation code completely revamped, following analysis of samson
observations for the vemap2 project.
2) includes an iterative algorithm for estimating vapor pressure and
radiation
3) par output no longer an option, since the old par algorithm is suspect
4) 2-day tmin smoothing no longer an option
5) solar constant now calculated as an analytical function of yearday,
instead of the monthly array of values in earlier versions. removes the
discontinuities between months. 
6) boxcar function now uses the previous n days of data, making it a
"pulled boxcar" instead of a "centered boxcar"

Changes from version 2.0 to version 3.1
Modified to include the following improvements to the original MTCLIM logic:
1) Improved vapor pressure calculation
2) Improved transmissivity calculation
3) Improved daylength calculation

Some other differences between version 3.1 and previous versions of MTCLIM:
1) No english units option
2) No total or average radiation option
3) No threshold radiation option
5) No pre-rainy days correction for transmissivity
6) No radiation correction to air temperatures
7) No LAI corrections to air temperatures
8) Only one precipitation station allowed
9) Some parameters formerly in initialization file are now in mtclim_const.h
*/

/*
--------------------------------------------------------------------------

UNITS:
Temperatures         degrees C
Temp. lapse rates    degrees C / 1000 m
Precipitation        cm / day
Vapor pressure       Pa (also for Vapor Pressure Deficit, VPD)
Radiation            W/m2, average over daylight period
Daylength            s, sunrise to sunset, flat horizons 
Elevation            m
Latitude             decimal degrees
Aspect               decimal degrees
Slope                decimal degrees
E/W horizons         decimal degrees

--------------------------------------------------------------------------
*/
/*
Code History
------------
Original code written by R.R. Nemani
Updated  4/ 1/1989 by J.C. Coughlan
Updated  6/ 1/1989 by J.C. Coughlan
Updated 12/23/1989 by Joe Glassy
Updated  3/ 2/1990 by Raymond Hunt
Updated  2/ 4/1991 by Raymond Hunt
Updated  1/ 4/1993 by Raymond Hunt (version 2.1)
Updated  3/26/1997 by Peter Thornton (version 3.0)
Updated  5/16/1997 P.T. corrected error in calculation of tmean and tday,
spotted by Mike White
Updated  7/14/1997 P.T. (version 3.1) replaced an older version of the vapor
pressure correction with a newer algorithm. The algorithm implemented in this
version corresponds to that found in:
Kimball et al., in press. An improved method for estimating surface humidity
from daily minimum temperature. Ag. For. Met., 1997 (in press).
Updated 7/28/1997 P.T. Finalized MTCLIM version 3.1
Updated 9/2/97 P.E.T. cosmetic change in boxcar()
Updated 5/7/98 P.E.T. Version 4.0 (radiation and humidity algorithms)
Updated 8/1/98 P.E.T. Version 4.1 radiation code matches manuscript submitted
	to Ag For Met in July 98.
Updated 4/20/99 PET Version 4.2 added slope and aspect corrections to radiation.

--------------------------------------------------------------------------
*/

/*
Files
-----

Parameter initialization file
Input meteorological data file
Output meteorological data file  (*.mtcout)
*/

/*
Example initialization file:

---top of init file-------------------------------------------------------
sample               (this entire line written to outfiles)
other comments       (this entire line discarded)

IOFILES                    Keyword, don't edit this line 
test.mtcin                 Name of input meteorological data file
test                       Prefix for output file
                      
CONTROL                    Keyword, don't edit this line
3                          (int) Number of header lines in input file
365                        (int) Number of days of data in input file
0                          (int) Dewpoint temperature input? (0=NO, 1=YES)
1                          (int) Humidity output flag        (0=VPD, 1=VP)
1                          (int) Year field in input file?   (0=NO, 1=YES)

PARAMETERS                 Keyword, don't edit this line
500.0                      (double) Base station elevation, meters
50.0                       (double) Base station annual precip isohyet, cm
48.0                       (double) Site latitude, degrees (- for south)
1500.0                     (double) Site elevation, meters
15.0                       (double) Site slope, degrees
180.0                      (double) Site aspect, degrees (0=N,90=E,180=S,270=W)
75.0                       (double) Site annual precip isohyet, cm
2.0                        (double) Site east horizon, degrees
5.0                        (double) Site west horizon, degrees
-6.0                       (double) Lapse rate for max temperature, deg C/1000m
-3.0                       (double) Lapse rate for min temperature, deg C/1000m

END                        Keyword, don't edit this line
---end of init file-------------------------------------------------------

*.init FILE INFO
For all lines, except the header and comment lines, the  parameter value on the
left can be followed by comments on the right, as long as there is whitespace
separating the value on the left from the following comment. Blank lines can be
inserted or deleted, but all keyword lines and parameter lines must be
included.  The order and number of non-blank lines in this file is crucial. The
keywords are there as a rudimentary quality check on the format of the
initialization file, and they ensure that the proper number of lines are read.
They DON'T ensure that the parameters are in the proper order.

NOTE: The dashed lines at the top and bottom of the example file shown above
are NOT part of the file.
*/

/*
INPUT FILE INFO
All input temperatures are in degrees Celcius, and precipitation is in 
centimeters.
Input data file format:
<some number of header lines (can be zero)>
<year (optional)> <yearday> <Tmax> <Tmin> <Tdew (optional)> <precipitation>
.
.
.

Example input data file... without year field, and without dewpoint temperature
---start of input data file --------------
This is a header line, which is discarded
101 16.0 2.0 1.2
102 17.0 3.0  0.1
103 16.5 5.2  0.0
104 20.1 6.4  0.0
---end of input data file ----------------
*/

/*
OUTPUT FILE INFO
The output file is created by appending ".mtcout" to the output filename
prefix specified in the initialization file. Existing files are overwritten,
so be careful to rename files between runs if you want to save the results, 
and for safety, don't use ".mtcout" as the extension for the input data file.
*/

/*
******************************
Input and Output CONTROL FLAGS
******************************
There is a flag in the initialization file that controls the input of dewpoint
temperature information. If your input file does not contain dewpoint data,
set this flag to 0. Otherwise, if you have dewpoint temperature information 
in your input file, set this flag to 1, and be sure that the dewpoint
temperature field is between the tmin and prcp fields, as specified under
the heading "INPUT FILE INFO", above.

There is another flag in the initialization file that controls the output of
humidity information. When set to 0, humidity output is the vapor pressure
deficit from the saturated vapor pressure at the daylight average temperature
to the ambient vapor pressure, in Pa. When the flag is set to 1, the output is
simply the  ambient vapor pressure in Pa.  By using the ambient vapor pressure
(flag  set to 1), you can avoid possible errors related to a difference between
the temperature chosen for saturation calculations (in this case tday) and the
actual temperature at any given time of day.

Another flag controls the treatment of a year-field in the input and output
files. When the flag is set (1), it is assumed that the first field in the
input file contains the year, which can be any integer between -32000 and
32000 (approx). This field is simply copied line-for-line into the output
file, where it preceeds the standard yearday output field.

Example output data file... (using VPD output, no PAR, no year field)
---start of output data file -----------------------------------
MTCLIM v3.1 OUTPUT FILE : Mon Mar 24 14:50:51 1997
MTCLIM version 3.1 testing
  yday    Tmax    Tmin    Tday    prcp      VPD     srad  daylen
       (deg C) (deg C) (deg C)    (cm)     (Pa)  (W m-2)     (s)
   101   10.00   -0.50    7.39    1.80   439.51   379.95   47672
   102   11.00    1.10    7.67    0.15   387.27   364.95   47880
   103   10.50    2.80    6.84    0.00   243.78   302.15   48087
   104   14.10    3.40    9.29    0.00   390.67   390.22   48293
---end of output data file -------------------------------------

*/

/*********************
**                  **
**  START OF CODE   **
**                  **
*********************/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <vicNl.h>

#include <mtclim42_vic.h>   /* physical constants */

static char vcid[] = "$Id: mtclim42_vic.c,v 4.1.2.3 2004/06/15 21:10:43 tbohn Exp $";

/****************************
 **                         **
 **    START OF FUNCTION    **
 **       DEFINITIONS       **
 **                         **
 ****************************/

/* data_alloc() allocates space for input and output data arrays */
int data_alloc(const control_struct *ctrl, data_struct *data)
{
  int ok=1;
  int ndays;
  
  ndays = ctrl->ndays;
  
  if (ok && ctrl->inyear && !(data->year = (int*) malloc(ndays * sizeof(int)))) {
    printf("Error allocating for year array\n");
    ok=0;
  } 
  if (ok && !(data->yday = (int*) malloc(ndays * sizeof(int)))) {
    printf("Error allocating for yearday array\n");
    ok=0;
  } 
  if (ok && !(data->tmax = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for tmax array\n");
    ok=0;
  }
  if (ok && !(data->tmin = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for tmin array\n");
    ok=0;
  }
  if (ok && !(data->prcp = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for prcp array\n");
    ok=0;
  }
  if (ok && ctrl->indewpt && 
      !(data->tdew = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for input humidity array\n");
    ok=0;
  }
  if (ok && !(data->s_tmax = (double*) malloc(ndays * sizeof(double))))	{
    printf("Error allocating for site Tmax array\n");
    ok=0;
  }
  if (ok && !(data->s_tmin = (double*) malloc(ndays * sizeof(double))))	{
    printf("Error allocating for site Tmin array\n");
    ok=0;
  }
  if (ok && !(data->s_tday = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for site Tday array\n");
    ok=0;
  }
  if (ok && !(data->s_prcp = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for site prcp array\n");
    ok=0;
  }
  if (ok && !(data->s_hum = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for site VPD array\n");
    ok=0;
  }
  if (ok && !(data->s_srad = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for site radiation array\n");
    ok=0;
  }
  if (ok && !(data->s_dayl = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for site daylength array\n");
    ok=0;
  }
  /* start vic_change */
  if (ok && !(data->s_tskc = (double*) malloc(ndays * sizeof(double)))) {
    printf("Error allocating for site cloudiness array\n");
    ok=0;
  }
  /* end vic_change */
  return (!ok);
} /* end of data_alloc() */

/* calc_tair() calculates daily air temperatures */
int calc_tair(const control_struct *ctrl, const parameter_struct *p, 
	      data_struct *data)
{
  int ok=1;
  int i,ndays;
  double dz;
  double tmean, tmax, tmin;
  
  ndays = ctrl->ndays;
  /* calculate elevation difference in kilometers */
  dz = (p->site_elev - p->base_elev)/1000.0;
  
  /* apply lapse rate corrections to tmax and tmin */
  /* Since tmax lapse rate usually has a larger absolute value than tmin
     lapse rate, it is possible at high elevation sites for these corrections
     to result in tmin > tmax. Check for that occurrence and force
     tmin = corrected tmax - 0.5 deg C. */
  for (i=0 ; i<ndays ; i++) {
    /* lapse rate corrections */
    data->s_tmax[i] = tmax = data->tmax[i] + (dz * p->tmax_lr);
    data->s_tmin[i] = tmin = data->tmin[i] + (dz * p->tmin_lr);
    
    /* derived temperatures */
    tmean = (tmax + tmin)/2.0;
    data->s_tday[i] = ((tmax - tmean)*TDAYCOEF) + tmean;
  }