📄 read_soilparam_arc.c
字号:
/** Get Snow Surface Roughness **/
fscanf(soilparam,"%s",tmpstr);
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
temp.snow_rough = read_arcinfo_value(namestr,temp.lat,temp.lng);
/** Get Average Annual Precipitation **/
fscanf(soilparam,"%s",tmpstr);
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
temp.annual_prec = read_arcinfo_value(namestr,temp.lat,temp.lng);
/** Get Layer Percent Sand **/
for(layer=0;layer<options.Nlayer;layer++) {
fscanf(soilparam,"%s",tmpstr);
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
sand[layer] = read_arcinfo_value(namestr,temp.lat,temp.lng);
}
/** Get Layer Percent Clay **/
for(layer=0;layer<options.Nlayer;layer++) {
fscanf(soilparam,"%s",tmpstr);
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
clay[layer] = read_arcinfo_value(namestr,temp.lat,temp.lng);
}
/** Get Layer Saturated Hydraulic Conductivity **/
for(layer=0;layer<options.Nlayer;layer++) {
fscanf(soilparam,"%s",tmpstr);
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
temp.Ksat[layer] = read_arcinfo_value(namestr,temp.lat,temp.lng);
}
/** Get Layer Soil Moisture Diffusion Parameter **/
for(layer=0;layer<options.Nlayer;layer++) {
fscanf(soilparam,"%s",tmpstr);
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
temp.phi_s[layer] = read_arcinfo_value(namestr,temp.lat,temp.lng);
}
/** Get Initial Layer Moisture **/
for(layer=0;layer<options.Nlayer;layer++) {
fscanf(soilparam,"%s",tmpstr);
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
temp.init_moist[layer] = read_arcinfo_value(namestr,temp.lat,temp.lng);
}
/** Get Layer Thickness **/
for(layer=0;layer<options.Nlayer;layer++) {
fscanf(soilparam,"%s",tmpstr);
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
temp.depth[layer] = read_arcinfo_value(namestr,temp.lat,temp.lng);
temp.depth[layer] = (float)(int)(temp.depth[layer] * 1000 + 0.5) / 1000;
if(temp.depth[layer] < MINSOILDEPTH) {
fprintf(stderr,"Model will not function with layer depth %f < %f m.\n",
temp.depth[layer],MINSOILDEPTH);
exit(0);
}
}
if(temp.depth[0] > temp.depth[1]) {
sprintf(ErrStr,"ERROR: Model will not function with layer %i depth (%f m) < layer %i depth (%f m).\n",
0,temp.depth[0],1,temp.depth[1]);
nrerror(ErrStr);
}
/** Get Layer Bulk Density **/
for(layer=0;layer<options.Nlayer;layer++) {
fscanf(soilparam,"%s",tmpstr);
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
temp.bulk_density[layer] = read_arcinfo_value(namestr,temp.lat,temp.lng);
}
/** Get Layer Porosity **/
for(layer=0;layer<options.Nlayer;layer++) {
fscanf(soilparam,"%s",tmpstr);
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
temp.porosity[layer] = read_arcinfo_value(namestr,temp.lat,temp.lng);
}
/** Activate Frozen Soil for Grid Cell **/
fscanf(soilparam,"%s",tmpstr);
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
temp.FS_ACTIVE = (char)read_arcinfo_value(namestr,temp.lat,temp.lng);
if (options.COMPUTE_TREELINE && (fscanf(soilparam,"%s",tmpstr)) != EOF) {
/** Get Avg July Air Temperature **/
strcpy(namestr,soilparamdir);
strcat(namestr,tmpstr);
temp.avgJulyAirTemp = read_arcinfo_value(namestr,temp.lat,temp.lng);
options.JULY_TAVG_SUPPLIED = TRUE;
}
/*************************************************
if ARNO_PARAMS == TRUE then convert the baseflow
parameters d1, d2, d3, d4 to Ds, Dsmax, Ws, and c. JA
*************************************************/
if(options.ARNO_PARAMS) {
layer = options.Nlayer-1;
temp.Dsmax = temp.Dsmax *
pow((double)(1./(temp.max_moist[layer]-temp.Ws)), -temp.c) +
temp.Ds * temp.max_moist[layer];
temp.Ds = temp.Ds * temp.Ws / temp.Dsmax;
temp.Ws = temp.Ws/temp.max_moist[layer];
}
/*******************************************
Compute Soil Layer Properties
*******************************************/
sum_depth = 0.;
for(layer = 0; layer < options.Nlayer; layer++) {
sum_depth += temp.depth[layer];
temp.bulk_density[layer] *= 1000.;
temp.soil_density[layer] = temp.bulk_density[layer]
/ (1.0 - temp.porosity[layer]);
temp.quartz[layer] = sand[layer] / 100.;
temp.max_moist[layer] = temp.depth[layer] * temp.porosity[layer] * 1000.;
temp.bubble[layer] = exp(5.3396738 + 0.1845038*clay[layer]
- 2.48394546*temp.porosity[layer]
- 0.00213853*pow(clay[layer],2.)
- 0.04356349*sand[layer]*temp.porosity[layer]
- 0.61745089*clay[layer]*temp.porosity[layer]
+ 0.00143598*pow(sand[layer],2.)
* pow(temp.porosity[layer],2.)
- 0.00855375*pow(clay[layer],2.)
* pow(temp.porosity[layer],2.)
- 0.00001282*pow(sand[layer],2.)*clay[layer]
+ 0.00895359*pow(clay[layer],2.)*temp.porosity[layer]
- 0.00072472*pow(sand[layer],2.)*temp.porosity[layer]
+ 0.00000540*pow(clay[layer],2.)*sand[layer]
+ 0.50028060*pow(temp.porosity[layer],2.)*clay[layer]);
temp.expt[layer] = exp(-0.7842831 + 0.0177544*sand[layer]
- 1.062498*temp.porosity[layer]
- 0.00005304*pow(sand[layer],2.)
- 0.00273493*pow(clay[layer],2.)
+ 1.11134946*pow(temp.porosity[layer],2.)
- 0.03088295*sand[layer]*temp.porosity[layer]
+ 0.00026587*pow(sand[layer],2.)
* pow(temp.porosity[layer],2.)
- 0.00610522*pow(clay[layer],2.)
* pow(temp.porosity[layer],2.)
- 0.00000235*pow(sand[layer],2.)*clay[layer]
+ 0.00798746*pow(clay[layer],2.)*temp.porosity[layer]
- 0.00674491*pow(temp.porosity[layer],2.)*clay[layer]);
temp.expt[layer] = 2. / temp.expt[layer] + 3.;
temp.resid_moist[layer] = - 0.0182482 + 0.00087269 * sand[layer]
+ 0.00513488 * clay[layer]
+ 0.02939286 * temp.porosity[layer]
- 0.00015395 * pow(clay[layer],2.)
- 0.00108270 * sand[layer] * temp.porosity[layer]
- 0.00018233 * pow(clay[layer],2.)
* pow(temp.porosity[layer],2.)
+ 0.00030703 * pow(clay[layer],2.0)
* temp.porosity[layer]
- 0.00235840 * pow(temp.porosity[layer],2.)
* clay[layer];
/** Check for valid values of generated parameters **/
if(temp.bubble[layer]<1.36) {
fprintf(stderr,"WARNING: estimated bubbling pressure too low (%f), resetting to minimum value (%f).\n",temp.bubble[layer],1.36);
temp.bubble[layer] = 1.36;
}
if(temp.bubble[layer]>187.2) {
fprintf(stderr,"WARNING: estimated bubbling pressure too high (%f), resetting to maximum value (%f).\n",temp.bubble[layer],187.2);
temp.bubble[layer] = 187.2;
}
if(temp.expt[layer] < 2. / 1.090 + 3.) {
fprintf(stderr,"WARNING: estimated exponential (expt) too low (%f), resetting to minimum value (%f).\n", temp.expt[layer], 2. / 1.090 + 3.);
temp.expt[layer] = 2. / 1.090 + 3.;
}
if(temp.expt[layer] > 2. / 0.037 + 3.) {
fprintf(stderr,"WARNING: estimated exponential (expt) too high (%f), resetting to maximum value (%f).\n",temp.expt[layer], 2. / 0.037 + 3.);
temp.expt[layer] = 2. / 0.037 + 3.;
}
if(temp.resid_moist[layer] < -0.038) {
fprintf(stderr,"WARNING: estimated residual soil moisture too low (%f), resetting to minimum value (%f).\n",temp.resid_moist[layer],-0.038);
temp.resid_moist[layer] = -0.038;
}
if(temp.resid_moist[layer] > 0.205) {
fprintf(stderr,"WARNING: estimated residual soil moisture too high (%f), resetting to maximum value (%f).\n",temp.resid_moist[layer],0.205);
temp.resid_moist[layer] = 0.205;
}
tmp_bubble += temp.bubble[layer];
}
for(layer=0;layer<options.Nlayer;layer++) temp.bubble[layer] = tmp_bubble/3.;
/*******************************************
Validate Initial Soil Layer Moisture Content
*******************************************/
if (!options.INIT_STATE) { // only do this if we're not getting initial moisture from a state file
for(layer = 0; layer < options.Nlayer; layer++) {
if(temp.init_moist[layer] > temp.max_moist[layer]) {
fprintf(stderr,"Initial soil moisture (%f mm) is greater than the maximum moisture (%f mm) for layer %i.\n\tResetting soil moisture to maximum.\n",
temp.init_moist[layer], temp.max_moist[layer], layer);
temp.init_moist[layer] = temp.max_moist[layer];
}
if(temp.init_moist[layer] < temp.resid_moist[layer] * temp.depth[layer] * 1000.) {
fprintf(stderr,"Initial soil moisture (%f mm) is less than calculated residual moisture (%f mm) for layer %i.\n\tResetting soil moisture to residual moisture.\n",
temp.init_moist[layer], temp.resid_moist[layer] * temp.depth[layer] * 1000., layer);
temp.init_moist[layer] = temp.resid_moist[layer] * temp.depth[layer] * 1000.;
}
}
}
/**********************************************
Compute Maximum Infiltration for Upper Layers
**********************************************/
if(options.Nlayer==2)
temp.max_infil = (1.0+temp.b_infilt)*temp.max_moist[0];
else
temp.max_infil = (1.0+temp.b_infilt)*(temp.max_moist[0]
+ temp.max_moist[1]);
/****************************************************************
Compute Soil Layer Critical and Wilting Point Moisture Contents
****************************************************************/
for(layer=0;layer<options.Nlayer;layer++) {
temp.Wcr[layer] = Wcr_FRACT[layer] * temp.max_moist[layer];
temp.Wpwp[layer] = Wpwp_FRACT[layer] * temp.max_moist[layer];
if(temp.Wpwp[layer] > temp.Wcr[layer]) {
sprintf(ErrStr,"Calculated wilting point moisture (%f mm) is greater than calculated critical point moisture (%f mm) for layer %i.\n\tIn the soil parameter file, Wpwp_FRACT MUST be <= Wcr_FRACT.\n",
temp.Wpwp[layer], temp.Wcr[layer], layer);
nrerror(ErrStr);
}
if(temp.Wpwp[layer] < temp.resid_moist[layer] * temp.depth[layer] * 1000.) {
sprintf(ErrStr,"Calculated wilting point moisture (%f mm) is less than calculated residual moisture (%f mm) for layer %i.\n\tIn the soil parameter file, Wpwp_FRACT MUST be >= resid_moist / (1.0 - bulk_density/soil_density).\n",
temp.Wpwp[layer], temp.resid_moist[layer] * temp.depth[layer] * 1000., layer);
nrerror(ErrStr);
}
}
/*************************************************
Determine Central Longitude of Current Time Zone
*************************************************/
temp.time_zone_lng = off_gmt * 360./24.;
}
else RUN[0] = 0;
/* Allocate Layer - Node fraction array */
temp.layer_node_fract = (float **)malloc((options.Nlayer+1)*sizeof(float *));
for(layer=0;layer<=options.Nlayer;layer++)
temp.layer_node_fract[layer]
= (float *)malloc(options.Nnode*sizeof(float));
return temp;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -