Hi,
I have the following program. It worked great until I added the For loop at the begining, and modified if((i=0)&&(doy <=1)). Now the output looks good for a couple iterations, but then goes into this endless output of day nan nan nan ~0.035 nan nan nan nan. Any suggestions on how to fix this problem? Many thanks!
The code:
Code:
#define SMA 0.00036 /* DIM scaling parameter for SM curve */
#define SMB 105.0 /* DIM scaling parameter for SM curve */
#define FMET 0.59 /* DIM metabolic fraction of NPP; Potter et al. 1993, CHANGE FOR DIFFERENT BIOMES */
#define FSTR 0.70 /* DIM slow SOC fraction of structural SOC pool; Ise and Moorcroft 2006 */
#define KSTR 0.40 /* DIM proportional rate K for structural SOC decomposition; Ise and Moorcroft 2006 */
#define KSLW 0.0093 /* DIM proportional rate K for slow SOC decomposition; Ise and Moorcroft 2006 */
#define FRAUT 0.53 /* DIM proportion of GPP represented by autotrophic respiration. This value is taken from Choudhury (2001). The origional model used Gifford et al. 2005,
** /* value of 0.46. CHANGE FOR DIFFERENT BIOMES */
#define ME 0.45 /* DIM Microbial efficiency; McGill et al. 1981; Potter et al. 1993 */
#define KOPT 0.0301 /* Optimal decomposition rate (1/d) for soil metabolic fraction; Ise and Moorcroft 2006 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <unistd.h>
int main()
{
int doy;
doy=0;
double soil_t; /* surface soil temperature; deg K */
double soil_m; /* surface soil moisture; % saturation */
double gpp; /* GPP; g C m-2 d-1 */
double npp;
double rmet;
double rstr;
double rslw;
double npp_ann_init; /* Initial NPP pool; g C m-2 yr-1 */
double soil_t_av; /* long-term avg surface soil T; deg K */
double soil_m_av; /* long-term avg JJA soil moisture; % saturation */
double met_c_init; /* initial soil metabolic C pool; g C m-2; Ise and Moorcroft, 2006 */
double struct_c_init; /* initial structural soil C pool; g C m-2; Ise and Moorcroft, 2006 */
double slow_c_init; /* initial slow soil C pool; g C m-2; Ise and Moorcroft, 2006 */
double r_aut; /* autotrophic respiration; g C m-2 d-1 */
double r_het; /* heterotrophic respiration; g C m-2 d-1 */
double dCmet;
double dCstr;
double dCslw;
double cmet;
double cstr;
double cslw;
double r_tot; /* total respiration; g C m-2 d-1 */
double k_met_opt; /* optimal soil metabolic C decomposition rate; d-1 */
double k_met_av; /* average soil metabolic C decomposition rate; d-1 */
double k_met; /* actual soil metabolic C decomp rate; d-1 */
double k_struct; /* soil structural C decomp rate; d-1 */
double k_slow; /* soil slow C decomp rate; d-1 */
double tmult; /* DIM 0-1, soil temperature rate response curve */
double tmult_av; /* DIM 0-1, soil temperature rate response under avg conditions */
double wmult; /* DIM 0-1, soil moisture rate response curve */
double wmult_av; /* DIM 0-1, soil moisture rate response under avg conditions */
int flag; /* DIM 0.0 or 1.0 for spinup or daily run, respectively */
/* 1.0 open data input file */
FILE *datain;
int i;
i=1;
for (i = 0 ; i < 2; i++)
{
/* datain =fopen("/home/jennifer.watts/Carbon_Models/TCF_model_modified/Testing_Spinup/", "r"); */
datain =fopen("tcf_spinup.txt", "r");
if (datain == NULL)
{
fprintf(stderr, "Can't open input file!\n");
exit(-1);
}
while(fscanf(datain,"%d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf",
&flag, &doy, &gpp, &soil_t, &soil_m, &soil_t_av, &soil_m_av, &npp_ann_init, &met_c_init, &struct_c_init, &slow_c_init) != EOF)
{
doy++;
/* address inital values for year 1, day 1, of spin-up*/
if((i=0) && (doy <= 1 ))
{
cmet = met_c_init;
cstr = struct_c_init;
cslw = slow_c_init;
rmet = 0;
rstr = 0;
rslw = 0;
}
/* 3.0 Environmental Constraint Calculations*/
k_met_opt = KOPT;
/* calculate soil temperature rate response (Lloyd and Taylor 1994) */
tmult = exp(308.56 * ((1.0 / 66.02) - (1.0 / ((soil_t) - 227.13))));
/* calculate soil moisture response (Monsi 1969, Bunnell and Tait 1974, Oberbauer et al. 1991, 1996) */
wmult = SMA * (SMB * soil_m - pow(soil_m, 2.0));
/* 4.0 Calculate soil CO2 input and output; soil C pool values */
/* calculate decomposition rate constants */
k_met = k_met_opt * tmult * wmult;
k_struct = KSTR * k_met;
k_slow = KSLW * k_met;
rmet = (k_met * cmet * (1.0 - ME));
rstr = (k_struct * cstr * (1.0 - ME));
rslw = (k_slow * cslw * (1.0 - ME));
r_het = (rmet + rstr + rslw);
r_aut = (gpp * FRAUT);
r_tot = (r_het + r_aut);
npp = (gpp - r_aut);
dCmet = (FMET * npp - rmet);
dCstr = (((1 - FMET) * npp) - rstr - FSTR * cstr);
dCslw = ((FSTR * cstr) - rslw);
cmet += (cmet + dCmet);
cstr += (cstr + dCstr);
cslw += (cslw + dCslw);
/* print out daily updating of carbon pools for spinup run */
printf("%d %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f\n", doy, cmet, cstr, cslw, npp, r_tot, dCmet, dCstr, dCslw);
}
fclose(datain);
}
}