From 77c3b355c2101788f903068ebff324a930ccb42a Mon Sep 17 00:00:00 2001 From: Alomir Date: Fri, 6 Jun 2025 09:06:58 -0400 Subject: [PATCH 01/18] Removes code related to multi-site support --- Makefile | 4 +- docs/parameters.md | 84 +-- src/common/modelParams.c | 234 ++++++++ src/common/modelParams.h | 92 ++++ src/common/spatialParams.c | 776 --------------------------- src/common/spatialParams.h | 264 ---------- src/sipnet/events.c | 136 ++--- src/sipnet/events.h | 15 +- src/sipnet/frontend.c | 27 +- src/sipnet/sipnet.c | 1025 ++++++++++++++++-------------------- src/sipnet/sipnet.h | 78 ++- src/sipnet/sipnet.in | 4 +- 12 files changed, 920 insertions(+), 1819 deletions(-) create mode 100644 src/common/modelParams.c create mode 100644 src/common/modelParams.h delete mode 100644 src/common/spatialParams.c delete mode 100644 src/common/spatialParams.h diff --git a/Makefile b/Makefile index 02eeb7b60..90d5daf20 100755 --- a/Makefile +++ b/Makefile @@ -1,13 +1,13 @@ CC=gcc LD=gcc AR=ar -rs -CFLAGS=-Wall -Werror -g -Isrc +CFLAGS=-Wall -Werror -g -Isrc -Wno-c2x-extensions LIBLINKS=-lm LIB_DIR=./libs LDFLAGS=-L$(LIB_DIR) # Main executables -COMMON_CFILES:=util.c namelistInput.c spatialParams.c +COMMON_CFILES:=util.c namelistInput.c modelParams.c COMMON_CFILES:=$(addprefix src/common/, $(COMMON_CFILES)) COMMON_OFILES=$(COMMON_CFILES:.c=.o) diff --git a/docs/parameters.md b/docs/parameters.md index c44f46b1f..9764c0443 100644 --- a/docs/parameters.md +++ b/docs/parameters.md @@ -88,11 +88,13 @@ Run-time parameters can change from one run to the next, or when the model is st | 4 | $C_{\text{soil},0}$ | soilInit | Initial soil carbon | $\text{g C} \cdot \text{m}^{-2} \text{ ground area}$ | | | 5 | $W_{\text{litter},0}$ | litterWFracInit | | unitless | fraction of litterWHC | | 6 | $W_{\text{soil},0}$ | soilWFracInit | | unitless | fraction of soilWHC | -| | $N_{\text{soil},0}$ | | Initial soil organic nitrogen content | g N m$^{-2}$ | | -| | ${CH_4}_{\text{soil},0}$ | | Initial methane concentration in the soil | g C m$^{-2}$ | | -| | ${N_2O}_{\text{soil},0}$ | | Nitrous oxide concentration in the soil | g N m$^{-2}$ | | -| | $f_{\text{fine root},0}$ | fineRootFrac | Fraction of `plantWoodInit` allocated to initial fine root carbon pool | | | -| | $f_{\text{coarse root},0}$ | coarseRootFrac | Fraction of `plantWoodInit` allocated to initial coarse root carbon pool | | | +| 7 | $N_{\text{org, litter},0}$ | | Initial litter organic nitrogen content | g N m$^{-2}$ | | +| 8 | $N_{\text{org, soil},0}$ | | Initial soil organic nitrogen content | g N m$^{-2}$ | | +| 9 | $N_{\text{min},0}$ | | Initial soil mineral nitrogen content | g N m$^{-2}$ | +| 10 | ${CH_4}_{\text{soil},0}$ | | Initial methane concentration in the soil | g C m$^{-2}$ | | +| 11 | ${N_2O}_{\text{soil},0}$ | | Nitrous oxide concentration in the soil | g N m$^{-2}$ | | +| 12 | $f_{\text{fine root},0}$ | fineRootFrac | Fraction of `plantWoodInit` allocated to initial fine root carbon pool | | | +| 13 | $f_{\text{coarse root},0}$ | coarseRootFrac | Fraction of `plantWoodInit` allocated to initial coarse root carbon pool | | | ### Autotrophic respiration parameters @@ -389,25 +391,25 @@ aMaxFrac 0.85 0 0.66 0.86 0.005 ### Climate -For each step of the model, for each location, the following inputs are needed. These are provided in a file named `.clim` with the following columns: +For each step of the model, the following inputs are needed. These are provided in a file named `.clim` with the following columns: | col | parameter | description | units | notes | -| -- | ----------- | --------------------------------------------------- | ---------------------------------------- | ----------------------------------------------------------------------------------------------------------- | -| 1 | loc | spatial location index | | maps to param-spatial file, can be 0 for a single site, as when used by PEcAn | -| 2 | year | year of start of this timestep | | integer, e.g. 2010| -| 3 | day | day of start of this timestep | | integer where 1 = Jan 1| -| 4 | time | time of start of this timestep | hours after midnight | e.g. noon = 12.0, midnight = 0.0, can be a fraction | -| 5 | length | length of this timestep | days | variable-length timesteps allowed, typically not used | -| 6 | tair | avg. air temp for this time step| degrees Celsius | | -| 7 | tsoil | average soil temperature for this time step | degrees Celsius| can be estimated from Tair | -| 8 | par | average photosynthetically active radiation (PAR) for this time step | $\text{Einsteins} \cdot m^{-2} \text{ground area} \cdot \text{time step}^{-1}$ | input is in Einsteins \* m^-2 ground area, summed over entire time step | -| 9 | precip | total precip. for this time step| cm | input is in mm; water equivilant - either rain or snow | -| 10 | vpd | average vapor pressure deficit | kPa | input is in Pa, can be calculated from air temperature and relative humidity.| -| 11 | vpdSoil | average vapor pressure deficit between soil and air | kPa | input is in Pa ; differs from vpd in that saturation vapor pressure is calculated using Tsoil rather than Tair | -| 12 | vPress | average vapor pressure in canopy airspace | kPa | input is in Pa | -| 13 | wspd | avg. wind speed | m/s | | -| 14 | soilWetness | fractional soil wetness | unitless (0-1) | $f_\text{WHC}$; Used if `MODEL_WATER=0`; if `MODEL_WATER=1`, soil wetness is simulated| - +|-----| ----------- | --------------------------------------------------- | ---------------------------------------- | ----------------------------------------------------------------------------------------------------------- | +| 1 | year | year of start of this timestep | | integer, e.g. 2010| +| 2 | day | day of start of this timestep | | integer where 1 = Jan 1| +| 3 | time | time of start of this timestep | hours after midnight | e.g. noon = 12.0, midnight = 0.0, can be a fraction | +| 4 | length | length of this timestep | days | variable-length timesteps allowed, typically not used | +| 5 | tair | avg. air temp for this time step| degrees Celsius | | +| 6 | tsoil | average soil temperature for this time step | degrees Celsius| can be estimated from Tair | +| 7 | par | average photosynthetically active radiation (PAR) for this time step | $\text{Einsteins} \cdot m^{-2} \text{ground area} \cdot \text{time step}^{-1}$ | input is in Einsteins \* m^-2 ground area, summed over entire time step | +| 8 | precip | total precip. for this time step| cm | input is in mm; water equivilant - either rain or snow | +| 9 | vpd | average vapor pressure deficit | kPa | input is in Pa, can be calculated from air temperature and relative humidity.| +| 10 | vpdSoil | average vapor pressure deficit between soil and air | kPa | input is in Pa ; differs from vpd in that saturation vapor pressure is calculated using Tsoil rather than Tair | +| 11 | vPress | average vapor pressure in canopy airspace | kPa | input is in Pa | +| 12 | wspd | avg. wind speed | m/s | | +| 13 | soilWetness | fractional soil wetness | unitless (0-1) | $f_\text{WHC}$; Used if `MODEL_WATER=0`; if `MODEL_WATER=1`, soil wetness is simulated| + +Note: An older format for this file included location as the first column. #### Example `sipnet.clim` file: Column names are not used, but are: diff --git a/src/common/modelParams.c b/src/common/modelParams.c new file mode 100644 index 000000000..9e9dca857 --- /dev/null +++ b/src/common/modelParams.c @@ -0,0 +1,234 @@ +#include +#include +#include +#include +#include "modelParams.h" +#include "util.h" + +// Private/helper functions: not defined in modelParams.h + +// return 1 if all maxParams parameters have been initialized, 0 if not +int allParamsInitialized(ModelParams *ModelParams) { + if (ModelParams->numParams == ModelParams->maxParams) { + // all parameters have been initialized + return 1; + } + + return 0; +} + +// Set array[0..length-1] to all be equal to value +void setAll(double *array, int length, double value) { + int i; + + for (i = 0; i < length; i++) { + array[i] = value; + } +} + +// Checks to make sure that all parameters with isRequired=true have been read +// If not, kills program +// Writes out names of all parameters that weren't read (even if not required, +// as a warning message) +void checkAllRead(ModelParams *ModelParams) { + int i; + int okay; + OneModelParam *param; + + okay = 1; // so far so good + for (i = 0; i < ModelParams->numParams; i++) { + param = &(ModelParams->params[i]); + if (param->value == NULL) { + if (param->isRequired) { // should have been read but wasn't! + okay = 0; + printf("ERROR: Didn't read required parameter %s\n", param->name); + } else { + printf("WARNING: Didn't read parameter %s (not flagged as required, so " + "continuing)\n", + param->name); + } + } + } + + if (!okay) { + exit(1); + } +} + +/*************************************************/ + +// Public functions: defined in ModelParams.h + +// allocate space for a new ModelParams structure, return a pointer to it +ModelParams *newModelParams(int maxParams) { + ModelParams *modelParams; + + modelParams = (ModelParams *)malloc(sizeof(ModelParams)); + modelParams->maxParams = maxParams; + modelParams->numParams = 0; // for now, anyway + modelParams->numParamsRead = 0; // for now, anyway + + // allocate space for vectors within this structure: + modelParams->params = + (OneModelParam *)malloc(maxParams * sizeof(OneModelParam)); + modelParams->readIndices = (int *)malloc(maxParams * sizeof(int)); // the + // biggest + // it + // could + // be + + return modelParams; +} + +void initializeOneModelParam(ModelParams *modelParams, char *name, + double *externalLoc, int isRequired) { + int paramIndex; // index of next uninitialized parameter + OneModelParam *param; // a pointer to the next uninitialized parameter + + if (allParamsInitialized(modelParams)) { // all parameters have been + // initialized! + printf("Error trying to initialize %s: have already initialized all %d " + "parameters\n", + name, modelParams->maxParams); + printf("Check value of maxParams passed into newModelParams function\n"); + exit(1); + } + + // otherwise, get the index of the next uninitialized parameter + paramIndex = modelParams->numParams; + // and set param to point to it for easier access + param = &(modelParams->params[paramIndex]); + + strcpy(param->name, name); + param->value = externalLoc; + param->isRequired = isRequired; + + modelParams->numParams++; +} + +void checkParamFormat(char *line, const char *sep) { + // strtok modifies string, so we need a copy + char lineCopy[256]; + strcpy(lineCopy, line); + int numParams = 1; + char *par = strtok(lineCopy, sep); + while (par != NULL) { + par = strtok(NULL, sep); + ++numParams; + } + if (numParams > 2) { + printf("WARNING: extra columns in .param file are being ignored\n"); + } +} + +void readModelParams(ModelParams *modelParams, FILE *paramFile) { + const char *SEPARATORS = " \t\n\r"; // characters that can separate values in + // parameter files + const char *COMMENT_CHARS = "!"; // comment characters (ignore everything + // after this on a line) + + char line[256]; + char pName[MODEL_PARAM_MAXNAME]; // parameter name + int paramIndex; + OneModelParam *param; // a pointer to a single parameter, for easier access + char strValue[32]; // before we know whether value is a number or "*" + double value; + char *errc; + int isComment; + + // Check for old-style (spatial param) format on first line containing params + int formatChecked = 0; + + while (fgets(line, sizeof(line), paramFile) != NULL) { // while not EOF or + // error + // remove trailing comments: + isComment = stripComment(line, COMMENT_CHARS); + + if (!isComment) { // if this isn't just a comment line or blank line + // Check for old spatial-param format and warn if appropriate + if (!formatChecked) { + checkParamFormat(line, SEPARATORS); + formatChecked = 1; + } + + // tokenize line: + strcpy(pName, strtok(line, SEPARATORS)); // copy first token into pName + + strcpy(strValue, strtok(NULL, SEPARATORS)); + value = strtod(strValue, &errc); + paramIndex = locateParam(modelParams, pName); + + if (paramIndex == -1) { // not found + printf("WARNING: Ignoring parameter %s: this parameter wasn't " + "initialized in the code\n", + pName); + } else if (valueSet(modelParams, paramIndex)) { + printf("Error reading parameter file: read %s, but this parameter has " + "already been set\n", + pName); + exit(1); + } else { // otherwise, we're good to go + // set param to point to the appropriate parameter, for easier access + param = &(modelParams->params[paramIndex]); + *(param->value) = value; + param->isRead = 1; + // mark this as the next parameter read from file: + modelParams->readIndices[modelParams->numParamsRead] = paramIndex; + modelParams->numParamsRead++; + } // else (no errors in reading this line from parameter file) + } // if !isComment + } // while not EOF or error + + // check for error in reading: + if (ferror(paramFile)) { + printf("Error reading file in readModelParams\n"); + printf("ferror = %d\n", ferror(paramFile)); + exit(1); + } + + checkAllRead(modelParams); // terminate program if some required parameters + // weren't read +} // readModelParams + +// Return numParams, the actual number of parameters that have been +// initialized with initializeOneModelParam +int getnumParams(ModelParams *modelParams) { return modelParams->numParams; } + +// Return number of parameters that have been read in from file +int getNumParamsRead(ModelParams *modelParams) { + return modelParams->numParamsRead; +} + +// Find parameter with given name in the parameters vector +// If found, return index in vector, otherwise return -1 +int locateParam(ModelParams *modelParams, char *name) { + int i; + int found; + + i = 0; + found = 0; + while (i < modelParams->numParams && !found) { + if (strcmpIgnoreCase(name, modelParams->params[i].name) == 0) { + found = 1; + } + i++; + } + + if (found) { + return (i - 1); + } + + return -1; +} + +// Return 1 if parameter i has had its value set, 0 otherwise +int valueSet(ModelParams *modelParams, int i) { + return modelParams->params[i].isRead; +} + +void deleteModelParams(ModelParams *modelParams) { + // Free space used by ModelParams structure itself: + free(modelParams->params); + free(modelParams->readIndices); + free(modelParams); +} diff --git a/src/common/modelParams.h b/src/common/modelParams.h new file mode 100644 index 000000000..8fad9e5f1 --- /dev/null +++ b/src/common/modelParams.h @@ -0,0 +1,92 @@ +#ifndef MODEL_PARAMS_H +#define MODEL_PARAMS_H + +#include + +#define MODEL_PARAM_MAXNAME 64 + +// struct to hold a single param +typedef struct OneModelParamStruct { + char name[MODEL_PARAM_MAXNAME]; // name of parameter + int isRequired; // 0 or 1; 1 indicates that we should terminate run if this + // parameter isn't specified in input file + double *value; // a pointer to this param's (external location) value + int isRead; // whether this param has been read +} OneModelParam; + +typedef struct ModelParamsStruct { + OneModelParam *params; // vector of parameters, dynamically-allocated + // (stored in order in which they were initialized) + int maxParams; // maximum number of parameters + int numParams; // actual number of parameters (tracks # that have been + // initialized with initializeOneSpatialParam) + int numParamsRead; // tracks number of parameters that have been read in from + // file + int *readIndices; // indices (in parameters vector) of the parameters read + // from file, in the order in which they were read (this + // will allow us to output parameters in the same order + // later) +} ModelParams; + +// allocate space for a new spatialParams structure, return a pointer to it +ModelParams *newModelParams(int maxParameters); + +/* Initialize next parameter: + Set name of parameter equal to "name", + and set parameter's value pointer to point to the location given by + "externalLoc" parameter (required to be non-NULL) + Also set parameter's isRequired value: if true, then we'll terminate + execution if this parameter is not read from file +*/ +void initializeOneModelParam(ModelParams *params, char *name, + double *externalLoc, int isRequired); + +/* Read all parameters from file, setting all values in Params structure + appropriately + + Structure of paramFile: + name value + + The order of the parameters in paramFile does not matter. + + ! is a comment character in paramFile: anything after a ! on a line is + ignored + + PRE: paramFile is open and file pointer points to start of file + + NOTE: for compatibility with the prior format, the following structure is + also acceptable: name value changeable min max sigma When encountered, + the extra columns are ignored and a warning is generated. + + */ +void readModelParams(ModelParams *params, FILE *paramFile); + +// Return numParameters, the actual number of parameters that have been +// initialized with initializeOneSpatialParam +int getNumModelParams(ModelParams *params); + +// Return number of parameters that have been read in from file +int getNumModelParamsRead(ModelParams *params); + +// Find parameter with given name in the parameters vector +// If found, return index in vector, otherwise return -1 +int locateParam(ModelParams *params, char *name); + +// Return 1 if parameter i has had its value set, 0 otherwise +int valueSet(ModelParams *params, int i); + +/* Check to see whether value is within allowable range of parameter i + Return 1 if okay, 0 if bad +*/ +int checkParam(ModelParams *params, int i, double value); + +/* Return value of parameter i + PRE: 0 <= i < modelParams->numParameters +*/ +double getParam(ModelParams *params, int i); + +// Clean up: deallocate spatialParams and any other dynamically-allocated +// pointers that need deallocating +void deleteModelParams(ModelParams *params); + +#endif diff --git a/src/common/spatialParams.c b/src/common/spatialParams.c deleted file mode 100644 index 776c36317..000000000 --- a/src/common/spatialParams.c +++ /dev/null @@ -1,776 +0,0 @@ -/* spatialParams: structure and functions to track and manipulate parameters - that may vary spatially - - Author: Bill Sacks - - Creation date: 8/13/04 -*/ - -#include -#include -#include -#include -#include "spatialParams.h" // includes definition of SpatialParams structure, as well as OneSpatialParam structure -#include "util.h" - -// Private/helper functions: not defined in spatialParams.h - -// return 1 if all maxParameters parameters have been initialized, 0 if not -int allParamsInitialized(SpatialParams *spatialParams) { - if (spatialParams->numParameters == spatialParams->maxParameters) { - // all parameters have been initialized - return 1; - } - - return 0; -} - -// Set array[0..length-1] to all be equal to value -void setAll(double *array, int length, double value) { - int i; - - for (i = 0; i < length; i++) { - array[i] = value; - } -} - -// Checks to make sure that all parameters with isRequired=true have been read -// If not, kills program -// Writes out names of all parameters that weren't read (even if not required, -// as a warning message) -void checkAllRead(SpatialParams *spatialParams) { - int i; - int okay; - OneSpatialParam *param; - - okay = 1; // so far so good - for (i = 0; i < spatialParams->numParameters; i++) { - param = &(spatialParams->parameters[i]); - if (param->value == NULL) { - if (param->isRequired) { // should have been read but wasn't! - okay = 0; - printf("ERROR: Didn't read required parameter %s\n", param->name); - } else { - printf("WARNING: Didn't read parameter %s (not flagged as required, so " - "continuing)\n", - param->name); - } - } - } - - if (!okay) { - exit(1); - } -} - -/* If isSpatial = 0, return array[0] (ignore loc) - Otherwise, return array[loc] -*/ -double getPossiblySpatial(double *array, int loc, int isSpatial) { - if (isSpatial) { - return array[loc]; - } - - return array[0]; -} - -/* If isSpatial = 0, set array[0] = value (ignore loc) - Otherwise, if loc = -1, set array[0..length-1] = value - Otherwise, set array[loc] = value -*/ -void setPossiblySpatial(double *array, int loc, double value, int isSpatial, - int length) { - if (isSpatial) { - if (loc == -1) { - setAll(array, length, value); - } else { - array[loc] = value; - } - } else { - // non-spatial - array[0] = value; - } -} - -/* Read array of spatial values from file, for a single parameter - Values are returned in the array spatialValues, which must be at least of - length numLocs spatialParamFile is the (already-open) file of spatial - parameters - - the file pointer should be pointing to the start of the line we want to - read pName gives the name of the parameter; we'll abort if the next parameter - does not have the name given by pName numLocs is the number of values to read - */ -void readSpatialValues(double *spatialValues, FILE *spatialParamFile, - char pName[PARAM_MAXNAME], int numLocs) { - char pNameSpatial[PARAM_MAXNAME]; - int status; - int i; - - // read name from spatialParamFile, make sure it's what we expect: - status = fscanf(spatialParamFile, "%s", pNameSpatial); - if (status == EOF || status == 0) { // error reading - printf("Error trying to read %s from spatialParamFile\n", pName); - exit(1); - } else if (strcmpIgnoreCase(pName, pNameSpatial) != 0) { // compare, ignore - // case - printf("Error: read %s from spatialParamFile, expected %s\n", pNameSpatial, - pName); - exit(1); - } - - // now read values and assign param->value, guess, best, and knob: - for (i = 0; i < numLocs; i++) { - status = fscanf(spatialParamFile, "%lf", &(spatialValues[i])); - if (status == EOF || status == 0) { // error reading - printf("Error: did not find enough values in spatialParamFile for %s\n", - pName); - exit(1); - } - } -} - -/*************************************************/ - -// Public functions: defined in spatialParams.h - -// allocate space for a new spatialParams structure, return a pointer to it -SpatialParams *newSpatialParams(int maxParameters, int numLocs) { - SpatialParams *spatialParams; - - spatialParams = (SpatialParams *)malloc(sizeof(SpatialParams)); - spatialParams->maxParameters = maxParameters; - spatialParams->numParameters = 0; // for now, anyway - spatialParams->numParamsRead = 0; // for now, anyway - spatialParams->numLocs = numLocs; - spatialParams->numChangeableParams = 0; // for now, anyway - - // allocate space for vectors within this structure: - spatialParams->parameters = - (OneSpatialParam *)malloc(maxParameters * sizeof(OneSpatialParam)); - spatialParams->readIndices = - (int *)malloc(maxParameters * sizeof(int)); // the biggest it could have - // to be - spatialParams->changeableParamIndices = - (int *)malloc(maxParameters * sizeof(int)); // the biggest it could have - // to be - - return spatialParams; -} - -/* Initialize next spatial parameter: - Set name of parameter equal to "name", - and set parameter's externalLoc pointer to point to the location given by - "externalLoc" parameter Also set parameter's isRequired value: if true, then - we'll terminate execution if this parameter is not read from file Note: - acceptable to have externalLoc = NULL, but then it won't be assigned by - loadSpatialParams -*/ -void initializeOneSpatialParam(SpatialParams *spatialParams, char *name, - double *externalLoc, int isRequired) { - int paramIndex; // index of next uninitialized parameter - OneSpatialParam *param; // a pointer to the next uninitialized parameter - - if (allParamsInitialized(spatialParams)) { // all parameters have been - // initialized! - printf("Error trying to initialize %s: have already initialized all %d " - "parameters\n", - name, spatialParams->maxParameters); - printf( - "Check value of maxParameters passed into newSpatialParams function\n"); - exit(1); - } - - // otherwise, get the index of the next uninitialized parameter - paramIndex = spatialParams->numParameters; - // and set param to point to it for easier access - param = &(spatialParams->parameters[paramIndex]); - - strcpy(param->name, name); - param->externalLoc = externalLoc; - param->isRequired = isRequired; - - // set value pointer to be NULL (this will mark unread parameter) - param->value = NULL; - - spatialParams->numParameters++; -} - -/* Read all parameters from file, setting all values in spatialParams structure - appropriately - - Structure of paramFile: - name value changeable min max sigma - where value is a number if parameter is non-spatial, or * if parameter is - spatial; changeable is 0 or 1 (boolean) If parameter is spatial, look in - spatialParamFile for value, which has the following structure: name value - [value value...] where the number of values is equal to the number of - spatial locations - - The order of the parameters in paramFile does not matter, - but parameters must be specified in the same order in paramFile and - spatialParamFile - - ! is a comment character in paramFile: anything after a ! on a line is - ignored note, though, that comments are not currently allowed in - spatialParamFile - - PRE: paramFile is open and file pointer points to start of file - spatialParamFile is open and file pointer points to 2nd line (after the - numLocs line) - */ -void readSpatialParams(SpatialParams *spatialParams, FILE *paramFile, - FILE *spatialParamFile) { - const char *SEPARATORS = " \t\n\r"; // characters that can separate values in - // parameter files - const char *COMMENT_CHARS = "!"; // comment characters (ignore everything - // after this on a line) - - char line[256]; - char pName[PARAM_MAXNAME]; // parameter name - int paramIndex; - OneSpatialParam *param; // a pointer to a single parameter, for easier access - char strValue[32]; // before we know whether value is a number or "*" - double value, min, max, sigma; - double *spatialValues; // temporary storage for values read from spatial - // param file - int changeable, spatiallyVarying; // is the parameter changeable? is it - // spatially-varying? - char *errc; - int isComment; - int numLocs, i; - - numLocs = spatialParams->numLocs; // for those parameters that are - // spatially-varying - spatialValues = (double *)malloc(numLocs * sizeof(double)); - - while (fgets(line, sizeof(line), paramFile) != NULL) { // while not EOF or - // error - // remove trailing comments: - isComment = stripComment(line, COMMENT_CHARS); - - if (!isComment) { // if this isn't just a comment line or blank line - // tokenize line: - strcpy(pName, strtok(line, SEPARATORS)); // copy first token into pName - strcpy(strValue, strtok(NULL, SEPARATORS)); // copy next token into - // strValue; wait until later - // to figure out if it's "*" - // or a number - changeable = strtol(strtok(NULL, SEPARATORS), &errc, 0); - min = strtod(strtok(NULL, SEPARATORS), &errc); - max = strtod(strtok(NULL, SEPARATORS), &errc); - sigma = strtod(strtok(NULL, SEPARATORS), &errc); - - // now we need to see if we read in an actual value, or a "*" (if the - // latter, it's a spatially-varying parameter) note that we read this - // before doing the error checking on paramIndex: - // that way, even if we ignore this parameter, we've still read (and - // ignored) the corresponding line in the spatial param file, if it's a - // spatially-varying parameter - if (strcmp(strValue, "*") == 0) { // we read "*": this is a - // spatially-varying parameter - readSpatialValues(spatialValues, spatialParamFile, pName, numLocs); - spatiallyVarying = 1; - } else { - spatiallyVarying = 0; - } - - paramIndex = locateParam(spatialParams, pName); - - if (paramIndex == -1) { // not found - printf("WARNING: Ignoring parameter %s: this parameter wasn't " - "initialized in the code\n", - pName); - } else if (valueSet(spatialParams, paramIndex)) { - printf("Error reading parameter file: read %s, but this parameter has " - "already been set\n", - pName); - exit(1); - } else { // otherwise, we're good to go - - // set param to point to the appropriate parameter, for easier access - param = &(spatialParams->parameters[paramIndex]); - - // mark this as the next parameter read from file: - spatialParams->readIndices[spatialParams->numParamsRead] = paramIndex; - spatialParams->numParamsRead++; - - // fill the new spatialParam structure with changeable, min, max, and - // sigma - param->isChangeable = changeable; - param->min = min; - param->max = max; - param->sigma = sigma; - - if (changeable) { // we need to update info on changeable params in - // spatialParams - spatialParams - ->changeableParamIndices[spatialParams->numChangeableParams] = - paramIndex; - // before we change it, spatialParams->numChangeableParams is the - // index of the first free spot in the changeableParamIndices vector - spatialParams->numChangeableParams++; - } - - if (spatiallyVarying) { - param->numLocs = numLocs; - - // allocate space: - param->value = (double *)malloc(numLocs * sizeof(double)); - param->guess = (double *)malloc(numLocs * sizeof(double)); - param->best = (double *)malloc(numLocs * sizeof(double)); - param->knob = (double *)malloc(numLocs * sizeof(double)); - - // now assign param->value, guess, best, and knob (these were read - // above, in call to readSpatialValues): - for (i = 0; i < numLocs; i++) { - // assign value, guess and best to all be the guess value initially: - param->value[i] = param->guess[i] = param->best[i] = - spatialValues[i]; - param->knob[i] = 0; // assign knob to be 0 initially - } // for i - } // if spatially-varying - else { // non-spatially-varying parameter - param->numLocs = 0; // signifies non-spatially-varying - - // allocate space for a single value in each of param->value, guess, - // best and knob: - param->value = (double *)malloc(sizeof(double)); - param->guess = (double *)malloc(sizeof(double)); - param->best = (double *)malloc(sizeof(double)); - param->knob = (double *)malloc(sizeof(double)); - - // assign value, guess and best to all be the guess value initially: - value = strtod(strValue, &errc); // convert string value to double - param->value[0] = param->guess[0] = param->best[0] = value; - param->knob[0] = 0; // assign knob to be 0 initially - } // else non-spatially-varying - } // else (no errors in reading this line from parameter file) - } // if !isComment - } // while not EOF or error - - // check for error in reading: - if (ferror(paramFile)) { - printf("Error reading file in readSpatialParams\n"); - printf("ferror = %d\n", ferror(paramFile)); - exit(1); - } - - checkAllRead(spatialParams); // terminate program if some required parameters - // weren't read - - free(spatialValues); -} // readSpatialParams - -// Return numParameters, the actual number of parameters that have been -// initialized with initializeOneSpatialParam -int getNumParameters(SpatialParams *spatialParams) { - return spatialParams->numParameters; -} - -// Return number of parameters that have been read in from file -int getNumParamsRead(SpatialParams *spatialParams) { - return spatialParams->numParamsRead; -} - -// Find parameter with given name in the parameters vector -// If found, return index in vector, otherwise return -1 -int locateParam(SpatialParams *spatialParams, char *name) { - int i; - int found; - - i = 0; - found = 0; - while (i < spatialParams->numParameters && !found) { - if (strcmpIgnoreCase(name, spatialParams->parameters[i].name) == 0) { - found = 1; - } - i++; - } - - if (found) { - return (i - 1); - } - - return -1; -} - -// Return 1 if parameter i has had its value set, 0 otherwise -int valueSet(SpatialParams *spatialParams, int i) { - return (spatialParams->parameters[i].value != NULL); -} - -/* Return 1 if parameter i varies spatially, 0 if not - PRE: 0 <= i < spatialParams->numParameters -*/ -int isSpatial(SpatialParams *spatialParams, int i) { - return (spatialParams->parameters[i].numLocs > 0); -} - -/* Return 1 if parameter i is changeable, 0 if not - PRE: 0 <= i < spatialParams->numParameters -*/ -int isChangeable(SpatialParams *spatialParams, int i) { - return spatialParams->parameters[i].isChangeable; -} - -/* Check to see whether value is within allowable range of parameter i - Return 1 if okay, 0 if bad -*/ -int checkSpatialParam(SpatialParams *spatialParams, int i, double value) { - return (value >= getSpatialParamMin(spatialParams, i) && - value <= getSpatialParamMax(spatialParams, i)); -} - -/* Return value of parameter i, location loc (if this parameter is non-spatial, - ignore loc) PRE: 0 <= i < spatialParams->numParameters if this parameter is - spatial, then 0 <= loc < spatialParams->numLocs -*/ -double getSpatialParam(SpatialParams *spatialParams, int i, int loc) { - return getPossiblySpatial(spatialParams->parameters[i].value, loc, - isSpatial(spatialParams, i)); -} - -/* Set parameter i, location loc to have given value - If parameter is non-spatial, ignore loc - If loc = -1, set this parameter at ALL LOCATIONS to have given value - PRE: 0 <= i < spatialParams->numParameters - if this parameter is spatial, then -1 <= loc < spatialParams->numLocs -*/ -void setSpatialParam(SpatialParams *spatialParams, int i, int loc, - double value) { - setPossiblySpatial(spatialParams->parameters[i].value, loc, value, - isSpatial(spatialParams, i), spatialParams->numLocs); -} - -/* Return guess value of parameter i, location loc (if this parameter is - non-spatial, ignore loc) PRE: 0 <= i < spatialParams->numParameters if this - parameter is spatial, then 0 <= loc < spatialParams->numLocs -*/ -double getSpatialParamGuess(SpatialParams *spatialParams, int i, int loc) { - return getPossiblySpatial(spatialParams->parameters[i].guess, loc, - isSpatial(spatialParams, i)); -} - -/* Return min. value of range of parameter i (note: this is constant across - space) PRE: 0 <= i < spatialParams->numParameters -*/ -double getSpatialParamMin(SpatialParams *spatialParams, int i) { - return spatialParams->parameters[i].min; -} - -/* Return max. value of range of parameter i (note: this is constant across - space) PRE: 0 <= i < spatialParams->numParameters -*/ -double getSpatialParamMax(SpatialParams *spatialParams, int i) { - return spatialParams->parameters[i].max; -} - -/* Return sigma value of parameter i (note: this is constant across space) - PRE: 0 <= i < spatialParams->numParameters -*/ -double getSpatialParamSigma(SpatialParams *spatialParams, int i) { - return spatialParams->parameters[i].sigma; -} - -/* Return knob value of parameter i, location loc (if this parameter is - non-spatial, ignore loc) PRE: 0 <= i < spatialParams->numParameters if this - parameter is spatial, then 0 <= loc < spatialParams->numLocs -*/ -double getSpatialParamKnob(SpatialParams *spatialParams, int i, int loc) { - return getPossiblySpatial(spatialParams->parameters[i].knob, loc, - isSpatial(spatialParams, i)); -} - -/* Set "knob" of parameter i, location loc to have given value - If parameter is non-spatial, ignore loc - If loc = -1, set knob value of this parameter at ALL LOCATIONS to have given - value PRE: 0 <= i < spatialParams->numParameters if this parameter is - spatial, then -1 <= loc < spatialParams->numLocs -*/ -void setSpatialParamKnob(SpatialParams *spatialParams, int i, int loc, - double value) { - setPossiblySpatial(spatialParams->parameters[i].knob, loc, value, - isSpatial(spatialParams, i), spatialParams->numLocs); -} - -/* Return best value of parameter i, location loc (if this parameter is - non-spatial, ignore loc) PRE: 0 <= i < spatialParams->numParameters if this - parameter is spatial, then 0 <= loc < spatialParams->numLocs -*/ -double getSpatialParamBest(SpatialParams *spatialParams, int i, int loc) { - return getPossiblySpatial(spatialParams->parameters[i].best, loc, - isSpatial(spatialParams, i)); -} - -/* Set best value of parameter i, location loc to have given value - If parameter is non-spatial, ignore loc - If loc = -1, set best value of this parameter at ALL LOCATIONS to have given - value PRE: 0 <= i < spatialParams->numParameters if this parameter is - spatial, then -1 <= loc < spatialParams->numLocs -*/ -void setSpatialParamBest(SpatialParams *spatialParams, int i, int loc, - double value) { - setPossiblySpatial(spatialParams->parameters[i].best, loc, value, - isSpatial(spatialParams, i), spatialParams->numLocs); -} - -/* Set best value of all parameters at given location to be equal to their - current value (Note: ignores parameters that were never read in) If loc = -1, - set best value of all parameters at ALL locations to be equal to their - current value -*/ -void setAllSpatialParamBests(SpatialParams *spatialParams, int loc) { - int numParams; - int firstLoc, lastLoc; - int i, currLoc; - double value; - - numParams = spatialParams->numParameters; - if (loc == -1) { // set best at all locations - firstLoc = 0; - lastLoc = spatialParams->numLocs - 1; - } else { - // set best at a single location - firstLoc = lastLoc = loc; - } - - for (i = 0; i < numParams; i++) { - if (valueSet(spatialParams, i)) { // if valueSet is false, ignore this - // parameter (it was never read in) - if (isSpatial(spatialParams, i)) { - for (currLoc = firstLoc; currLoc <= lastLoc; currLoc++) { - value = getSpatialParam(spatialParams, i, currLoc); - setSpatialParamBest(spatialParams, i, currLoc, value); - } - } else { // non-spatial - value = getSpatialParam(spatialParams, i, 0); - setSpatialParamBest(spatialParams, i, 0, value); - } - } - } -} - -/* Return the index (into spatialParams->parameters) of a randomly-chosen - spatial parameter that is changeable - PRE: random number generator has been seeded -*/ -int randomChangeableSpatialParam(SpatialParams *spatialParams) { - int rnd; - - rnd = (int)floor( - spatialParams->numChangeableParams * - (rand() * 1.0 / (RAND_MAX + 1.0))); // 0 <= rnd < - // spatialParams->numChangeableParams - - return spatialParams->changeableParamIndices[rnd]; -} - -/* Load spatial parameters into memory locations pointed to by paramPtrs - For spatially-varying parameters, load value at location given by loc - (Purpose: load all parameters into a localized structure for a model run, to - make the model run more efficient) - - PRE: spatialParams->parameters is loaded with param. values (particularly, at - given location) externalLoc pointers have been set for each parameter 0 <= loc - < spatialParams->numLocs -*/ -void loadSpatialParams(SpatialParams *spatialParams, int loc) { - int i, numParams; - double value; - - numParams = spatialParams->numParameters; - for (i = 0; i < numParams; i++) { - // externalLoc will equal NULL if it was assigned to NULL in - // readOneSpatialParam - if (valueSet(spatialParams, i) && - spatialParams->parameters[i].externalLoc != NULL) { - value = getSpatialParam(spatialParams, i, loc); - // copy value into location pointed to by this externalLoc - *(spatialParams->parameters[i].externalLoc) = value; - } - } -} - -/* If randomReset = 0, set values of all parameters to be equal to guess values - If randomReset non-zero, set parameter values to be somewhere (chosen uniform - randomly) between min and max - - Note that non-changeable parameters will still be set to their guess - values, though (Note: ignores parameters that were never read in) (If - randomReset non-zero, random number generator must have been seeded) - Set best values of all parameters to be equal to current values Set knobs of - all parameters to be equal to knob argument -*/ -void resetSpatialParams(SpatialParams *spatialParams, double knob, - int randomReset) { - int numParams, numLocs; - int i, loc; - double value; - - numParams = spatialParams->numParameters; - numLocs = spatialParams->numLocs; - - for (i = 0; i < numParams; i++) { - if (valueSet(spatialParams, i)) { // if valueSet is false, ignore this - // parameter (it was never read in) - if (isSpatial(spatialParams, i)) { - for (loc = 0; loc < numLocs; loc++) { - if ((randomReset == 0) || !isChangeable(spatialParams, i)) { - // non-random: use guess value - value = getSpatialParamGuess(spatialParams, i, loc); - } else { - // random: set parameter value to be somewhere between min and max - value = getSpatialParamMin(spatialParams, i) + - (getSpatialParamMax(spatialParams, i) - - getSpatialParamMin(spatialParams, i)) * - ((float)rand() / (float)RAND_MAX); - } - setSpatialParam(spatialParams, i, loc, value); - setSpatialParamBest(spatialParams, i, loc, value); - setSpatialParamKnob(spatialParams, i, loc, knob); - } - } else { // non-spatial - if ((randomReset == 0) || !isChangeable(spatialParams, i)) { - // non-random: use guess value - value = getSpatialParamGuess(spatialParams, i, 0); - } else { - // random: set parameter value to be somewhere between min and max - value = getSpatialParamMin(spatialParams, i) + - (getSpatialParamMax(spatialParams, i) - - getSpatialParamMin(spatialParams, i)) * - ((float)rand() / (float)RAND_MAX); - } - setSpatialParam(spatialParams, i, 0, value); - setSpatialParamBest(spatialParams, i, 0, value); - setSpatialParamKnob(spatialParams, i, 0, knob); - } // else non-spatial - } // if valueSet - } // for i -} // resetSpatialParams - -/* Write best parameter values, and other parameter info, to files - (only write parameters that were read in, and in the same order that they - were read in) Note that other parameter info won't have changed since the - read - we just copy that over Write all parameter info to file with name - *paramFile, and spatially-varying parameter values to file with name - *spatialParamFile See readSpatialParams for file formats -*/ -void writeBestSpatialParams(SpatialParams *spatialParams, char *paramFile, - char *spatialParamFile) { - FILE *paramF, *spatialParamF; - int numParamsRead, numLocs; - int i, j, index; - - paramF = openFile(paramFile, "w"); - spatialParamF = openFile(spatialParamFile, "w"); - - numParamsRead = getNumParamsRead(spatialParams); - numLocs = spatialParams->numLocs; - - fprintf(spatialParamF, "%d\n", numLocs); - - for (i = 0; i < numParamsRead; i++) { // loop through parameters - index = spatialParams->readIndices[i]; // allows us to just consider the - // parameters that were read in, in - // the proper order - fprintf(paramF, "%s\t", spatialParams->parameters[index].name); - if (isSpatial(spatialParams, index)) { - fprintf(paramF, "*\t"); // write "*" in place of value, then write values - // to spatialParamFile - fprintf(spatialParamF, "%s", spatialParams->parameters[index].name); - for (j = 0; j < numLocs; j++) { - // loop through locations - fprintf(spatialParamF, "\t%f", - getSpatialParamBest(spatialParams, index, j)); - } - fprintf(spatialParamF, "\n"); - } else { - // non-spatial: write value - fprintf(paramF, "%f\t", getSpatialParamBest(spatialParams, index, 0)); - } - - // now write other info to file (other info is always non-spatial): - fprintf(paramF, "%d\t%f\t%f\t%f\n", isChangeable(spatialParams, index), - getSpatialParamMin(spatialParams, index), - getSpatialParamMax(spatialParams, index), - getSpatialParamSigma(spatialParams, index)); - } - - fclose(paramF); - fclose(spatialParamF); -} - -/* Write name, guess, min & max, value, best and knob of each changeable param - to file (only write parameters that were read in, and in the same order that - they were read in) If loc >= 0, write info for that location; if loc = -1, - write value at location 0, and append a * to spatial parameter values PRE: - outF is open for writing -*/ -void writeChangeableParamInfo(SpatialParams *spatialParams, int loc, - FILE *outF) { - int i, index, numParamsRead, count; - int theLoc; // location that we're printing info for - char appendage[2]; - - if (loc == -1) { - theLoc = 0; - } else { - theLoc = loc; - } - - numParamsRead = getNumParamsRead(spatialParams); - count = 0; // count of changeable parameters - fprintf(outF, "Changeable Parameters:\n"); - for (i = 0; i < numParamsRead; i++) { - // allows us to just consider the parameters that were read in, in - // the proper order - index = spatialParams->readIndices[i]; - if (isChangeable(spatialParams, index)) { - if (loc == -1 && isSpatial(spatialParams, index)) { - // we'll append a * to spatially-varying values - strcpy(appendage, "*"); - } else { - strcpy(appendage, ""); - } - - fprintf(outF, "[%d] %s:\t", count, spatialParams->parameters[index].name); - fprintf(outF, "Value = %f%s\t", - getSpatialParam(spatialParams, index, theLoc), appendage); - fprintf(outF, "Best = %f%s\t", - getSpatialParamBest(spatialParams, index, theLoc), appendage); - fprintf(outF, "Knob = %f%s\t", - getSpatialParamKnob(spatialParams, index, theLoc), appendage); - fprintf(outF, "Prior: %f%s ", - getSpatialParamGuess(spatialParams, index, theLoc), appendage); - fprintf(outF, " [%f, %f]\n", getSpatialParamMin(spatialParams, index), - getSpatialParamMax(spatialParams, index)); - - count++; - } - } -} - -// Clean up: deallocate spatialParams and any other dynamically-allocated -// pointers that need deallocating -void deleteSpatialParams(SpatialParams *spatialParams) { - int numParameters, i; - OneSpatialParam *param; - - // first loop through individual parameters, de-allocating space they use: - numParameters = spatialParams->numParameters; - for (i = 0; i < numParameters; i++) { - param = &(spatialParams->parameters[i]); - if (param->value != NULL) { // this parameter has been used - so space has - // been allocated for its various values - free(param->value); - free(param->guess); - free(param->best); - free(param->knob); - } - } - - // now free space used by spatialParams structure itself: - free(spatialParams->parameters); - free(spatialParams->readIndices); - free(spatialParams->changeableParamIndices); - free(spatialParams); -} diff --git a/src/common/spatialParams.h b/src/common/spatialParams.h deleted file mode 100644 index 6ac8e8659..000000000 --- a/src/common/spatialParams.h +++ /dev/null @@ -1,264 +0,0 @@ -// header file for spatialParams.c -// includes definition of SpatialParams structure - -#ifndef SPATIAL_PARAMS_H -#define SPATIAL_PARAMS_H - -#include - -#define PARAM_MAXNAME 64 - -// struct to hold a single (possibly) spatially-varying param -typedef struct OneSpatialParamStruct { - char name[PARAM_MAXNAME]; // name of parameter - int isRequired; // 0 or 1; 1 indicates that we should terminate run if this - // parameter isn't specified in input file - int numLocs; // 0 if this parameter is non-spatial - int isChangeable; // 0 or 1 - double *value; /* if non-spatial, points to a single double - if spatial, points to a vector of doubles - before it's set, value = NULL */ - double *externalLoc; // a pointer to a (non-spatial) copy of this parameter - // value elsewhere (e.g. local to the model) - - // items that we care about when doing parameter optimization: - // parameter info read from file: - double *guess; // possibly spatial; if non-changeable, ptr will be void - double min; // left-hand side of allowable range - double max; // right-hand side of allowable range - double sigma; // parameter uncertainty - - // variables related to optimization: - double *best; // best value found so far; possibly patial; if non-changeable, - // ptr will be void - double *knob; /* a "knob" that can be twisted to adjust this parameter's - optimization (can be used however the individual optimization scheme - wants) possibly spatial; if non-changeable, ptr will be void */ - -} OneSpatialParam; - -typedef struct SpatialParamsStruct { - OneSpatialParam *parameters; // vector of parameters, dynamically-allocated - // (stored in order in which they were - // initialized) - int maxParameters; // maximum number of parameters - int numParameters; // actual number of parameters (tracks # that have been - // initialized with initializeOneSpatialParam) - int numParamsRead; // tracks number of parameters that have been read in from - // file - int *readIndices; // indices (in parameters vector) of the parameters read - // from file, in the order in which they were read (this - // will allow us to output parameters in the same order - // later) - int numChangeableParams; // number of changeable params (in an optimization) - int *changeableParamIndices; // indices (in parameters vector) of the - // changeable params (length of this vector is - // numChangeableParams) - int numLocs; // number of spatial locations -} SpatialParams; - -// allocate space for a new spatialParams structure, return a pointer to it -SpatialParams *newSpatialParams(int maxParameters, int numLocs); - -/* Initialize next spatial parameter: - Set name of parameter equal to "name", - and set parameter's externalLoc pointer to point to the location given by - "externalLoc" parameter - Also set parameter's isRequired value: if true, then we'll terminate - execution if this parameter is not read from file - Note: acceptable to have externalLoc = NULL, but then it won't be assigned by - loadSpatialParams -*/ -void initializeOneSpatialParam(SpatialParams *spatialParams, char *name, - double *externalLoc, int isRequired); - -/* Read all parameters from file, setting all values in spatialParams structure - appropriately - - Structure of paramFile: - name value changeable min max sigma - where - value is a number if parameter is non-spatial, or * if parameter is spatial; - changeable is 0 or 1 (boolean) - - If parameter is spatial, look in spatialParamFile for value, which has the - following structure: name value [value value...] where the number of values - is equal to the number of spatial locations - - The order of the parameters in paramFile does not matter, - but parameters must be specified in the same order in paramFile and - spatialParamFile - - ! is a comment character in paramFile: anything after a ! on a line is - ignored note, though, that comments are not currently allowed in - spatialParamFile - - PRE: paramFile is open and file pointer points to start of file - spatialParamFile is open and file pointer points to 2nd line (after the - numLocs line) - */ -void readSpatialParams(SpatialParams *spatialParams, FILE *paramFile, - FILE *spatialParamFile); - -// Return numParameters, the actual number of parameters that have been -// initialized with initializeOneSpatialParam -int getNumParameters(SpatialParams *spatialParams); - -// Return number of parameters that have been read in from file -int getNumParamsRead(SpatialParams *spatialParams); - -// Find parameter with given name in the parameters vector -// If found, return index in vector, otherwise return -1 -int locateParam(SpatialParams *spatialParams, char *name); - -// Return 1 if parameter i has had its value set, 0 otherwise -int valueSet(SpatialParams *spatialParams, int i); - -/* Return 1 if parameter i varies spatially, 0 if not - PRE: 0 <= i < spatialParams->numParameters -*/ -int isSpatial(SpatialParams *spatialParams, int i); - -/* Return 1 if parameter i is changeable, 0 if not - PRE: 0 <= i < spatialParams->numParameters -*/ -int isChangeable(SpatialParams *spatialParams, int i); - -/* Check to see whether value is within allowable range of parameter i - Return 1 if okay, 0 if bad -*/ -int checkSpatialParam(SpatialParams *spatialParams, int i, double value); - -/* Return value of parameter i, location loc (if this parameter is - non-spatial, ignore loc) - PRE: 0 <= i < spatialParams->numParameters - if this parameter is spatial, then 0 <= loc < spatialParams->numLocs -*/ -double getSpatialParam(SpatialParams *spatialParams, int i, int loc); - -/* Set parameter i, location loc to have given value - If parameter is non-spatial, ignore loc - If loc = -1, set this parameter at ALL LOCATIONS to have given value - PRE: 0 <= i < spatialParams->numParameters - if this parameter is spatial, then -1 <= loc < spatialParams->numLocs -*/ -void setSpatialParam(SpatialParams *spatialParams, int i, int loc, - double value); - -/* Return guess value of parameter i, location loc (if this parameter is - non-spatial, ignore loc) - PRE: 0 <= i < spatialParams->numParameters - if this parameter is spatial, then 0 <= loc < spatialParams->numLocs -*/ -double getSpatialParamGuess(SpatialParams *spatialParams, int i, int loc); - -/* Return min. value of range of parameter i (note: this is constant across - space) PRE: 0 <= i < spatialParams->numParameters -*/ -double getSpatialParamMin(SpatialParams *spatialParams, int i); - -/* Return max. value of range of parameter i (note: this is constant across - space) PRE: 0 <= i < spatialParams->numParameters -*/ -double getSpatialParamMax(SpatialParams *spatialParams, int i); - -/* Return sigma value of parameter i (note: this is constant across space) - PRE: 0 <= i < spatialParams->numParameters -*/ -double getSpatialParamSigma(SpatialParams *spatialParams, int i); - -/* Return knob value of parameter i, location loc (if this parameter is - non-spatial, ignore loc) - PRE: 0 <= i < spatialParams->numParameters - if this parameter is spatial, then 0 <= loc < spatialParams->numLocs -*/ -double getSpatialParamKnob(SpatialParams *spatialParams, int i, int loc); - -/* Set "knob" of parameter i, location loc to have given value - If parameter is non-spatial, ignore loc - If loc = -1, set knob value of this parameter at ALL LOCATIONS to have given - value PRE: 0 <= i < spatialParams->numParameters if this parameter is - spatial, then -1 <= loc < spatialParams->numLocs -*/ -void setSpatialParamKnob(SpatialParams *spatialParams, int i, int loc, - double value); - -/* Return best value of parameter i, location loc (if this parameter is - non-spatial, ignore loc) - PRE: 0 <= i < spatialParams->numParameters - if this parameter is spatial, then 0 <= loc < spatialParams->numLocs -*/ -double getSpatialParamBest(SpatialParams *spatialParams, int i, int loc); - -/* Set best value of parameter i, location loc to have given value - If parameter is non-spatial, ignore loc - If loc = -1, set best value of this parameter at ALL LOCATIONS to have given - value - PRE: 0 <= i < spatialParams->numParameters if this parameter is spatial, - then -1 <= loc < spatialParams->numLocs -*/ -void setSpatialParamBest(SpatialParams *spatialParams, int i, int loc, - double value); - -/* Set best value of all parameters at given location to be equal to their - current value (Note: ignores parameters that were never read in) - If loc = -1, set best value of all parameters at ALL locations to be equal - to their current value -*/ -void setAllSpatialParamBests(SpatialParams *spatialParams, int loc); - -/* Return the index (into spatialParams->parameters) of a randomly-chosen - spatial parameter that is changeable - PRE: random number generator has been seeded -*/ -int randomChangeableSpatialParam(SpatialParams *spatialParams); - -/* Load spatial parameters into memory locations pointed to by paramPtrs - For spatially-varying parameters, load value at location given by loc - (Purpose: load all parameters into a localized structure for a model run, to - make the model run more efficient) - - PRE: spatialParams->parameters is loaded with param. values (particularly, at - given location) externalLoc pointers have been set for each parameter 0 <= loc - < spatialParams->numLocs -*/ -void loadSpatialParams(SpatialParams *spatialParams, int loc); - -/* If randomReset = 0, set values of all parameters to be equal to guess values - If randomReset non-zero, set parameter values to be somewhere (chosen - uniform randomly) between min and max - - Note that non-changeable parameters will still be set to their guess - values, though - (Note: ignores parameters that were never read in) - (If randomReset non-zero, random number generator must have been seeded) - Set best values of all parameters to be equal to current values - Set knobs of all parameters to be equal to knob argument -*/ -void resetSpatialParams(SpatialParams *spatialParams, double knob, - int randomReset); - -/* Write best parameter values, and other parameter info, to files - (only write parameters that were read in, and in the same order that they - were read in) - Note that other parameter info won't have changed since the read - we just - copy that over - Write all parameter info to file with name *paramFile, and spatially-varying - parameter values to file with name *spatialParamFile - See readOneSpatialParam for file formats -*/ -void writeBestSpatialParams(SpatialParams *spatialParams, char *paramFile, - char *spatialParamFile); - -/* Write name, guess, min & max, value, best and knob of each changeable param - to file If loc >= 0, write info for that location; if loc = -1, write value - at location 0, and append a * to spatial parameter values - PRE: outF is open for writing -*/ -void writeChangeableParamInfo(SpatialParams *spatialParams, int loc, - FILE *outF); - -// Clean up: deallocate spatialParams and any other dynamically-allocated -// pointers that need deallocating -void deleteSpatialParams(SpatialParams *spatialParams); - -#endif diff --git a/src/sipnet/events.c b/src/sipnet/events.c index 813d6ee49..03bf82f84 100644 --- a/src/sipnet/events.c +++ b/src/sipnet/events.c @@ -4,13 +4,12 @@ * @brief Handles reading, parsing, and storing agronomic events for SIPNET simulations. * * The `events.in` file specifies agronomic events with the following columns: - * | col | parameter | description | units | - * |-----|-------------|--------------------------------------------|------------------| - * | 1 | loc | spatial location index | | - * | 2 | year | year of start of this timestep | | - * | 3 | day | day of start of this timestep | Day of year | - * | 4 | event_type | type of event | (e.g., "irrig") | - * | 5...n| event_param| parameters specific to the event type | (varies) | + * | col | parameter | description | units | + * |------|-------------|--------------------------------------------|------------------| + * | 1 | year | year of start of this timestep | | + * | 2 | day | day of start of this timestep | Day of year | + * | 3 | event_type | type of event | (e.g., "irrig") | + * | 4...n| event_param | parameters specific to the event type | (varies) | * * * Event types include: @@ -33,10 +32,9 @@ void printEvent(EventNode *event); -EventNode *createEventNode(int loc, int year, int day, int eventType, +EventNode *createEventNode(int year, int day, int eventType, const char *eventParamsStr) { EventNode *event = (EventNode *)malloc(sizeof(EventNode)); - event->loc = loc; event->year = year; event->day = day; event->type = eventType; @@ -45,18 +43,18 @@ EventNode *createEventNode(int loc, int year, int day, int eventType, case HARVEST: { double fracRA, fracRB, fracTA, fracTB; HarvestParams *params = (HarvestParams *)malloc(sizeof(HarvestParams)); - int numRead = sscanf(eventParamsStr, "%lf %lf %lf %lf", &fracRA, &fracRB, - &fracTA, &fracTB); + int numRead = + sscanf(eventParamsStr, // NOLINT + "%lf %lf %lf %lf", &fracRA, &fracRB, &fracTA, &fracTB); if (numRead != NUM_HARVEST_PARAMS) { - printf("Error parsing Harvest params for loc %d year %d day %d\n", loc, - year, day); + printf("Error parsing Harvest params for year %d day %d\n", year, day); exit(EXIT_CODE_INPUT_FILE_ERROR); } // Validate the params if ((fracRA + fracTA > 1) || (fracRB + fracTB > 1)) { - printf("Invalid harvest event for loc %d year %d day %d; above and " - "below must each add to 1 or less", - loc, year, day); + printf("Invalid harvest event for year %d day %d; above and below must " + "each add to 1 or less", + year, day); exit(EXIT_CODE_BAD_PARAMETER_VALUE); } params->fractionRemovedAbove = fracRA; @@ -70,10 +68,11 @@ EventNode *createEventNode(int loc, int year, int day, int eventType, int method; IrrigationParams *params = (IrrigationParams *)malloc(sizeof(IrrigationParams)); - int numRead = sscanf(eventParamsStr, "%lf %d", &amountAdded, &method); + int numRead = sscanf(eventParamsStr, // NOLINT + "%lf %d", &amountAdded, &method); if (numRead != NUM_IRRIGATION_PARAMS) { - printf("Error parsing Irrigation params for loc %d year %d day %d\n", - loc, year, day); + printf("Error parsing Irrigation params for year %d day %d\n", year, + day); exit(EXIT_CODE_INPUT_FILE_ERROR); } params->amountAdded = amountAdded; @@ -87,10 +86,11 @@ EventNode *createEventNode(int loc, int year, int day, int eventType, // double nh4_no3_frac; FertilizationParams *params = (FertilizationParams *)malloc(sizeof(FertilizationParams)); - int numRead = sscanf(eventParamsStr, "%lf %lf %lf", &orgN, &orgC, &minN); + int numRead = sscanf(eventParamsStr, // NOLINT + "%lf %lf %lf", &orgN, &orgC, &minN); if (numRead != NUM_FERTILIZATION_PARAMS) { - printf("Error parsing Fertilization params for loc %d year %d day %d\n", - loc, year, day); + printf("Error parsing Fertilization params for year %d day %d\n", year, + day); exit(EXIT_CODE_INPUT_FILE_ERROR); } // scanf(eventParamsStr, "%lf %lf %lf %lf", &org_N, &org_C, &min_N, @@ -104,11 +104,11 @@ EventNode *createEventNode(int loc, int year, int day, int eventType, case PLANTING: { double leafC, woodC, fineRootC, coarseRootC; PlantingParams *params = (PlantingParams *)malloc(sizeof(PlantingParams)); - int numRead = sscanf(eventParamsStr, "%lf %lf %lf %lf", &leafC, &woodC, - &fineRootC, &coarseRootC); + int numRead = + sscanf(eventParamsStr, // NOLINT + "%lf %lf %lf %lf", &leafC, &woodC, &fineRootC, &coarseRootC); if (numRead != NUM_PLANTING_PARAMS) { - printf("Error parsing Planting params for loc %d year %d day %d\n", loc, - year, day); + printf("Error parsing Planting params for year %d day %d\n", year, day); exit(EXIT_CODE_INPUT_FILE_ERROR); } params->leafC = leafC; @@ -120,11 +120,10 @@ EventNode *createEventNode(int loc, int year, int day, int eventType, case TILLAGE: { double fracLT, somDM, litterDM; TillageParams *params = (TillageParams *)malloc(sizeof(TillageParams)); - int numRead = - sscanf(eventParamsStr, "%lf %lf %lf", &fracLT, &somDM, &litterDM); + int numRead = sscanf(eventParamsStr, // NOLINT + "%lf %lf %lf", &fracLT, &somDM, &litterDM); if (numRead != NUM_TILLAGE_PARAMS) { - printf("Error parsing Tillage params for loc %d year %d day %d\n", loc, - year, day); + printf("Error parsing Tillage params for year %d day %d\n", year, day); exit(EXIT_CODE_INPUT_FILE_ERROR); } params->fractionLitterTransferred = fracLT; @@ -157,18 +156,16 @@ event_type_t getEventType(const char *eventTypeStr) { return UNKNOWN_EVENT; } -EventNode **readEventData(char *eventFile, int numLocs) { - int loc, year, day, eventType; - int currLoc, currYear, currDay; +EventNode *readEventData(char *eventFile) { + int year, day, eventType; + int currYear, currDay; int EVENT_LINE_SIZE = 1024; int numBytes; char *eventParamsStr; char eventTypeStr[20]; char line[EVENT_LINE_SIZE]; EventNode *curr, *next; - - EventNode **events = - (EventNode **)calloc(sizeof(EventNode *), numLocs * sizeof(EventNode *)); + EventNode *events = NULL; // Check for a non-empty file if (access(eventFile, F_OK) != 0) { @@ -187,8 +184,8 @@ EventNode **readEventData(char *eventFile, int numLocs) { return events; } - int numRead = sscanf(line, "%d %d %d %s %n", &loc, &year, &day, eventTypeStr, - &numBytes); + int numRead = sscanf(line, // NOLINT + "%d %d %s %n", &year, &day, eventTypeStr, &numBytes); if (numRead != NUM_EVENT_CORE_PARAMS) { printf("Error reading event file: bad data on first line\n"); exit(EXIT_CODE_INPUT_FILE_ERROR); @@ -201,21 +198,20 @@ EventNode **readEventData(char *eventFile, int numLocs) { exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM); } - next = createEventNode(loc, year, day, eventType, eventParamsStr); - events[loc] = next; - currLoc = loc; + events = createEventNode(year, day, eventType, eventParamsStr); + next = events; currYear = year; currDay = day; while (fgets(line, EVENT_LINE_SIZE, in) != NULL) { // We have another event curr = next; - numRead = sscanf(line, "%d %d %d %s %n", &loc, &year, &day, eventTypeStr, - &numBytes); + numRead = sscanf(line, "%d %d %s %n", // NOLINT + &year, &day, eventTypeStr, &numBytes); if (numRead != 4) { - printf("Error reading event file: bad data on line after loc %d year %d " - "day %d\n", - currLoc, currYear, currDay); + printf( + "Error reading event file: bad data on line after year %d day %d\n", + currYear, currDay); exit(EXIT_CODE_INPUT_FILE_ERROR); } eventParamsStr = line + numBytes; @@ -226,35 +222,16 @@ EventNode **readEventData(char *eventFile, int numLocs) { exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM); } - // make sure location and time are non-decreasing - if (loc < currLoc) { - printf("Error reading event file: was reading location %d, trying to " - "read location %d\n", - currLoc, loc); - printf("Event records for a given location should be contiguous, and " - "locations should be in ascending order\n"); - exit(EXIT_CODE_INPUT_FILE_ERROR); - } - if ((loc == currLoc) && - ((year < currYear) || ((year == currYear) && (day < currDay)))) { - printf("Error reading event file: for location %d, last event was at " - "(%d, %d) ", - currLoc, currYear, currDay); + if ((year < currYear) || ((year == currYear) && (day < currDay))) { + printf("Error reading event file: last event was at (%d, %d) ", currYear, + currDay); printf("next event is at (%d, %d)\n", year, day); - printf("Event records for a given location should be in time-ascending " - "order\n"); + printf("Event records must be in time-ascending order\n"); exit(EXIT_CODE_INPUT_FILE_ERROR); } - next = createEventNode(loc, year, day, eventType, eventParamsStr); - if (currLoc == loc) { - // Same location, add the new event to this location's list - curr->nextEvent = next; - } else { - // New location, update location and start a new list - currLoc = loc; - events[currLoc] = next; - } + next = createEventNode(year, day, eventType, eventParamsStr); + curr->nextEvent = next; currYear = year; currDay = day; } @@ -267,24 +244,23 @@ void printEvent(EventNode *event) { if (event == NULL) { return; } - int loc = event->loc; int year = event->year; int day = event->day; switch (event->type) { case IRRIGATION: - printf("IRRIGATION at loc %d on %d %d, ", loc, year, day); + printf("IRRIGATION on %d %d, ", year, day); IrrigationParams *const iParams = (IrrigationParams *)event->eventParams; printf("with params: amount added %4.2f\n", iParams->amountAdded); break; case FERTILIZATION: - printf("FERTILIZATION at loc %d on %d %d, ", loc, year, day); + printf("FERTILIZATION on %d %d, ", year, day); FertilizationParams *const fParams = (FertilizationParams *)event->eventParams; printf("with params: org N %4.2f, org C %4.2f, min N %4.2f\n", fParams->orgN, fParams->orgC, fParams->minN); break; case PLANTING: - printf("PLANTING at loc %d on %d %d, ", loc, year, day); + printf("PLANTING on %d %d, ", year, day); PlantingParams *const pParams = (PlantingParams *)event->eventParams; printf("with params: leaf C %4.2f, wood C %4.2f, fine root C %4.2f, " "coarse root C %4.2f\n", @@ -292,7 +268,7 @@ void printEvent(EventNode *event) { pParams->coarseRootC); break; case TILLAGE: - printf("TILLAGE at loc %d on %d %d, ", loc, year, day); + printf("TILLAGE at on %d %d, ", year, day); TillageParams *const tParams = (TillageParams *)event->eventParams; printf("with params: frac litter transferred %4.2f, som decomp modifier " "%4.2f, litter decomp modifier %4.2f\n", @@ -300,7 +276,7 @@ void printEvent(EventNode *event) { tParams->litterDecompModifier); break; case HARVEST: - printf("HARVEST at loc %d on %d %d, ", loc, year, day); + printf("HARVEST on %d %d, ", year, day); HarvestParams *const hParams = (HarvestParams *)event->eventParams; printf("with params: frac removed above %4.2f, frac removed below %4.2f, " "frac transferred above %4.2f, frac transferred below %4.2f\n", @@ -336,7 +312,7 @@ FILE *openEventOutFile(int printHeader) { if (printHeader) { // Use format string analogous to the one in writeEventOut for // better alignment (won't be perfect, but definitely better) - fprintf(eventFile, "%3s %4s %3s %-7s %s", "loc", "year", "day", "type", + fprintf(eventFile, "%4s %3s %-7s %s", "year", "day", "type", "param_name=delta[,param_name=delta,...]\n"); } return eventFile; @@ -347,10 +323,10 @@ void writeEventOut(FILE *out, EventNode *event, int numParams, ...) { int ind = 0; // Spec: - // loc year day event_type [,=,...] + // year day event_type [,=,...] // Standard prefix for all - fprintf(out, "%-3d %4d %3d %-7s ", event->loc, event->year, event->day, + fprintf(out, "%4d %3d %-7s ", event->year, event->day, eventTypeToString(event->type)); // Variable output per event type va_start(args, numParams); diff --git a/src/sipnet/events.h b/src/sipnet/events.h index 90669253e..b2fd9f8de 100644 --- a/src/sipnet/events.h +++ b/src/sipnet/events.h @@ -42,7 +42,6 @@ typedef struct FertilizationParams { double orgN; double orgC; double minN; - // double nh4_no3_frac; for two-pool version } FertilizationParams; #define NUM_PLANTING_PARAMS 4 @@ -60,11 +59,11 @@ typedef struct TillageParams { double litterDecompModifier; } TillageParams; -#define NUM_EVENT_CORE_PARAMS 4 +#define NUM_EVENT_CORE_PARAMS 3 typedef struct EventNode EventNode; struct EventNode { event_type_t type; - int loc, year, day; + int year, day; void *eventParams; EventNode *nextEvent; }; @@ -80,12 +79,10 @@ const char *eventTypeToString(event_type_t type); /*! * Read event data from input filename (canonically events.in) * - * Format: returned data is structured as an array of EventNode pointers, - * indexed by location. Each element of the array is the first event for that - * location (or null if there are no events). It is assumed that the events are - * ordered first by location and then by year and day. + * Format: returned data is structured as an linked list of EventNode pointers. + * It is assumed that the events are ordered by year and day. */ -EventNode **readEventData(char *eventFile, int numLocs); +EventNode *readEventData(char *eventFile); /*! * Open the event output file and optionally write a header row @@ -103,7 +100,7 @@ FILE *openEventOutFile(int printHeader); * * Output format: * - * loc year day event_type \=\[,\=\,...] + * year day event_type \=\[,\=\,...] * * \param out FILE pointer to events.out * \param event Pointer to event node diff --git a/src/sipnet/frontend.c b/src/sipnet/frontend.c index 13ea743ff..4c90a2282 100644 --- a/src/sipnet/frontend.c +++ b/src/sipnet/frontend.c @@ -11,7 +11,7 @@ #include "common/exitCodes.h" #include "common/namelistInput.h" -#include "common/spatialParams.h" +#include "common/modelParams.h" #include "common/util.h" #include "events.h" @@ -27,8 +27,6 @@ #define INPUT_FILE "sipnet.in" #define DO_MAIN_OUTPUT 1 #define DO_SINGLE_OUTPUTS 0 -// Default is run at all locations -#define LOC (-1) void checkRuntype(NamelistInputs *namelist, const char *inputFile) { NamelistInputItem *inputItem = locateNamelistInputItem(namelist, "RUNTYPE"); @@ -61,11 +59,14 @@ int main(int argc, char *argv[]) { char inputFile[INPUT_MAXNAME] = INPUT_FILE; NamelistInputs *namelistInputs; + // for compatibility + int loc; + FILE *out; int option; // reading in optional arguments - SpatialParams *spatialParams; // the parameters used in the model (possibly - // spatially-varying) + ModelParams *modelParams; // the parameters used in the model (possibly + // spatially-varying) OutputItems *outputItems; // structure to hold information for output to // single-variable files (if doSingleOutputs is // true) @@ -76,10 +77,6 @@ int main(int argc, char *argv[]) { // variables? int doSingleOutputs = DO_SINGLE_OUTPUTS; // do we do extra outputting of // single-variable files? - int loc = LOC; // location to run at (set through optional -l argument) - int numLocs; // read in initModel - int *steps; // number of time steps in each location - int printHeader = HEADER; char fileName[FILE_MAXNAME]; @@ -114,7 +111,7 @@ int main(int argc, char *argv[]) { namelistInputs = newNamelistInputs(); addNamelistInputItem(namelistInputs, "RUNTYPE", STRING_TYPE, runtype, RUNTYPE_MAXNAME); addNamelistInputItem(namelistInputs, "FILENAME", STRING_TYPE, fileName, FILE_MAXNAME); - addNamelistInputItem(namelistInputs, "LOCATION", INT_TYPE, &loc,0); + addNamelistInputItem(namelistInputs, "LOCATION", INT_TYPE, &loc, 0); addNamelistInputItem(namelistInputs, "DO_MAIN_OUTPUT", INT_TYPE, &doMainOutput, 0); addNamelistInputItem(namelistInputs, "DO_SINGLE_OUTPUTS", INT_TYPE, &doSingleOutputs, 0); addNamelistInputItem(namelistInputs, "PRINT_HEADER", INT_TYPE, &printHeader,0); @@ -135,10 +132,10 @@ int main(int argc, char *argv[]) { strcat(paramFile, ".param"); strcpy(climFile, fileName); strcat(climFile, ".clim"); - numLocs = initModel(&spatialParams, &steps, paramFile, climFile); + initModel(&modelParams, paramFile, climFile); #if EVENT_HANDLER - initEvents(EVENT_IN_FILE, numLocs, printHeader); + initEvents(EVENT_IN_FILE, printHeader); #endif if (doSingleOutputs) { @@ -157,18 +154,16 @@ int main(int argc, char *argv[]) { out = NULL; } - runModelOutput(out, outputItems, printHeader, spatialParams, loc); + runModelOutput(out, outputItems, printHeader); if (doMainOutput) { fclose(out); } - cleanupModel(numLocs); - deleteSpatialParams(spatialParams); + cleanupModel(); if (outputItems != NULL) { deleteOutputItems(outputItems); } - free(steps); return 0; } diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 444fd70bb..600f4be08 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -212,7 +212,7 @@ struct ClimateVars { ClimateNode *nextClim; }; -#define NUM_CLIM_FILE_COLS 14 +#define NUM_CLIM_FILE_COLS 13 // model parameters which can change from one run to the next // these include initializations of state @@ -537,8 +537,7 @@ typedef struct PhenologyTrackersStruct { /* variables to track each year's // made global so they can be initialized in a separate function // so they only have to be initialized once -static ClimateNode **firstClimates; // a vector of pointers to first climates - // of each point in space +static ClimateNode *firstClimate; // pointers to first climate // more global variables: // these are global to increase efficiency (avoid lots of parameter passing) @@ -558,44 +557,30 @@ static ClimateNode *climate; // current climate static Fluxes fluxes; #if EVENT_HANDLER -static EventNode **events; -static EventNode *locEvent; // current location event list +static EventNode *events; +static EventNode *event; static FILE *eventOutFile; #endif -/* Read climate file into linked lists, - make firstClimates be a vector where each element is a pointer to the head of - a list corresponding to one spatial location - - return an array containing the number of time steps in each location - (dynamically allocated with malloc) - - numLocs = number of spatial locations: there should be one set of entries in - climFile for each location, or a location can be skipped, and we'll use - climate at location 0 for that location, too - - format of climFile: - each line represents one time step, with the following format: - - location year day time intervalLength tair tsoil par precip vpd vpdSoil - vPress wspd soilWetness +/*! + * Read climate file into linked list + * + * Each line of the climate file represents one time step, with the following + * format: + * year day time intervalLength tair tsoil par precip vpd vpdSoil vPress wspd + soilWetness + * + * NOTE: there should be NO blank lines in the file. - NOTE: there should be NO blank lines in file, even between different spatial - locations; all entries for a given location should be contiguous, and - locations should be in ascending order. There may be locations for which - there is no climate data, in which case we use climate from location 0 -*/ -int *readClimData(char *climFile, int numLocs) { + * @param climFile Name of climate file + */ +void readClimData(const char *climFile) { FILE *in; ClimateNode *curr, *next; - int loc, year, day; + int year, day; int lastYear = -1; double time, length; // time in hours, length in days (or fraction of day) double tair, tsoil, par, precip, vpd, vpdSoil, vPress, wspd, soilWetness; - int currLoc; - int count; - int i; - int *steps; // # of time steps in each location #if GDD double thisGdd; // growing degree days of this time step @@ -604,49 +589,49 @@ int *readClimData(char *climFile, int numLocs) { int status; // status of the read - steps = (int *)malloc(numLocs * sizeof(int)); - for (i = 0; i < numLocs; i++) { - steps[i] = 0; // 0 will denote nothing read - } + // for format check + int firstLoc, dummyLoc; + int expectedNumCols = NUM_CLIM_FILE_COLS; + ssize_t numCharsRead; + char *firstLine; + int hasLoc = 0; in = openFile(climFile, "r"); - status = fscanf(in, "%d %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", - &loc, &year, &day, &time, &length, &tair, &tsoil, &par, - &precip, &vpd, &vpdSoil, &vPress, &wspd, &soilWetness); - if (status == EOF) { + // Check format of first line to see if location is still specified (we will + // ignore it if so) + numCharsRead = getline(&firstLine, 0, in); + if (numCharsRead == -1) { // EOF printf("Error: no climate data in %s\n", climFile); exit(EXIT_CODE_INPUT_FILE_ERROR); } + status = sscanf(firstLine, // NOLINT + "%d %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", + &firstLoc, &year, &day, &time, &length, &tair, &tsoil, &par, + &precip, &vpd, &vpdSoil, &vPress, &wspd, &soilWetness); + if (status == expectedNumCols + 1) { + printf("WARNING: ignoring location column in %s\n", climFile); + ++expectedNumCols; + hasLoc = 1; + } else { + status = fscanf(in, // NOLINT + "%d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", &year, + &day, &time, &length, &tair, &tsoil, &par, &precip, &vpd, + &vpdSoil, &vPress, &wspd, &soilWetness); + } + free(firstLine); - if (status != NUM_CLIM_FILE_COLS) { + if (status != expectedNumCols) { printf("Error reading climate file: bad data on first line\n"); exit(EXIT_CODE_INPUT_FILE_ERROR); } - if (loc != 0) { - printf("Error reading climate file: first location must be loc. 0\n"); - exit(EXIT_CODE_INPUT_FILE_ERROR); - } + firstClimate = (ClimateNode *)malloc(sizeof(ClimateNode)); + next = firstClimate; - firstClimates = - (ClimateNode **)malloc(numLocs * sizeof(ClimateNode *)); // vector of - // heads - for (currLoc = 0; currLoc < numLocs; currLoc++) { - firstClimates[currLoc] = NULL; // default is null - } - // if firstClimates[i] stays null for some i, that means we didn't read any - // climate data for location i, so use climate from location 0 - - currLoc = loc; - // allocate space for head - firstClimates[currLoc] = (ClimateNode *)malloc(sizeof(ClimateNode)); - next = firstClimates[currLoc]; - count = 0; while (status != EOF) { // we have another day's climate curr = next; - count++; // # of time steps in this location curr->year = year; curr->day = day; @@ -688,223 +673,172 @@ int *readClimData(char *climFile, int numLocs) { lastYear = year; - status = fscanf(in, "%d %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", - &loc, &year, &day, &time, &length, &tair, &tsoil, &par, - &precip, &vpd, &vpdSoil, &vPress, &wspd, &soilWetness); - + if (hasLoc) { + status = + fscanf(in, // NOLINT + "%d %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", + &dummyLoc, &year, &day, &time, &length, &tair, &tsoil, &par, + &precip, &vpd, &vpdSoil, &vPress, &wspd, &soilWetness); + } else { + status = fscanf(in, // NOLINT + "%d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", + &year, &day, &time, &length, &tair, &tsoil, &par, &precip, + &vpd, &vpdSoil, &vPress, &wspd, &soilWetness); + } if (status != EOF) { // we have another climate record - check new location, compare with old // location (currLoc), make sure new location is valid, and act // accordingly if (status != NUM_CLIM_FILE_COLS) { - printf("Error reading climate file: bad data near loc %d year %d" - " day %d\n", - loc, year, day); + printf("Error reading climate file: bad data near year %d day %d\n", + year, day); exit(EXIT_CODE_INPUT_FILE_ERROR); } - - if (loc == currLoc) { - // still reading climate records from the same place: add a node at end - // of linked list - next = (ClimateNode *)malloc(sizeof(ClimateNode)); - curr->nextClim = next; - // set this down here rather than at top of loop so head treated the - // same as rest of list - } else if (loc >= numLocs) { // (loc == numLocs is an error since we use - // 0-indexing) - printf("Error reading climate file: trying to read location %d, but " - "numLocs = %d\n", - loc, numLocs); - exit(1); - } else if (loc > currLoc) { - // we've advanced to the next location (note: possible that we skipped - // some locations: this is OK) - curr->nextClim = NULL; // terminate last linked list - steps[currLoc] = count; // record the number of time steps for last - // location - lastYear = -1; - count = 0; // reset count of # of time steps -#if GDD - gdd = 0; // reset growing degree days -#endif - currLoc = loc; - firstClimates[currLoc] = - (ClimateNode *)malloc(sizeof(ClimateNode)); // allocate space for - // head - next = firstClimates[currLoc]; // we'll start writing to next - // linked list - } else { // loc < currLoc - printf("Error reading climate file: was reading location %d, trying to " - "read location %d\n", - currLoc, loc); - printf("Climate records for a given location should be contiguous, and " - "locations should be in ascending order\n"); - exit(1); + // Check for older file with multiple locations - that's an error now + if (dummyLoc != firstLoc) { + printf("Error reading climate file: multiple locations not " + "supported\n"); + exit(EXIT_CODE_INPUT_FILE_ERROR); } + + // still reading climate records from the same place: add a node at end + // of linked list + next = (ClimateNode *)malloc(sizeof(ClimateNode)); + curr->nextClim = next; + // set this down here rather than at top of loop so head treated the + // same as rest of list } else { // status == EOF - no more records curr->nextClim = NULL; // terminate this last linked list - steps[currLoc] = count; // record the number of time steps for last - // location } - } // end while fclose(in); - - for (i = 0; i < numLocs; i++) { - if (steps[i] == 0) { // nothing read for this location - steps[i] = steps[0]; // this location will duplicate location 0 - } - } - return steps; } -// allocate & initialize spatialParamsPtr (a pointer to a SpatialParams pointer -// to allow for allocation) read in parameter file, put parameters (and other -// info) in spatialParams return numLocs (read from first line of spatialParams -// file) - -// note that the last argument in the "initializeOneSpatialParam" function -// indicates whether the given parameter is required -// 1 -> must be in param file, 0 -> optional -// parameters that are required for the model to run should be flagged as 1 to -// allow for error checking if a parameter is only required with certain -// options, this argument can be set to a boolean value that evaluates to true -// iff these options are enabled for maximum convenience, all parameters can be -// flagged as 0 (i.e. optional), but you are taking your life into your own -// hands if you do so -int readParamData(SpatialParams **spatialParamsPtr, char *paramFile, - char *spatialParamFile) { - FILE *paramF, *spatialParamF; - SpatialParams *spatialParams; // a local variable to prevent lots of - // unnecessary dereferences - int numLocs; - +/*! + * Read initial model parameter values from param file + * + * Note that the last argument in the "initializeOneModelParam" function + * indicates whether the given parameter is required: + * 1 -> must be in param file + * 0 -> optional + * + * @param modelParamsPtr ModelParams struct, will be alloc'd here + * @param paramFile Name of parameter file + */ +void readParamData(ModelParams **modelParamsPtr, const char *paramFile) { + FILE *paramF; + ModelParams *modelParams; // to prevent lots of + // unnecessary dereferences paramF = openFile(paramFile, "r"); - spatialParamF = openFile(spatialParamFile, "r"); - - int numRead = fscanf(spatialParamF, "%d", &numLocs); - - if (numLocs < 1 || numRead != 1) { - printf("Error: numLocs must be >= 1: read %d\n", numLocs); - exit(1); - } - *spatialParamsPtr = newSpatialParams(NUM_PARAMS, numLocs); - spatialParams = *spatialParamsPtr; // to prevent lots of unnecessary - // dereferences + *modelParamsPtr = newModelParams(NUM_PARAMS); + modelParams = *modelParamsPtr; // to prevent lots of unnecessary dereferences // clang-format off // NOLINTBEGIN - initializeOneSpatialParam(spatialParams, "plantWoodInit", &(params.plantWoodInit), 1); - initializeOneSpatialParam(spatialParams, "laiInit", &(params.laiInit), 1); - initializeOneSpatialParam(spatialParams, "litterInit", &(params.litterInit), 1); - initializeOneSpatialParam(spatialParams, "soilInit", &(params.soilInit), 1); - initializeOneSpatialParam(spatialParams, "litterWFracInit", &(params.litterWFracInit), 1); - initializeOneSpatialParam(spatialParams, "soilWFracInit", &(params.soilWFracInit), 1); - initializeOneSpatialParam(spatialParams, "snowInit", &(params.snowInit), 1); - initializeOneSpatialParam(spatialParams, "aMax", &(params.aMax), 1); - initializeOneSpatialParam(spatialParams, "aMaxFrac", &(params.aMaxFrac), 1); - initializeOneSpatialParam(spatialParams, "baseFolRespFrac", &(params.baseFolRespFrac), 1); - - initializeOneSpatialParam(spatialParams, "psnTMin", &(params.psnTMin), 1); - initializeOneSpatialParam(spatialParams, "psnTOpt", &(params.psnTOpt), 1); - initializeOneSpatialParam(spatialParams, "vegRespQ10", &(params.vegRespQ10), 1); - initializeOneSpatialParam(spatialParams, "growthRespFrac", &(params.growthRespFrac), GROWTH_RESP); - initializeOneSpatialParam(spatialParams, "frozenSoilFolREff", &(params.frozenSoilFolREff), 1); - initializeOneSpatialParam(spatialParams, "frozenSoilThreshold", &(params.frozenSoilThreshold), 1); - initializeOneSpatialParam(spatialParams, "dVpdSlope", &(params.dVpdSlope), 1); - initializeOneSpatialParam(spatialParams, "dVpdExp", &(params.dVpdExp), 1); - initializeOneSpatialParam(spatialParams, "halfSatPar", &(params.halfSatPar), 1); - initializeOneSpatialParam(spatialParams, "attenuation", &(params.attenuation), 1); - - initializeOneSpatialParam(spatialParams, "leafOnDay", &(params.leafOnDay), !((GDD) || (SOIL_PHENOL))); - initializeOneSpatialParam(spatialParams, "gddLeafOn", &(params.gddLeafOn), GDD); - initializeOneSpatialParam(spatialParams, "soilTempLeafOn", &(params.soilTempLeafOn), SOIL_PHENOL); - initializeOneSpatialParam(spatialParams, "leafOffDay", &(params.leafOffDay), 1); - initializeOneSpatialParam(spatialParams, "leafGrowth", &(params.leafGrowth), 1); - initializeOneSpatialParam(spatialParams, "fracLeafFall", &(params.fracLeafFall), 1); - initializeOneSpatialParam(spatialParams, "leafAllocation", &(params.leafAllocation), 1); - initializeOneSpatialParam(spatialParams, "leafTurnoverRate", &(params.leafTurnoverRate), 1); - initializeOneSpatialParam(spatialParams, "baseVegResp", &(params.baseVegResp), 1); - initializeOneSpatialParam(spatialParams, "litterBreakdownRate", &(params.litterBreakdownRate), LITTER_POOL); - - initializeOneSpatialParam(spatialParams, "fracLitterRespired", &(params.fracLitterRespired), LITTER_POOL); - initializeOneSpatialParam(spatialParams, "baseSoilResp", &(params.baseSoilResp), 1); - initializeOneSpatialParam(spatialParams, "baseSoilRespCold", &(params.baseSoilRespCold), SEASONAL_R_SOIL); - initializeOneSpatialParam(spatialParams, "soilRespQ10", &(params.soilRespQ10), 1); - initializeOneSpatialParam(spatialParams, "soilRespQ10Cold", &(params.soilRespQ10Cold), SEASONAL_R_SOIL); - initializeOneSpatialParam(spatialParams, "coldSoilThreshold", &(params.coldSoilThreshold), SEASONAL_R_SOIL); - - initializeOneSpatialParam(spatialParams, "E0", &(params.E0), LLOYD_TAYLOR); - initializeOneSpatialParam(spatialParams, "T0", &(params.T0), LLOYD_TAYLOR); - initializeOneSpatialParam(spatialParams, "soilRespMoistEffect", &(params.soilRespMoistEffect), ((WATER_HRESP) && !(DAYCENT_WATER_HRESP))); - initializeOneSpatialParam(spatialParams, "waterRemoveFrac", &(params.waterRemoveFrac), 1); - initializeOneSpatialParam(spatialParams, "frozenSoilEff", &(params.frozenSoilEff), 1); - initializeOneSpatialParam(spatialParams, "wueConst", &(params.wueConst), 1); - initializeOneSpatialParam(spatialParams, "litterWHC", &(params.litterWHC), 1); - initializeOneSpatialParam(spatialParams, "soilWHC", &(params.soilWHC), 1); - initializeOneSpatialParam(spatialParams, "immedEvapFrac", &(params.immedEvapFrac), COMPLEX_WATER); - initializeOneSpatialParam(spatialParams, "fastFlowFrac", &(params.fastFlowFrac), COMPLEX_WATER); - initializeOneSpatialParam(spatialParams, "leafPoolDepth", &(params.leafPoolDepth), LEAF_WATER); - - initializeOneSpatialParam(spatialParams, "snowMelt", &(params.snowMelt), SNOW); - initializeOneSpatialParam(spatialParams, "litWaterDrainRate", &(params.litWaterDrainRate), LITTER_WATER_DRAINAGE); - initializeOneSpatialParam(spatialParams, "rdConst", &(params.rdConst), (COMPLEX_WATER) || (PENMAN_MONTEITH_TRANS)); - initializeOneSpatialParam(spatialParams, "rSoilConst1", &(params.rSoilConst1), COMPLEX_WATER); - initializeOneSpatialParam(spatialParams, "rSoilConst2", &(params.rSoilConst2), COMPLEX_WATER); - initializeOneSpatialParam(spatialParams, "leafCSpWt", &(params.leafCSpWt), 1); - initializeOneSpatialParam(spatialParams, "cFracLeaf", &(params.cFracLeaf), 1); - initializeOneSpatialParam(spatialParams, "woodTurnoverRate", &(params.woodTurnoverRate), 1); - initializeOneSpatialParam(spatialParams, "qualityLeaf", &(params.qualityLeaf), SOIL_QUALITY); - initializeOneSpatialParam(spatialParams, "qualityWood", &(params.qualityWood), SOIL_QUALITY); - - initializeOneSpatialParam(spatialParams, "efficiency", &(params.efficiency), (SOIL_QUALITY) || (MICROBES)); - initializeOneSpatialParam(spatialParams, "maxIngestionRate", &(params.maxIngestionRate), (SOIL_QUALITY) || (MICROBES)); - initializeOneSpatialParam(spatialParams, "halfSatIngestion", &(params.halfSatIngestion), MICROBES); - initializeOneSpatialParam(spatialParams, "totNitrogen", &(params.totNitrogen), STOICHIOMETRY); - initializeOneSpatialParam(spatialParams, "microbeNC", &(params.microbeNC), STOICHIOMETRY); - initializeOneSpatialParam(spatialParams, "microbeInit", &(params.microbeInit), (SOIL_QUALITY) || (MICROBES)); - initializeOneSpatialParam(spatialParams, "fineRootFrac", &(params.fineRootFrac), ROOTS); - initializeOneSpatialParam(spatialParams, "coarseRootFrac", &(params.coarseRootFrac), ROOTS); - - initializeOneSpatialParam(spatialParams, "fineRootAllocation", &(params.fineRootAllocation), ROOTS); - initializeOneSpatialParam(spatialParams, "woodAllocation", &(params.woodAllocation), ROOTS); - - initializeOneSpatialParam(spatialParams, "fineRootExudation", &(params.fineRootExudation), ROOTS); - initializeOneSpatialParam(spatialParams, "coarseRootExudation", &(params.coarseRootExudation), ROOTS); - initializeOneSpatialParam(spatialParams, "fineRootTurnoverRate", &(params.fineRootTurnoverRate), ROOTS); - initializeOneSpatialParam(spatialParams, "coarseRootTurnoverRate", &(params.coarseRootTurnoverRate), ROOTS); - initializeOneSpatialParam(spatialParams, "baseFineRootResp", &(params.baseFineRootResp), ROOTS); - initializeOneSpatialParam(spatialParams, "baseCoarseRootResp", &(params.baseCoarseRootResp), ROOTS); - initializeOneSpatialParam(spatialParams, "fineRootQ10", &(params.fineRootQ10), ROOTS); - initializeOneSpatialParam(spatialParams, "coarseRootQ10", &(params.coarseRootQ10), ROOTS); - - initializeOneSpatialParam(spatialParams, "baseMicrobeResp", &(params.baseMicrobeResp), MICROBES); - initializeOneSpatialParam(spatialParams, "microbeQ10", &(params.microbeQ10), MICROBES); - initializeOneSpatialParam(spatialParams, "microbePulseEff", &(params.microbePulseEff), (ROOTS) && (MICROBES) ); - initializeOneSpatialParam(spatialParams, "m_ballBerry", &(params.m_ballBerry), 1); + initializeOneModelParam(modelParams, "plantWoodInit", &(params.plantWoodInit), 1); + initializeOneModelParam(modelParams, "laiInit", &(params.laiInit), 1); + initializeOneModelParam(modelParams, "litterInit", &(params.litterInit), 1); + initializeOneModelParam(modelParams, "soilInit", &(params.soilInit), 1); + initializeOneModelParam(modelParams, "litterWFracInit", &(params.litterWFracInit), 1); + initializeOneModelParam(modelParams, "soilWFracInit", &(params.soilWFracInit), 1); + initializeOneModelParam(modelParams, "snowInit", &(params.snowInit), 1); + initializeOneModelParam(modelParams, "aMax", &(params.aMax), 1); + initializeOneModelParam(modelParams, "aMaxFrac", &(params.aMaxFrac), 1); + initializeOneModelParam(modelParams, "baseFolRespFrac", &(params.baseFolRespFrac), 1); + + initializeOneModelParam(modelParams, "psnTMin", &(params.psnTMin), 1); + initializeOneModelParam(modelParams, "psnTOpt", &(params.psnTOpt), 1); + initializeOneModelParam(modelParams, "vegRespQ10", &(params.vegRespQ10), 1); + initializeOneModelParam(modelParams, "growthRespFrac", &(params.growthRespFrac), GROWTH_RESP); + initializeOneModelParam(modelParams, "frozenSoilFolREff", &(params.frozenSoilFolREff), 1); + initializeOneModelParam(modelParams, "frozenSoilThreshold", &(params.frozenSoilThreshold), 1); + initializeOneModelParam(modelParams, "dVpdSlope", &(params.dVpdSlope), 1); + initializeOneModelParam(modelParams, "dVpdExp", &(params.dVpdExp), 1); + initializeOneModelParam(modelParams, "halfSatPar", &(params.halfSatPar), 1); + initializeOneModelParam(modelParams, "attenuation", &(params.attenuation), 1); + + initializeOneModelParam(modelParams, "leafOnDay", &(params.leafOnDay), !((GDD) || (SOIL_PHENOL))); + initializeOneModelParam(modelParams, "gddLeafOn", &(params.gddLeafOn), GDD); + initializeOneModelParam(modelParams, "soilTempLeafOn", &(params.soilTempLeafOn), SOIL_PHENOL); + initializeOneModelParam(modelParams, "leafOffDay", &(params.leafOffDay), 1); + initializeOneModelParam(modelParams, "leafGrowth", &(params.leafGrowth), 1); + initializeOneModelParam(modelParams, "fracLeafFall", &(params.fracLeafFall), 1); + initializeOneModelParam(modelParams, "leafAllocation", &(params.leafAllocation), 1); + initializeOneModelParam(modelParams, "leafTurnoverRate", &(params.leafTurnoverRate), 1); + initializeOneModelParam(modelParams, "baseVegResp", &(params.baseVegResp), 1); + initializeOneModelParam(modelParams, "litterBreakdownRate", &(params.litterBreakdownRate), LITTER_POOL); + + initializeOneModelParam(modelParams, "fracLitterRespired", &(params.fracLitterRespired), LITTER_POOL); + initializeOneModelParam(modelParams, "baseSoilResp", &(params.baseSoilResp), 1); + initializeOneModelParam(modelParams, "baseSoilRespCold", &(params.baseSoilRespCold), SEASONAL_R_SOIL); + initializeOneModelParam(modelParams, "soilRespQ10", &(params.soilRespQ10), 1); + initializeOneModelParam(modelParams, "soilRespQ10Cold", &(params.soilRespQ10Cold), SEASONAL_R_SOIL); + initializeOneModelParam(modelParams, "coldSoilThreshold", &(params.coldSoilThreshold), SEASONAL_R_SOIL); + + initializeOneModelParam(modelParams, "E0", &(params.E0), LLOYD_TAYLOR); + initializeOneModelParam(modelParams, "T0", &(params.T0), LLOYD_TAYLOR); + initializeOneModelParam(modelParams, "soilRespMoistEffect", &(params.soilRespMoistEffect), ((WATER_HRESP) && !(DAYCENT_WATER_HRESP))); + initializeOneModelParam(modelParams, "waterRemoveFrac", &(params.waterRemoveFrac), 1); + initializeOneModelParam(modelParams, "frozenSoilEff", &(params.frozenSoilEff), 1); + initializeOneModelParam(modelParams, "wueConst", &(params.wueConst), 1); + initializeOneModelParam(modelParams, "litterWHC", &(params.litterWHC), 1); + initializeOneModelParam(modelParams, "soilWHC", &(params.soilWHC), 1); + initializeOneModelParam(modelParams, "immedEvapFrac", &(params.immedEvapFrac), COMPLEX_WATER); + initializeOneModelParam(modelParams, "fastFlowFrac", &(params.fastFlowFrac), COMPLEX_WATER); + initializeOneModelParam(modelParams, "leafPoolDepth", &(params.leafPoolDepth), LEAF_WATER); + + initializeOneModelParam(modelParams, "snowMelt", &(params.snowMelt), SNOW); + initializeOneModelParam(modelParams, "litWaterDrainRate", &(params.litWaterDrainRate), LITTER_WATER_DRAINAGE); + initializeOneModelParam(modelParams, "rdConst", &(params.rdConst), (COMPLEX_WATER) || (PENMAN_MONTEITH_TRANS)); + initializeOneModelParam(modelParams, "rSoilConst1", &(params.rSoilConst1), COMPLEX_WATER); + initializeOneModelParam(modelParams, "rSoilConst2", &(params.rSoilConst2), COMPLEX_WATER); + initializeOneModelParam(modelParams, "leafCSpWt", &(params.leafCSpWt), 1); + initializeOneModelParam(modelParams, "cFracLeaf", &(params.cFracLeaf), 1); + initializeOneModelParam(modelParams, "woodTurnoverRate", &(params.woodTurnoverRate), 1); + initializeOneModelParam(modelParams, "qualityLeaf", &(params.qualityLeaf), SOIL_QUALITY); + initializeOneModelParam(modelParams, "qualityWood", &(params.qualityWood), SOIL_QUALITY); + + initializeOneModelParam(modelParams, "efficiency", &(params.efficiency), (SOIL_QUALITY) || (MICROBES)); + initializeOneModelParam(modelParams, "maxIngestionRate", &(params.maxIngestionRate), (SOIL_QUALITY) || (MICROBES)); + initializeOneModelParam(modelParams, "halfSatIngestion", &(params.halfSatIngestion), MICROBES); + initializeOneModelParam(modelParams, "totNitrogen", &(params.totNitrogen), STOICHIOMETRY); + initializeOneModelParam(modelParams, "microbeNC", &(params.microbeNC), STOICHIOMETRY); + initializeOneModelParam(modelParams, "microbeInit", &(params.microbeInit), (SOIL_QUALITY) || (MICROBES)); + initializeOneModelParam(modelParams, "fineRootFrac", &(params.fineRootFrac), ROOTS); + initializeOneModelParam(modelParams, "coarseRootFrac", &(params.coarseRootFrac), ROOTS); + + initializeOneModelParam(modelParams, "fineRootAllocation", &(params.fineRootAllocation), ROOTS); + initializeOneModelParam(modelParams, "woodAllocation", &(params.woodAllocation), ROOTS); + + initializeOneModelParam(modelParams, "fineRootExudation", &(params.fineRootExudation), ROOTS); + initializeOneModelParam(modelParams, "coarseRootExudation", &(params.coarseRootExudation), ROOTS); + initializeOneModelParam(modelParams, "fineRootTurnoverRate", &(params.fineRootTurnoverRate), ROOTS); + initializeOneModelParam(modelParams, "coarseRootTurnoverRate", &(params.coarseRootTurnoverRate), ROOTS); + initializeOneModelParam(modelParams, "baseFineRootResp", &(params.baseFineRootResp), ROOTS); + initializeOneModelParam(modelParams, "baseCoarseRootResp", &(params.baseCoarseRootResp), ROOTS); + initializeOneModelParam(modelParams, "fineRootQ10", &(params.fineRootQ10), ROOTS); + initializeOneModelParam(modelParams, "coarseRootQ10", &(params.coarseRootQ10), ROOTS); + + initializeOneModelParam(modelParams, "baseMicrobeResp", &(params.baseMicrobeResp), MICROBES); + initializeOneModelParam(modelParams, "microbeQ10", &(params.microbeQ10), MICROBES); + initializeOneModelParam(modelParams, "microbePulseEff", &(params.microbePulseEff), (ROOTS) && (MICROBES) ); + initializeOneModelParam(modelParams, "m_ballBerry", &(params.m_ballBerry), 1); // NOLINTEND // clang-format on - readSpatialParams(spatialParams, paramF, spatialParamF); + readModelParams(modelParams, paramF); fclose(paramF); - fclose(spatialParamF); - - return numLocs; } -// pre: out is open for writing -// print header line to output file - -// I wish someone who write a .header file to automatically output the correct -// header - -// Not only that ...I'd like the options used in the model run to be added to -// the file; - +/*! + * Print header row to output file + * + * @param out File pointer for output + */ void outputHeader(FILE *out) { fprintf(out, "Notes: (PlantWoodC, PlantLeafC, Soil and Litter in g C/m^2; " "Water and Snow in cm; SoilWetness is fraction of WHC;\n"); @@ -929,12 +863,17 @@ void outputHeader(FILE *out) { "evapotranspiration fluxestranspiration fPAR\n"); } -// pre: out is open for writing -// print current state to output file -void outputState(FILE *out, int loc, int year, int day, double time) { +/*! + * Print current state values to output file + * @param out File pointer for output + * @param year + * @param day + * @param time + */ +void outputState(FILE *out, int year, int day, double time) { - fprintf(out, "%8d %4d %3d %5.2f %8.2f %8.2f ", loc, year, day, time, - envi.plantWoodC, envi.plantLeafC); + fprintf(out, "%4d %3d %5.2f %8.2f %8.2f ", year, day, time, envi.plantWoodC, + envi.plantLeafC); #if SOIL_MULTIPOOL int counter; @@ -969,41 +908,28 @@ void outputState(FILE *out, int loc, int year, int day, double time) { } // de-allocate space used for climate linked list -void freeClimateList(int numLocs) { +void freeClimateList() { ClimateNode *curr, *prev; - int loc; - - for (loc = 0; loc < numLocs; loc++) { // loop through firstClimates, - // deallocating each linked list - curr = firstClimates[loc]; - while (curr != NULL) { - prev = curr; - curr = curr->nextClim; - free(prev); - } + + curr = firstClimate; + while (curr != NULL) { + prev = curr; + curr = curr->nextClim; + free(prev); } - // and finally deallocate the vector itself: - free(firstClimates); } #if EVENT_HANDLER // de-allocate space used for events linked list -void freeEventList(int numLocs) { +void freeEventList() { EventNode *curr, *prev; - int loc; - - for (loc = 0; loc < numLocs; loc++) { - // loop through events, deallocating each linked list - curr = events[loc]; - while (curr != NULL) { - prev = curr; - curr = curr->nextEvent; - free(prev); - } - } - // and finally deallocate the vector itself: - free(events); + curr = events; + while (curr != NULL) { + prev = curr; + curr = curr->nextEvent; + free(prev); + } } #endif @@ -1012,7 +938,7 @@ void freeEventList(int numLocs) { // rather than returning a value, they have as parameters the variable(s) which // they modify so a single function can modify multiple variables -/** +/*! * @brief Compute canopy light effect using Simpson's rule. * * Similar to light attenuation in PnET, first calculate light @@ -1131,7 +1057,7 @@ void calcLightEff(double *lightEff, double lai, double par) { // day^-1) and base foliar respiration without temp, water, etc. (g C * m^-2 // ground area * day^-1) void potPsn(double *potGrossPsn, double *baseFolResp, double lai, double tair, - double vpd, double par, int day) { + double vpd, double par, int /*day*/) { double grossAMax; // maximum possible gross respiration (nmol CO2 * g^-1 leaf // * sec^-1) double dTemp, dVpd, lightEff; // decrease in photosynth. due to temp, vpd and @@ -1162,95 +1088,97 @@ void potPsn(double *potGrossPsn, double *baseFolResp, double lai, double tair, *potGrossPsn = grossAMax * dTemp * dVpd * lightEff * conversion; *baseFolResp = respPerGram * conversion; // do foliar resp. even if no // photosynthesis in this time step - - // printf("%f %f %f %f %f %f ", climate->length, grossAMax, conversion, dTemp, - // dVpd, lightEff); } -void moisture_bwb(double *trans, double *dWater, double potGrossPsn, double vpd, - double vPress, double plantLeafC, double leafCSpWt, - double soilWater) { - /*moisture_bwb: Calculates moisture use by using a Ball Woodrow Berry - * Instead of A (net photosynthesis) we are driving the Ball Woodrow Berry - * equation with potGrossPsn units: g C * m^-2 ground area * day^-1 - * formulation which estimates gs from A */ - double potTrans; // potential transpiration in the absense of plant water - // stress (cm H20 * day^-1) - double removableWater; - double gs_canopy; // gs stomatal conductance... canopy level - double CO2_stom = 380; - double vPress_sat; // Saturating Vapor Press (sum of VPD and vPress) - double RH_pcent; // Relative Humidity - the ratio of vPress to saturating - // vapor pressure - double lai_int; // calculate lai from current plantLeafC and the leafCSpwt */ - lai_int = envi.plantLeafC / params.leafCSpWt; - /* We do not require lai since potGrossPsn is in units of - * g C * m^-2 ground area * day^-1 */ - vPress_sat = climate->vPress + climate->vpd; // calculate saturating vapor - // pressure - RH_pcent = climate->vPress / vPress_sat; // calculate relative humidity for - // ball berry - - // double wue; // water use efficiency, in mg CO2 fixed * g^-1 H20 transpired - - if (potGrossPsn < TINY) { // avoid divide by 0 - *trans = 0.0; // no photosynthesis -> no transpiration - *dWater = 1; // dWater doesn't matter, since we don't have any - // photosynthesis - } else { - /*Ball Berry Equation - * unit issues: - * potGrossPsn should be converted to - * micro mol c02 per m^2 leaf area - * there are 1/12 mol of carbon in a gram - * potGrossPsn/12*1000*LAI = umol carbon dioxide per - * m^2 - * */ - - gs_canopy = params.m_ballBerry * potGrossPsn / 12 * 1000 * lai_int * - (RH_pcent / CO2_stom); - - // gcan = m*A*RelHum/CO2; - // mol / m^2 leaf area - - potTrans = (gs_canopy * climate->vpd) / lai_int * 20 / 1000; - /* gs_canopy*climate_>vpd is the transpiration rate mol/m^2 leaf area per - * day dividing by lai_init converts to mol per m^2 ground area there are - * 20g water per mol multiply by the density of water 1cm^3/g multipy by - * 100cm and divide by 10^6cm^3 - * - * must convert to cm3 per cm2 land area. - */ - removableWater = soilWater * params.waterRemoveFrac; - if (climate->tsoil < params.frozenSoilThreshold) { - // frozen soil - less or no water available - /* frozen soil effect: fraction of water available if soil is frozen - * (assume amt. of water avail. w/ frozen soil scales linearly with amt. - * of water avail. in thawed soil) */ - removableWater *= params.frozenSoilEff; - } - if (removableWater >= potTrans) { - *trans = potTrans; - } else { - *trans = removableWater; - } - -#if WATER_PSN // we're modeling water stress - *dWater = *trans / potTrans; // from PnET: equivalent to setting DWATER_MAX - // = 1 -#else // WATER_PSN = 0 - if (climate->tsoil < params.frozenSoilThreshold && - params.frozenSoilEff == 0) - // (note: can't have partial shutdown of psn with frozen soil if WATER_PSN - // = 0) - *dWater = 0; // still allow total shut down of psn. if soil is frozen - else // either soil is thawed, or frozenSoilEff > 0 - *dWater = 1; // no water stress, even if *trans/potTrans < 1 -#endif // WATER_PSN - // printf("Remove %f potT %f dW %f vpd %f \n", removableWater, potTrans, - // *dWater, climate->vpd); - } -} +// Unused - remove? +// void moisture_bwb(double *trans, double *dWater, double potGrossPsn, double +// vpd, +// double vPress, double plantLeafC, double leafCSpWt, +// double soilWater) { +// /*moisture_bwb: Calculates moisture use by using a Ball Woodrow Berry +// * Instead of A (net photosynthesis) we are driving the Ball Woodrow Berry +// * equation with potGrossPsn units: g C * m^-2 ground area * day^-1 +// * formulation which estimates gs from A */ +// double potTrans; // potential transpiration in the absense of plant water +// // stress (cm H20 * day^-1) +// double removableWater; +// double gs_canopy; // gs stomatal conductance... canopy level +// double CO2_stom = 380; +// double vPress_sat; // Saturating Vapor Press (sum of VPD and vPress) +// double RH_pcent; // Relative Humidity - the ratio of vPress to saturating +// // vapor pressure +// double lai_int; // calculate lai from current plantLeafC and the leafCSpwt +// */ lai_int = envi.plantLeafC / params.leafCSpWt; +// /* We do not require lai since potGrossPsn is in units of +// * g C * m^-2 ground area * day^-1 */ +// vPress_sat = climate->vPress + climate->vpd; // calculate saturating vapor +// // pressure +// RH_pcent = climate->vPress / vPress_sat; // calculate relative humidity for +// // ball berry +// +// // double wue; // water use efficiency, in mg CO2 fixed * g^-1 H20 +// transpired +// +// if (potGrossPsn < TINY) { // avoid divide by 0 +// *trans = 0.0; // no photosynthesis -> no transpiration +// *dWater = 1; // dWater doesn't matter, since we don't have any +// // photosynthesis +// } else { +// /*Ball Berry Equation +// * unit issues: +// * potGrossPsn should be converted to +// * micro mol c02 per m^2 leaf area +// * there are 1/12 mol of carbon in a gram +// * potGrossPsn/12*1000*LAI = umol carbon dioxide per +// * m^2 +// * */ +// +// gs_canopy = params.m_ballBerry * potGrossPsn / 12 * 1000 * lai_int * +// (RH_pcent / CO2_stom); +// +// // gcan = m*A*RelHum/CO2; +// // mol / m^2 leaf area +// +// potTrans = (gs_canopy * climate->vpd) / lai_int * 20 / 1000; +// /* gs_canopy*climate_>vpd is the transpiration rate mol/m^2 leaf area per +// * day dividing by lai_init converts to mol per m^2 ground area there are +// * 20g water per mol multiply by the density of water 1cm^3/g multipy by +// * 100cm and divide by 10^6cm^3 +// * +// * must convert to cm3 per cm2 land area. +// */ +// removableWater = soilWater * params.waterRemoveFrac; +// if (climate->tsoil < params.frozenSoilThreshold) { +// // frozen soil - less or no water available +// /* frozen soil effect: fraction of water available if soil is frozen +// * (assume amt. of water avail. w/ frozen soil scales linearly with amt. +// * of water avail. in thawed soil) */ +// removableWater *= params.frozenSoilEff; +// } +// if (removableWater >= potTrans) { +// *trans = potTrans; +// } else { +// *trans = removableWater; +// } +// +// #if WATER_PSN // we're modeling water stress +// *dWater = *trans / potTrans; // from PnET: equivalent to setting +// DWATER_MAX +// // = 1 +// #else // WATER_PSN = 0 +// if (climate->tsoil < params.frozenSoilThreshold && +// params.frozenSoilEff == 0) +// // (note: can't have partial shutdown of psn with frozen soil if +// WATER_PSN +// // = 0) +// *dWater = 0; // still allow total shut down of psn. if soil is frozen +// else // either soil is thawed, or frozenSoilEff > 0 +// *dWater = 1; // no water stress, even if *trans/potTrans < 1 +// #endif // WATER_PSN +// // printf("Remove %f potT %f dW %f vpd %f \n", removableWater, potTrans, +// // *dWater, climate->vpd); +// } +//} /* it would be preferable to include CO2 concentration as an input variable in * the climate file @@ -1264,95 +1192,101 @@ void moisture_bwb(double *trans, double *dWater, double potGrossPsn, double vpd, * we will assume this is 370 and I may edit the Climate input file to include * CO2 in future. BALL_BERRY_m is hardcoded into the equation at 3.89 */ -// printf("%d\n", gs_canopy); - // Penman Monteith method of estimating transiration and water use // Started Nov 2006 - coding commenced Nov 20th -void moisture_pm(double *trans, double *dWater, double potGrossPsn, double vpd, - double soilWater) { - double potTrans; // potential transpiration in the absense of plant water - // stress (cm H20 * day^-1) - double removableWater; - // not used double wue; // water use efficiency, in mg CO2 fixed * g^-1 H20 - // transpired - double gapfraction = 17; // gap fraction = 17% - the Penman Monteith Equation - // doesn't appear to be sensitive to gapfraction; - double rCanConst = 30.8; // Rcan = rCanConst/windspeed - 30.8 ; - // is the median rCanConst to achieve the measured Rcan of 8.5 s/m; - double DELTA; // Delta is calculated as d Vsat/ d Tair; - - /* - * required constants and other values - * Rn // Net radiation // read in from clim file - * G // soil heat flux // read in from clim file - * RHO // DENSITY OF (DRY) AIR (KG M-3) - * CP // SPECIFIC HEAT OF AIR AT CONST PRESSURE (J KG-1 K-1) - * VPD // VAPOR PRESSURE DEFICIT (PA) - * rb - BULK AERODYNAMIC RESISTANCE (S/M) - * rc - CANOPY AERODYNAMIC RESISTANCE (S/M) - * LAMDA //LATENT HEAT OF VAPORIZATION OF WATER (J Kg-1) - 2450 @20C - * GAMMA // PSYCHROMETRIC CONSTANT Kg M-2 S - * DELTA slope of the saturation vapour pressure vs Temperature relationship - */ - if (potGrossPsn < TINY) { // avoid divide by 0 - *trans = 0.0; // no photosynthesis -> no transpiration - *dWater = 1; // dWater doesn't matter, since we don't have any - // photosynthesis - } else { - DELTA = (2508.3 / ((climate->tair + 237.3) * (climate->tair + 237.3)) * - exp((17.3 * climate->tair) / (climate->tair + 237.3))); - potTrans = - (DELTA * ((1 - gapfraction / 100) * climate->par - climate->tsoil) + - (RHO * CP * vpd) / (rCanConst / climate->wspd)) / - (DELTA + GAMMA * (1 + (params.rdConst / climate->wspd) / - (rCanConst / climate->wspd))); - - /* - * the aerodynamic resistance - of the canopy can be calculated as follows: - * rc = [ln((Zm-d)/Zom)*ln((Zh-d)/Zoh)] / k*k*wspd ; - * - * rc aerodynamic resistance [s m-1], - * Zm height of wind measurements [m], - * Zh height of humidity measurements [m], - * d zero plane displacement height [m], - * Zom roughness length governing momentum transfer [m], - * Zoh roughness length governing transfer of heat and vapour [m], - * k von Karman's constant, 0.41 [-], - * wspd wind speed at height z [m s-1]. - */ - removableWater = soilWater * params.waterRemoveFrac; - if (climate->tsoil < params.frozenSoilThreshold) { - // frozen soil - less or no water available - removableWater *= params.frozenSoilEff; - } - /* frozen soil effect: fraction of water available if soil is frozen (assume - * amt. of water avail. w/ frozen soil scales linearly with amt. of water - * avail. in thawed soil) */ - if (removableWater >= potTrans) { - *trans = potTrans; - } else { - *trans = removableWater; - } - -#if WATER_PSN // we're modeling water stress - *dWater = *trans / potTrans; // from PnET: equivalent to setting DWATER_MAX - // = 1 -#else // WATER_PSN = 0 - if (climate->tsoil < params.frozenSoilThreshold && - params.frozenSoilEff == 0) - // (note: can't have partial shutdown of psn with frozen soil if WATER_PSN - // = 0) - *dWater = 0; // still allow total shut down of psn. if soil is frozen - else // either soil is thawed, or frozenSoilEff > 0 - *dWater = 1; // no water stress, even if *trans/potTrans < 1 -#endif // WATER_PSN - // printf("PotGrossPsn: %f dWater %f potTrans %f\n", potGrossPsn, *dWater, - // potTrans); - } - - // printf("%f\n", *dWater); -} +// Unused - remove? +// void moisture_pm(double *trans, double *dWater, double potGrossPsn, double +// vpd, +// double soilWater) { +// double potTrans; // potential transpiration in the absense of plant water +// // stress (cm H20 * day^-1) +// double removableWater; +// // not used double wue; // water use efficiency, in mg CO2 fixed * g^-1 +// H20 +// // transpired +// double gapfraction = 17; // gap fraction = 17% - the Penman Monteith +// Equation +// // doesn't appear to be sensitive to gapfraction; +// double rCanConst = 30.8; // Rcan = rCanConst/windspeed - 30.8 ; +// // is the median rCanConst to achieve the measured Rcan of 8.5 s/m; +// double DELTA; // Delta is calculated as d Vsat/ d Tair; +// +// /* +// * required constants and other values +// * Rn // Net radiation // read in from clim file +// * G // soil heat flux // read in from clim file +// * RHO // DENSITY OF (DRY) AIR (KG M-3) +// * CP // SPECIFIC HEAT OF AIR AT CONST PRESSURE (J KG-1 K-1) +// * VPD // VAPOR PRESSURE DEFICIT (PA) +// * rb - BULK AERODYNAMIC RESISTANCE (S/M) +// * rc - CANOPY AERODYNAMIC RESISTANCE (S/M) +// * LAMDA //LATENT HEAT OF VAPORIZATION OF WATER (J Kg-1) - 2450 @20C +// * GAMMA // PSYCHROMETRIC CONSTANT Kg M-2 S +// * DELTA slope of the saturation vapour pressure vs Temperature relationship +// */ +// if (potGrossPsn < TINY) { // avoid divide by 0 +// *trans = 0.0; // no photosynthesis -> no transpiration +// *dWater = 1; // dWater doesn't matter, since we don't have any +// // photosynthesis +// } else { +// DELTA = (2508.3 / ((climate->tair + 237.3) * (climate->tair + 237.3)) * +// exp((17.3 * climate->tair) / (climate->tair + 237.3))); +// potTrans = +// (DELTA * ((1 - gapfraction / 100) * climate->par - climate->tsoil) + +// (RHO * CP * vpd) / (rCanConst / climate->wspd)) / +// (DELTA + GAMMA * (1 + (params.rdConst / climate->wspd) / +// (rCanConst / climate->wspd))); +// +// /* +// * the aerodynamic resistance - of the canopy can be calculated as +// follows: +// * rc = [ln((Zm-d)/Zom)*ln((Zh-d)/Zoh)] / k*k*wspd ; +// * +// * rc aerodynamic resistance [s m-1], +// * Zm height of wind measurements [m], +// * Zh height of humidity measurements [m], +// * d zero plane displacement height [m], +// * Zom roughness length governing momentum transfer [m], +// * Zoh roughness length governing transfer of heat and vapour [m], +// * k von Karman's constant, 0.41 [-], +// * wspd wind speed at height z [m s-1]. +// */ +// removableWater = soilWater * params.waterRemoveFrac; +// if (climate->tsoil < params.frozenSoilThreshold) { +// // frozen soil - less or no water available +// removableWater *= params.frozenSoilEff; +// } +// /* frozen soil effect: fraction of water available if soil is frozen +// (assume +// * amt. of water avail. w/ frozen soil scales linearly with amt. of water +// * avail. in thawed soil) */ +// if (removableWater >= potTrans) { +// *trans = potTrans; +// } else { +// *trans = removableWater; +// } +// +// #if WATER_PSN // we're modeling water stress +// *dWater = *trans / potTrans; // from PnET: equivalent to setting +// DWATER_MAX +// // = 1 +// #else // WATER_PSN = 0 +// if (climate->tsoil < params.frozenSoilThreshold && +// params.frozenSoilEff == 0) +// // (note: can't have partial shutdown of psn with frozen soil if +// WATER_PSN +// // = 0) +// *dWater = 0; // still allow total shut down of psn. if soil is frozen +// else // either soil is thawed, or frozenSoilEff > 0 +// *dWater = 1; // no water stress, even if *trans/potTrans < 1 +// #endif // WATER_PSN +// // printf("PotGrossPsn: %f dWater %f potTrans %f\n", potGrossPsn, *dWater, +// // potTrans); +// } +// +// // printf("%f\n", *dWater); +//} // calculate transpiration (cm H20 * day^-1) // and dWater (factor between 0 and 1) @@ -1806,7 +1740,7 @@ void calcRootResp(double *rootResp, double respQ10, double baseRate, // calculate foliar resp., wood maint. resp. and growth resp., all in g C * m^-2 // ground area * day^-1 growth resp. modeled in a very simple way void vegResp2(double *folResp, double *woodResp, double *growthResp, - double baseFolResp, double gpp) { + double baseFolResp, double /*gpp*/) { *folResp = baseFolResp * pow(params.vegRespQ10, (climate->tair - params.psnTOpt) / 10.0); if (climate->tsoil < params.frozenSoilThreshold) { @@ -1945,10 +1879,11 @@ void calcMaintenanceRespiration(double tsoil, double water, double whc) { #endif } -double microbeQualityEfficiency(double soilQuality) { - - return params.efficiency; // Efficiency an increasing function of quality -} +// Unused - remove? +// double microbeQualityEfficiency(double soilQuality) { +// +// return params.efficiency; // Efficiency an increasing function of quality +// } void microbeGrowth(void) { #if SOIL_MULTIPOOL @@ -2432,20 +2367,6 @@ void updateTrackers(double oldSoilWater) { trackers.yearlyNee = 0.0; lastYear = climate->year; - - // At start of 1999, reset cumulative trackers - // Note that this is only for one specific application: we don't usually - // want to do this - /* - if (climate->year == 1999) { - trackers.totGpp = 0.0; - trackers.totRtot = 0.0; - trackers.totRa = 0.0; - trackers.totRh = 0.0; - trackers.totNpp = 0.0; - trackers.totNee = 0.0; - } - */ } trackers.gpp = fluxes.photosynthesis * climate->length; @@ -2513,7 +2434,7 @@ void updateTrackers(double oldSoilWater) { * \brief Process events for current location/year/day * * For a given year and day (as determined by the global `climate` - * pointer), process all events listed in the global `locEvents` pointer for the + * pointer), process all events listed in the global `events` pointer for the * referenced location. * * For each event, modify state variables according to the model for that event, @@ -2525,7 +2446,7 @@ void processEvents(void) { // This should be in events.h/c, but with all the global state defined in this // file, let's leave it here for now. Maybe someday we will factor that out. - // If locEvent starts off NULL, this function will just fall through, as it + // If event starts off NULL, this function will just fall through, as it // should. const int climYear = climate->year; const int climDay = climate->day; @@ -2534,18 +2455,17 @@ void processEvents(void) { // be in chrono order. However, we need to check to make sure the current // event is not in the past, as that would indicate an event that did not have // a corresponding climate file record. - while (locEvent != NULL && locEvent->year <= climYear && - locEvent->day <= climDay) { - if (locEvent->year < climYear || locEvent->day < climDay) { - printf("Agronomic event found for loc: %d year: %d day: %d that does not " + while (event != NULL && event->year <= climYear && event->day <= climDay) { + if (event->year < climYear || event->day < climDay) { + printf("Agronomic event found for year: %d day: %d that does not " "have a corresponding record in the climate file\n", - locEvent->loc, locEvent->year, locEvent->day); + event->year, event->day); exit(EXIT_CODE_INPUT_FILE_ERROR); } - switch (locEvent->type) { + switch (event->type) { // Implementation TBD, as we enable the various event types case IRRIGATION: { - const IrrigationParams *irrParams = locEvent->eventParams; + const IrrigationParams *irrParams = event->eventParams; const double amount = irrParams->amountAdded; double soilAmount, evapAmount; if (irrParams->method == CANOPY) { @@ -2564,11 +2484,11 @@ void processEvents(void) { } fluxes.immedEvap += evapAmount; envi.soilWater += soilAmount; - writeEventOut(eventOutFile, locEvent, 2, "envi.soilWater", soilAmount, + writeEventOut(eventOutFile, event, 2, "envi.soilWater", soilAmount, "fluxes.immedEvap", evapAmount); } break; case PLANTING: { - const PlantingParams *plantParams = locEvent->eventParams; + const PlantingParams *plantParams = event->eventParams; const double leafC = plantParams->leafC; const double woodC = plantParams->woodC; const double fineRootC = plantParams->fineRootC; @@ -2582,13 +2502,13 @@ void processEvents(void) { // FUTURE: allocate to N pools - writeEventOut(eventOutFile, locEvent, 4, "envi.plantLeafC", leafC, + writeEventOut(eventOutFile, event, 4, "envi.plantLeafC", leafC, "envi.plantWoodC", woodC, "envi.fineRootC", fineRootC, "envi.coarseRootC", coarseRootC); } break; case HARVEST: { // Harvest can both remove biomass and move biomass to the litter pool - const HarvestParams *harvParams = locEvent->eventParams; + const HarvestParams *harvParams = event->eventParams; const double fracRA = harvParams->fractionRemovedAbove; const double fracTA = harvParams->fractionTransferredAbove; const double fracRB = harvParams->fractionRemovedBelow; @@ -2614,7 +2534,7 @@ void processEvents(void) { // FUTURE: move/remove biomass in N pools - writeEventOut(eventOutFile, locEvent, 5, "env.litter", litterAdd, + writeEventOut(eventOutFile, event, 5, "env.litter", litterAdd, "envi.plantLeafC", leafDelta, "envi.plantWoodC", woodDelta, "envi.fineRootC", fineDelta, "envi.coarseRootC", coarseDelta); @@ -2628,11 +2548,11 @@ void processEvents(void) { printf("Fertilization events not yet implemented\n"); break; default: - printf("Unknown event type (%d) in processEvents()\n", locEvent->type); + printf("Unknown event type (%d) in processEvents()\n", event->type); exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM); } - locEvent = locEvent->nextEvent; + event = event->nextEvent; } } #endif @@ -2758,11 +2678,7 @@ void initPhenologyTrackers(void) { // Setup model to run at given location (0-indexing: if only one location, loc // should be 0) -void setupModel(SpatialParams *spatialParams, int loc) { - - // load parameters into global param structure: spatialParams was told where - // to put values in readParamData - loadSpatialParams(spatialParams, loc); +void setupModel(void) { // a test: use constant (measured) soil respiration: // make it so soil resp. is 5.2 g C m-2 day-1 at 10 degrees C, @@ -2862,12 +2778,8 @@ void setupModel(SpatialParams *spatialParams, int loc) { envi.snow = params.snowInit; - if (firstClimates[loc] != NULL) { - climate = firstClimates[loc]; // set climate ptr to point to first climate - // record in this location - } else { // no climate data for this location - climate = firstClimates[0]; // use climate data from location 0 - } + climate = firstClimate; + initTrackers(); initPhenologyTrackers(); resetMeanTracker(meanNPP, 0); // initialize with mean NPP (over last @@ -2880,57 +2792,32 @@ void setupModel(SpatialParams *spatialParams, int loc) { // Setup events at given location #if EVENT_HANDLER -void setupEvents(int currLoc) { locEvent = events[currLoc]; } +void setupEvents() { event = events; } #endif -/* Do one run of the model using parameter values in spatialParams - If out != NULL, output results to out - If printHeader = 1, print a header for the output file, if 0 don't - If outputItems != NULL, do additional outputting as given by this structure - (1 variable per file) If loc == -1, then print currLoc as first item on each - line Run at spatial location given by loc (0-indexing) - or run everywhere if - loc = -1 Note: number of locations given in spatialParams -*/ -void runModelOutput(FILE *out, OutputItems *outputItems, int printHeader, - SpatialParams *spatialParams, int loc) { - int firstLoc, lastLoc, currLoc; - char label[64]; - +// See sipnet.h +void runModelOutput(FILE *out, OutputItems *outputItems, int printHeader) { if ((out != NULL) && printHeader) { outputHeader(out); } - if (loc == -1) { // run everywhere - firstLoc = 0; - lastLoc = spatialParams->numLocs - 1; - } else { // just run at one point - firstLoc = lastLoc = loc; - } - - for (currLoc = firstLoc; currLoc <= lastLoc; currLoc++) { - setupModel(spatialParams, currLoc); + setupModel(); #if EVENT_HANDLER - setupEvents(currLoc); + setupEvents(); #endif - if ((loc == -1) && (outputItems != NULL)) { // print the current location - // at the start of the line - sprintf(label, "%d", currLoc); - writeOutputItemLabels(outputItems, label); - } - while (climate != NULL) { - updateState(); - if (out != NULL) { - outputState(out, currLoc, climate->year, climate->day, climate->time); - } - if (outputItems != NULL) { - writeOutputItemValues(outputItems); - } - climate = climate->nextClim; + while (climate != NULL) { + updateState(); + if (out != NULL) { + outputState(out, climate->year, climate->day, climate->time); } if (outputItems != NULL) { - terminateOutputItemLines(outputItems); + writeOutputItemValues(outputItems); } + climate = climate->nextClim; + } + if (outputItems != NULL) { + terminateOutputItemLines(outputItems); } } @@ -2976,11 +2863,7 @@ void printModelComponents(FILE *out) { fprintf(out, "\n"); } -/* PRE: outputItems has been created with newOutputItems - - Setup outputItems structure - Each variable added will be output in a separate file ('*.varName') - */ +// See sipnet.h void setupOutputItems(OutputItems *outputItems) { addOutputItem(outputItems, "NEE", &(trackers.nee)); addOutputItem(outputItems, "NEE_cum", &(trackers.totNee)); @@ -2988,61 +2871,33 @@ void setupOutputItems(OutputItems *outputItems) { addOutputItem(outputItems, "GPP_cum", &(trackers.totGpp)); } -/* do initializations that only have to be done once for all model runs: - read in climate data and initial parameter values - parameter values get stored in spatialParams (along with other parameter - information), which gets allocated and initialized here (thus requiring that - spatialParams is passed in as a pointer to a pointer) number of time steps in - each location gets stored in steps vector, which gets dynamically allocated - with malloc (steps must be a pointer to a pointer so it can be malloc'ed) - - also set up pointers to different output data types - and setup meanNPP tracker - - initModel returns number of spatial locations - - paramFile is parameter data file - paramFile-spatial is file with parameter values of spatially-varying - parameters (1st line contains number of locations) climFile is climate data - file -*/ -int initModel(SpatialParams **spatialParams, int **steps, char *paramFile, - char *climFile) { - char spatialParamFile[256]; - int numLocs; - - strcpy(spatialParamFile, paramFile); - strcat(spatialParamFile, "-spatial"); - - numLocs = readParamData(spatialParams, paramFile, spatialParamFile); - *steps = readClimData(climFile, numLocs); +// See sipnet.h +void initModel(ModelParams **modelParams, const char *paramFile, + const char *climFile) { + readParamData(modelParams, paramFile); + readClimData(climFile); meanNPP = newMeanTracker(0, MEAN_NPP_DAYS, MEAN_NPP_MAX_ENTRIES); meanGPP = newMeanTracker(0, MEAN_GPP_SOIL_DAYS, MEAN_GPP_SOIL_MAX_ENTRIES); meanFPAR = newMeanTracker(0, MEAN_FPAR_DAYS, MEAN_FPAR_MAX_ENTRIES); - return numLocs; } -/* Do initialization of event data if event handling is turned on. - * Populates static event structs - */ #if EVENT_HANDLER -void initEvents(char *eventFile, int numLocs, int printHeader) { - events = readEventData(eventFile, numLocs); +// See sipnet.h +void initEvents(char *eventFile, int printHeader) { + events = readEventData(eventFile); eventOutFile = openEventOutFile(printHeader); } #endif -// call this when done running model: -// de-allocates space for climate and event linked lists -// (needs to know number of locations) -void cleanupModel(int numLocs) { - freeClimateList(numLocs); +// See sipnet.h +void cleanupModel() { + freeClimateList(); deallocateMeanTracker(meanNPP); deallocateMeanTracker(meanGPP); deallocateMeanTracker(meanFPAR); #if EVENT_HANDLER - freeEventList(numLocs); + freeEventList(); closeEventOutFile(eventOutFile); #endif } diff --git a/src/sipnet/sipnet.h b/src/sipnet/sipnet.h index 8f0b111cc..a2942e2bb 100644 --- a/src/sipnet/sipnet.h +++ b/src/sipnet/sipnet.h @@ -4,7 +4,7 @@ #define SIPNET_H #include -#include "common/spatialParams.h" +#include "common/modelParams.h" #include "outputItems.h" // write to file which model components are turned on @@ -12,61 +12,53 @@ // pre: out is open for writing void printModelComponents(FILE *out); -/* do initializations that only have to be done once for all model runs: - read in climate data and initial parameter values - parameter values get stored in spatialParams (along with other parameter - information), which gets allocated and initialized here (thus requiring that - spatialParams is passed in as a pointer to a pointer) number of time steps in - each location gets stored in steps vector, which gets dynamically allocated - with malloc (steps must be a pointer to a pointer so it can be malloc'ed) - - also set up pointers to different output data types - and setup meanNPP tracker - - initModel returns number of spatial locations - - paramFile is parameter data file - paramFile-spatial is file with parameter values of spatially-varying - parameters (1st line contains number of locations) - climFile is climate data file -*/ -int initModel(SpatialParams **spatialParams, int **steps, char *paramFile, - char *climFile); +/*! + * Do model initializations + * + * Read in initial parameter values and climate data. Also set up pointers to + * different output data types and setup meanNPP tracker. + * + * @param modelParams pointer to ModelParams struct, will be alloc'd here + * @param paramFile name of parameter file + * @param climFile name of climate file + */ +void initModel(ModelParams **modelParams, const char *paramFile, + const char *climFile); -/* +/*! * Read in event data for all the model runs * * Read in event data from a file with the following specification: * - one line per event - * - all events are ordered first by location (ascending) and then by year/day - * (ascending) + * - all events are ordered by year/day ascending * * @param eventFile Name of file containing event data - * @param numLocs Number of locations in the event file. */ -void initEvents(char *eventFile, int numLocs, int printHeader); - -// call this when done running model: -// de-allocates space for climate linked list -// (needs to know number of locations) -void cleanupModel(int numLocs); +void initEvents(char *eventFile, int printHeader); -/* Do one run of the model using parameter values in spatialParams - If out != NULL, output results to out - If printHeader = 1, print a header for the output file, if 0 don't - If outputItems != NULL, do additional outputting as given by this - structure (1 variable per file) - If loc == -1, then print currLoc as first item on each line - Run at spatial location given by loc (0-indexing) - or run everywhere if loc - = -1 Note: number of locations given in spatialParams -*/ -void runModelOutput(FILE *out, OutputItems *outputItems, int printHeader, - SpatialParams *spatialParams, int loc); +/*! + * Free allocated memory + * + * Call this when done running model + */ +void cleanupModel(void); -/* PRE: outputItems has been created with newOutputItems +/*! Run the model using parameter values in modelParams + * + * @param out File pointer for main output; can be null to suppress this + * @param outputItems OutputItems struct used for individual output param files + * Can be null to suppress this output + * @param printHeader Whether to print a header row in output and events.out + */ +void runModelOutput(FILE *out, OutputItems *outputItems, int printHeader); +/*! Setup outputItems structure + Each variable added will be output in a separate file ('*.varName') + + @param outputItems OutputItems struct to be filled out; must be created with + newOutputItems() */ void setupOutputItems(OutputItems *outputItems); diff --git a/src/sipnet/sipnet.in b/src/sipnet/sipnet.in index 8711a04e9..9c711360b 100644 --- a/src/sipnet/sipnet.in +++ b/src/sipnet/sipnet.in @@ -13,7 +13,7 @@ FILENAME = Sites/Niwot/niwot ! FILENAME.clim is the file of climate data for each time step LOCATION = 0 -! Location to run at (-1 means run at all locations) +! LOCATION is obsolete; this value is ignored DO_MAIN_OUTPUT = 1 ! If 1, do the primary output to FILENAME.out (time series of all output @@ -28,5 +28,3 @@ DO_SINGLE_OUTPUTS = 0 ! If 1, do extra outputs: one variable per file (e.g. FILENAME.NEE) ! If 0, don't do these extra outputs ! (These extra outputs are setup in sipnet.c : setupOutputItems) -! For a spatial run (LOCATION = -1), the first value on each line is -! the location, followed by the time series From 4f6e5bb7e4f6c96c13cc1ca011d83ab4bc5840ea Mon Sep 17 00:00:00 2001 From: Alomir Date: Fri, 6 Jun 2025 09:21:38 -0400 Subject: [PATCH 02/18] Updates for removal of multi-site support --- docs/CHANGELOG.md | 3 ++ docs/parameters.md | 108 +++++++++++++++++++++++---------------------- 2 files changed, 58 insertions(+), 53 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 94081d52f..41b96c03b 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -46,6 +46,9 @@ sections to include in release notes: - Removed many experimental sites, data, and executable code as part of reorg (#34, #37) - Removed obsolete run types senstest and montecarlo and associated code (#69, #76) - Removed obsolete estimate program and associated code (#70, #82) +- Removed multi-site support: + - Output files no longer have a location column + - Input climate files with location as the first column can still be read, though a warning will be thrown. Climate files with multiple location will error. ### Git SHA [TBD] diff --git a/docs/parameters.md b/docs/parameters.md index b873e15d8..a079238af 100644 --- a/docs/parameters.md +++ b/docs/parameters.md @@ -409,7 +409,10 @@ For each step of the model, the following inputs are needed. These are provided | 12 | wspd | avg. wind speed | m/s | | | 13 | soilWetness | fractional soil wetness | unitless (0-1) | $f_\text{WHC}$; Used if `MODEL_WATER=0`; if `MODEL_WATER=1`, soil wetness is simulated| -Note: An older format for this file included location as the first column. +Note: An older format for this file included location as the first column. Files with this older format can still be read by sipnet: +* SIPNET will print a warning indicating that it is ignoring the location column +* If there is more than one location specified in the file, SIPNET will error and halt + #### Example `sipnet.clim` file: Column names are not used, but are: @@ -442,18 +445,17 @@ For managed ecosystems, the following inputs are provided in a file named `event | col | parameter | description | units | notes | |-------|-------------|------------------------------------|-------------|---------------------------------------| -| 1 | loc | spatial location index | | maps to param-spatial file | -| 2 | year | year of start of this timestep | | e.g. 2010 | -| 3 | day | day of start of this timestep | Day of year | 1 = Jan 1 | -| 4 | event_type | type of event | | one of plant, harv, till, fert, irrig | -| 5...n | event_param | parameter associated with event | | see table below | +| 1 | year | year of start of this timestep | | e.g. 2010 | +| 2 | day | day of start of this timestep | Day of year | 1 = Jan 1 | +| 3 | event_type | type of event | | one of plant, harv, till, fert, irrig | +| 4...n | event_param | parameter associated with event | | see table below | - Agronomic events are stored in `events.in`, one event per row -- Events in the file are sorted first by location, and then chronologically +- Events in the file must be sorted chronologically - Events are specified by year and day (no hourly timestamp) -- It is assumed that there is one (or more) records in the climate file for each location/day that appears in the events file +- It is assumed that there is one (or more) records in the climate file for each year/day that appears in the events file - SIPNET will throw an error if it finds an event with no corresponding climate record -- Events are processed with the first climate record that occurs for the relevant location/day as an instantaneous one-time change +- Events are processed with the first climate record that occurs for the relevant year/day as an instantaneous one-time change - We may need events with duration later, spec TBD. Tillage is likely in this bucket. - The effects of an event are applied after fluxes are calculated for the current climate record; they are applied as a delta to one or more state variables, as required @@ -534,11 +536,11 @@ Notes: #### Example of `events.in` file: ``` -1 2022 40 irrig 5 1 # 5cm canopy irrigation on day 40 applied to soil -1 2022 40 fert 0 0 10 # fertilized with 10 g / m2 N_min on day 40 of 2022 -1 2022 35 till 0.2 0.3 # tilled on day 35, soil organic matter pool decomposition rate increases by 20% and soil litter pool decomposition rate increases by 30% -1 2022 50 plant 10 3 2 5 # plant emergence on day 50 with 10/3/2/4 g C / m2, respectively, added to the leaf/wood/fine root/coarse root pools -1 2022 250 harv 0.1 # harvest 10% of aboveground plant biomass on day 250 +2022 40 irrig 5 1 # 5cm canopy irrigation on day 40 applied to soil +2022 40 fert 0 0 10 # fertilized with 10 g / m2 N_min on day 40 of 2022 +2022 35 till 0.2 0.3 # tilled on day 35, soil organic matter pool decomposition rate increases by 20% and soil litter pool decomposition rate increases by 30% +2022 50 plant 10 3 2 5 # plant emergence on day 50 with 10/3/2/4 g C / m2, respectively, added to the leaf/wood/fine root/coarse root pools +2022 250 harv 0.1 # harvest 10% of aboveground plant biomass on day 250 ``` #### Events output @@ -546,57 +548,57 @@ Notes: SIPNET will create a file named `events.out` when event handling is enabled. -This file will have one row for each agronomic event that is processed. Each row lists location, -time, event type, and parameter name/value pairs for all state variables that the event +This file will have one row for each agronomic event that is processed. Each row lists year, +day, event type, and parameter name/value pairs for all state variables that the event affects. Example events.out file below, with header enabled for clarity. Note the delimiters: spaces up to the param-values pairs, commas separating PV pairs, and `=` between param and value. ``` -loc year day type param_name=delta[,param_name=delta,...] -0 2024 65 plant envi.plantLeafC=3.00,envi.plantWoodC=4.00,envi.fineRootC=5.00,envi.coarseRootC=6.00 -0 2024 70 irrig envi.soilWater=5.00 -0 2024 200 harv env.litter=5.46,envi.plantLeafC=-5.93,envi.plantWoodC=-4.75,envi.fineRootC=-3.73,envi.coarseRootC=-3.89 -1 2024 65 plant envi.plantLeafC=3.00,envi.plantWoodC=5.00,envi.fineRootC=7.00,envi.coarseRootC=9.00 -1 2024 70 irrig fluxes.immedEvap=2.50,envi.soilWater=2.50 -1 2024 200 harv env.litter=4.25,envi.plantLeafC=-1.39,envi.plantWoodC=-1.63,envi.fineRootC=-2.52,envi.coarseRootC=-2.97 +year day type param_name=delta[,param_name=delta,...] +2023 65 plant envi.plantLeafC=3.00,envi.plantWoodC=4.00,envi.fineRootC=5.00,envi.coarseRootC=6.00 +2023 70 irrig envi.soilWater=5.00 +2023 200 harv env.litter=5.46,envi.plantLeafC=-5.93,envi.plantWoodC=-4.75,envi.fineRootC=-3.73,envi.coarseRootC=-3.89 +2024 65 plant envi.plantLeafC=3.00,envi.plantWoodC=5.00,envi.fineRootC=7.00,envi.coarseRootC=9.00 +2024 70 irrig fluxes.immedEvap=2.50,envi.soilWater=2.50 +2024 200 harv env.litter=4.25,envi.plantLeafC=-1.39,envi.plantWoodC=-1.63,envi.fineRootC=-2.52,envi.coarseRootC=-2.97 ``` _Note: `events.out` logs all parameters changed by an event for debugging and testing purposes; For downstream analyses that only need the date and event type, `events.in` is equivalent and easier to parse._ + ## Outputs -| | Symbol | Parameter Name | Definition | Units | -| -- | ----------- | ------------------ | ------------------------------ | ----------- | -| 1 | |loc | spatial location index | | -| 2 | |year | year of start of this timestep | | -| 3 | |day | day of start of this timestep | | -| 4 | |time | time of start of this timestep | | -| 5 | |plantWoodC | carbon in wood | g C/m$^2$ | -| 6 | |plantLeafC | carbon in leaves | g C/m$^2$ | -| 7 | |soil | carbon in mineral soil | g C/m$^2$ | -| 8 | |microbeC | carbon in soil microbes | g C/m$^2$ | -| 9 | |coarseRootC | carbon in coarse roots | g C/m$^2$ | -| 10 | |fineRootC | carbon in fine roots | g C/m$^2$ | -| 11 | |litter | carbon in litter | g C/m$^2$ | -| 12 | |litterWater | moisture in litter layer | cm | -| 13 | |soilWater | moisture in soil | cm | -| 14 |$f_\text{WHC}$|soilWetnessFrac | moisture in soil as fraction | | -| 15 | |snow | snow water | cm | -| 16 | |npp | net primary production | g C/m$^2$ | -| 17 | |nee | net ecosystem production | g C/m$^2$ | -| 18 | |cumNEE | cumulative nee | g C/m$^2$ | -| 19 | $GPP$ |gpp | gross ecosystem production | g C/m$^2$ | -| 20 | $R_{A,\text{above}}$ |rAboveground | plant respiration above ground | g C/m$^2$ | -| 21 | $R_H$ |rSoil | soil respiration | g C/m$^2$ | -| 22 | $R_{A\text{, root}}$ |rRoot | root respiration | g C/m$^2$ | -| 23 | $R$ |rtot | total respiration | g C/m$^2$ | -| 24 | |fluxestranspiration | transpiration | cm | -| | $F^N_\text{vol}$ |fluxesn2o | Nitrous Oxide flux | g N/m$^2$ / timestep | -| | $F^C_{\text{CH}_4}$ |fluxesch4 | Methane Flux | g C/m$^2$ / timestep | -| | $F^N_\text{vol}$ |fluxesn2o | Nitrous Oxide flux | g N/m$^2$ / timestep | -| | $F^C_{\text{CH}_4}$ |fluxesch4 | Methane Flux | g C/m$^2$ / timestep | +| | Symbol | Parameter Name | Definition | Units | +|----|----------------------|---------------------|--------------------------------|----------------------| +| 1 | | year | year of start of this timestep | | +| 2 | | day | day of start of this timestep | | +| 3 | | time | time of start of this timestep | | +| 4 | | plantWoodC | carbon in wood | g C/m$^2$ | +| 5 | | plantLeafC | carbon in leaves | g C/m$^2$ | +| 6 | | soil | carbon in mineral soil | g C/m$^2$ | +| 7 | | microbeC | carbon in soil microbes | g C/m$^2$ | +| 8 | | coarseRootC | carbon in coarse roots | g C/m$^2$ | +| 9 | | fineRootC | carbon in fine roots | g C/m$^2$ | +| 10 | | litter | carbon in litter | g C/m$^2$ | +| 11 | | litterWater | moisture in litter layer | cm | +| 12 | | soilWater | moisture in soil | cm | +| 13 | $f_\text{WHC}$ | soilWetnessFrac | moisture in soil as fraction | | +| 14 | | snow | snow water | cm | +| 15 | | npp | net primary production | g C/m$^2$ | +| 16 | | nee | net ecosystem production | g C/m$^2$ | +| 17 | | cumNEE | cumulative nee | g C/m$^2$ | +| 18 | $GPP$ | gpp | gross ecosystem production | g C/m$^2$ | +| 19 | $R_{A,\text{above}}$ | rAboveground | plant respiration above ground | g C/m$^2$ | +| 20 | $R_H$ | rSoil | soil respiration | g C/m$^2$ | +| 21 | $R_{A\text{, root}}$ | rRoot | root respiration | g C/m$^2$ | +| 22 | $R$ | rtot | total respiration | g C/m$^2$ | +| 23 | | fluxestranspiration | transpiration | cm | +| 24 | $F^N_\text{vol}$ | fluxesn2o | Nitrous Oxide flux | g N/m$^2$ / timestep | +| 25 | $F^C_{\text{CH}_4}$ | fluxesch4 | Methane Flux | g C/m$^2$ / timestep | +| 26 | $F^N_\text{vol}$ | fluxesn2o | Nitrous Oxide flux | g N/m$^2$ / timestep | +| 27 | $F^C_{\text{CH}_4}$ | fluxesch4 | Methane Flux | g C/m$^2$ / timestep |