diff --git a/src/glm_flow.c b/src/glm_flow.c index 2e22b2c..bbd2283 100644 --- a/src/glm_flow.c +++ b/src/glm_flow.c @@ -405,7 +405,7 @@ void do_single_outflow(AED_REAL HeightOfOutflow, AED_REAL flow, OutflowDataType * Loop through all outflows and process - return the difference between * * total volume before and after. * ******************************************************************************/ -AED_REAL do_outflows(int jday) +AED_REAL do_outflows(int jday, AED_REAL day_fraction) { int i; AED_REAL DrawHeight = -1; //# Height of withdraw [m from bottom] @@ -461,12 +461,12 @@ AED_REAL do_outflows(int jday) if (seepage) { if (seepage_rate>zero) { // Darcy's Law used, so input rate is hydraulic conductivity (m/day) x hydraulic head - SeepDraw = seepage_rate * Lake[surfLayer].Height * Lake[surfLayer].LayerArea * 0.95; + SeepDraw = seepage_rate * Lake[surfLayer].Height * Lake[surfLayer].LayerArea * 0.95 * day_fraction; } else { // Constant seepage assumed, so input rate is dh (m/day) // 0.95 added since the effective area of seeping is probably // a bit less than max area of water innundation??? - SeepDraw = -seepage_rate * Lake[surfLayer].LayerArea * 0.95; + SeepDraw = -seepage_rate * Lake[surfLayer].LayerArea * 0.95 * day_fraction; } do_single_outflow(0., SeepDraw, NULL); } @@ -482,7 +482,7 @@ AED_REAL do_outflows(int jday) * level. Add in stack volumes when calculating overflow. * * After overflow, delete stack volumes from the structure. * ******************************************************************************/ -AED_REAL do_overflow(int jday) +AED_REAL do_overflow(int jday, AED_REAL day_fraction) { AED_REAL VolSum = Lake[surfLayer].Vol1; AED_REAL DrawHeight = 0.; @@ -501,7 +501,7 @@ AED_REAL do_overflow(int jday) AED_REAL ovfl_Q, ovfl_dz; ovfl_dz = MAX( Lake[surfLayer].Height - CrestHeight, zero ); - ovfl_Q = 2./3. * crest_factor * pow(2*g,0.5) * crest_width * pow(ovfl_dz,1.5); + ovfl_Q = 2./3. * crest_factor * pow(2*g,0.5) * crest_width * pow(ovfl_dz,1.5) * day_fraction * iSecsPerDay; ovfl_Q = MIN( (VolSum - VolAtCrest) , ovfl_Q ); do_single_outflow(CrestHeight, ovfl_Q , NULL); diff --git a/src/glm_flow.h b/src/glm_flow.h index f538cc4..6506762 100644 --- a/src/glm_flow.h +++ b/src/glm_flow.h @@ -32,8 +32,8 @@ #include "glm.h" -AED_REAL do_outflows(int jday); -AED_REAL do_overflow(int jday); +AED_REAL do_outflows(int jday, AED_REAL day_fraction); +AED_REAL do_overflow(int jday, AED_REAL day_fraction); AED_REAL do_inflows(void); #endif diff --git a/src/glm_init.c b/src/glm_init.c index 94b7d5f..0d7b9c9 100644 --- a/src/glm_init.c +++ b/src/glm_init.c @@ -1569,8 +1569,8 @@ void initialise_lake(int namlst) Lake[i].Salinity = the_sals[i]; } - if (the_heights[num_heights-1] > CrestHeight) { - fprintf(stderr, " ERROR: maximum height is greater than crest level\n"); + if (the_heights[num_heights-1] > MaxHeight) { + fprintf(stderr, " ERROR: initial first height of %f is greater than maximum height of %f \n", the_heights[num_heights-1], MaxHeight); exit(1); } num_depths = num_heights; diff --git a/src/glm_model.c b/src/glm_model.c index 2cb005b..eb82811 100644 --- a/src/glm_model.c +++ b/src/glm_model.c @@ -274,6 +274,7 @@ void do_model(int jstart, int nsave) AED_REAL Elev[MaxInf]; int jday, ntot, stepnum, stoptime; int i, j; + AED_REAL day_fraction; /*------------------------------------------------------------------------*/ @@ -294,8 +295,10 @@ void do_model(int jstart, int nsave) read_daily_met(jstart, &MetOld); MetData = MetOld; SWold = MetOld.ShortWave; - + jday = jstart - 1; + + write_output(jday, SecsPerDay, nsave, stepnum); /************************************************************************** * Loop over all days * **************************************************************************/ @@ -306,6 +309,7 @@ void do_model(int jstart, int nsave) //# If it is the last day, adjust the stop time for the day if necessary if (ntot == nDates) stoptime = stopTOD; if (stoptime == 0) break; + day_fraction = (stoptime - startTOD) / iSecsPerDay; //# Initialise daily values for volume & heat balance reporting (lake.csv) SurfData.dailyRain = 0.; SurfData.dailyEvap = 0.; @@ -326,7 +330,7 @@ void do_model(int jstart, int nsave) //# (stoptime - startTOD) allow for partial dates at the the beginning and end of //# simulation for (i = 0; i < NumInf; i++) { - Inflows[i].FlowRate = (FlowOld[i] + FlowNew[i]) / 2.0 * (stoptime - startTOD) ; + Inflows[i].FlowRate = (FlowOld[i] + FlowNew[i]) / 2.0 * day_fraction * iSecsPerDay; Inflows[i].TemInf = (TempOld[i] + TempNew[i]) / 2.0; Inflows[i].SalInf = (SaltOld[i] + SaltNew[i]) / 2.0; Inflows[i].SubmElev = Elev[i]; @@ -340,7 +344,7 @@ void do_model(int jstart, int nsave) read_daily_outflow(jday, NumOut, DrawNew); //# To get daily outflow (i.e. m3/day) times by the seconds in the current day for (i = 0; i < NumOut; i++) - Outflows[i].Draw = (DrawOld[i] + DrawNew[i]) / 2.0 * (stoptime - startTOD) ; + Outflows[i].Draw = (DrawOld[i] + DrawNew[i]) / 2.0 * day_fraction * iSecsPerDay ; read_daily_withdraw_temp(jday, &WithdrTempNew); WithdrawalTemp = (WithdrTempOld + WithdrTempNew) / 2.0; @@ -378,10 +382,10 @@ void do_model(int jstart, int nsave) SurfData.dailyInflow = do_inflows(); //# Do inflow for all streams //# Extract withdrawal from all offtakes - SurfData.dailyOutflow = do_outflows(jday); + SurfData.dailyOutflow = do_outflows(jday, day_fraction); //# Take care of any overflow - SurfData.dailyOverflow = do_overflow(jday); + SurfData.dailyOverflow = do_overflow(jday, day_fraction); //# Enforce layer limits check_layer_thickness(); @@ -395,14 +399,20 @@ void do_model(int jstart, int nsave) today = jday; #endif write_output(jday, SecsPerDay, nsave, stepnum); + write_diags(jday, calculate_lake_number()); + write_balance(jday); write_step += nsave; - if ( write_step > last_step ) write_step = last_step; #if PLOTS plotstep++; today = -1; #endif } + + if(ntot == nDates & stepnum < nsave){ + fprintf(stderr, " ERROR: NO netcdf output generated because nsave is less total number of time steps in simuluation\n"); + } + /********************************************************************** * End of daily calculations, Prepare for next day and return. * @@ -428,9 +438,6 @@ void do_model(int jstart, int nsave) printf(" Running day %8d, %4.2f%% of days complete%c", jday, ntot*100./nDates, EOLN); fflush(stdout); } - - write_diags(jday, calculate_lake_number()); - write_balance(jday); } //# do while (ntot < nDates) if (quiet < 2) { printf("\n"); fflush(stdout); } /*----------########### End of main daily loop ################-----------*/ @@ -447,6 +454,7 @@ void do_model_non_avg(int jstart, int nsave) AED_REAL SWold, SWnew, DailyKw, DailyEvap; int jday, ntot, stepnum, stoptime; int i, j; + AED_REAL day_fraction; /*************************************************************************** *CAB Note: these WQ arrays should be sized to Num_WQ_Vars not MaxVars, * @@ -465,8 +473,10 @@ void do_model_non_avg(int jstart, int nsave) stepnum = 0; stoptime = iSecsPerDay; SWold = 0.; - + jday = jstart - 1; + + write_output(jday, SecsPerDay, nsave, stepnum); /************************************************************************** * Loop over all days * **************************************************************************/ @@ -477,6 +487,7 @@ void do_model_non_avg(int jstart, int nsave) //# If it is the last day, adjust the stop time for the day if necessary if (ntot == nDates) stoptime = stopTOD; if (stoptime == 0) break; + day_fraction = (stoptime - startTOD) / iSecsPerDay; //# Initialise daily values for volume & heat balance reporting (lake.csv) SurfData.dailyRain = 0.; SurfData.dailyEvap = 0.; @@ -495,7 +506,7 @@ void do_model_non_avg(int jstart, int nsave) //# To get daily inflow (i.e. m3/day) times by SecsPerDay for (i = 0; i < NumInf; i++) { - Inflows[i].FlowRate = FlowNew[i] * (stoptime - startTOD) ; + Inflows[i].FlowRate = FlowNew[i] * day_fraction * iSecsPerDay; Inflows[i].TemInf = TempNew[i]; Inflows[i].SalInf = SaltNew[i]; Inflows[i].SubmElev = Elev[i]; @@ -511,7 +522,7 @@ void do_model_non_avg(int jstart, int nsave) //# (stoptime - startTOD) allow for partial dates at the the beginning and end of //# simulation for (i = 0; i < NumOut; i++) - Outflows[i].Draw = DrawNew[i] * (stoptime - startTOD) ; + Outflows[i].Draw = DrawNew[i] * day_fraction * iSecsPerDay ; read_daily_withdraw_temp(jday, &WithdrTempNew); WithdrawalTemp = WithdrTempNew; @@ -534,11 +545,11 @@ void do_model_non_avg(int jstart, int nsave) SurfData.dailyInflow = do_inflows(); if (Lake[surfLayer].Vol1 > zero) { - //# Extract withdrawal from all offtakes - SurfData.dailyOutflow = do_outflows(jday); + //# Extract withdrawal from all offtakes + SurfData.dailyOutflow = do_outflows(jday, day_fraction); - //# Take care of any overflow - SurfData.dailyOverflow = do_overflow(jday); + //# Take care of any overflow + SurfData.dailyOverflow = do_overflow(jday, day_fraction); } //# Enforce layer limits @@ -552,13 +563,20 @@ void do_model_non_avg(int jstart, int nsave) today = jday; #endif write_output(jday, SecsPerDay, nsave, stepnum); + write_diags(jday, calculate_lake_number()); + write_balance(jday); write_step += nsave; - if ( write_step > last_step ) write_step = last_step; + + //if ( write_step > last_step ) write_step = last_step; #if PLOTS plotstep++; today = -1; #endif } + + if(ntot == nDates & stepnum < nsave){ + fprintf(stderr, " ERROR: NO netcdf output generated because nsave is less total number of time steps in simuluation\n"); + } /********************************************************************** * End of daily calculations, Prepare for next day and return. * @@ -574,9 +592,6 @@ void do_model_non_avg(int jstart, int nsave) printf(" Running day %8d, %4.2f%% of days complete%c", jday, ntot*100./nDates, EOLN); fflush(stdout); } - - write_diags(jday, calculate_lake_number()); - write_balance(jday); } //# do while (ntot < nDates) if (quiet < 2) { printf("\n"); fflush(stdout); } /*----------########### End of main daily loop ################-----------*/ @@ -599,7 +614,8 @@ void do_model_coupled(int step_start, int step_end, AED_REAL WQNew[MaxInf * MaxVars]; int jday, ntot, stepnum, stoptime, cDays; int i, j; - + AED_REAL day_fraction; + /*------------------------------------------------------------------------*/ memset(WQNew, 0, sizeof(AED_REAL)*MaxInf*MaxVars); @@ -611,8 +627,12 @@ void do_model_coupled(int step_start, int step_end, stoptime = iSecsPerDay; SWold = 0.; + + cDays = step_end - step_start + 1; jday = step_start - 1; + + write_output(jday, SecsPerDay, nsave, stepnum); /************************************************************************** * Loop over all days * **************************************************************************/ @@ -623,6 +643,7 @@ void do_model_coupled(int step_start, int step_end, //# If it is the last day, adjust the stop time for the day if necessary if (ntot == nDates) stoptime = stopTOD; if (stoptime == 0) break; + day_fraction = (stoptime - startTOD) / iSecsPerDay; //# Initialise daily values for volume & heat balance reporting (lake.csv) SurfData.dailyRain = 0.; SurfData.dailyEvap = 0.; @@ -640,7 +661,7 @@ void do_model_coupled(int step_start, int step_end, //# (stoptime - startTOD) allow for partial dates at the the beginning and end of //# simulation for (i = 0; i < NumInf; i++) { - Inflows[i].FlowRate = FlowNew[i] * (stoptime - startTOD); + Inflows[i].FlowRate = FlowNew[i] * day_fraction * iSecsPerDay; // Inflows[i].TemInf = TempNew[i]; // Inflows[i].SalInf = SaltNew[i]; for (j = 0; j < Num_WQ_Vars; j++) { @@ -652,7 +673,7 @@ void do_model_coupled(int step_start, int step_end, // read_daily_outflow(jday, NumOut, DrawNew); //# To get daily outflow (i.e. m3/day) times by SecsPerDay for (i = 0; i < NumOut; i++) - Outflows[i].Draw = DrawNew[i] * (stoptime - startTOD); + Outflows[i].Draw = DrawNew[i] * day_fraction * iSecsPerDay; // read_daily_withdraw_temp(jday, &WithdrTempNew); // WithdrawalTemp = WithdrTempNew; @@ -675,11 +696,11 @@ void do_model_coupled(int step_start, int step_end, SurfData.dailyInflow = do_inflows(); if (Lake[surfLayer].Vol1 > zero) { - //# Extract withdrawal from all offtakes - SurfData.dailyOutflow = do_outflows(jday); + //# Extract withdrawal from all offtakes + SurfData.dailyOutflow = do_outflows(jday, day_fraction); - //# Take care of any overflow - SurfData.dailyOverflow = do_overflow(jday); + //# Take care of any overflow + SurfData.dailyOverflow = do_overflow(jday, day_fraction); } //# Enforce layer limits @@ -693,31 +714,39 @@ void do_model_coupled(int step_start, int step_end, today = jday; #endif write_output(jday, SecsPerDay, nsave, stepnum); + write_diags(jday, calculate_lake_number()); + write_balance(jday); write_step += nsave; - if ( write_step > last_step ) write_step = last_step; + + //if ( write_step > last_step ) write_step = last_step; #if PLOTS plotstep++; today = -1; + + #endif + } + + if(ntot == nDates & stepnum < nsave){ + fprintf(stderr, " ERROR: NO netcdf output generated because nsave is less total number of time steps in simuluation\n"); } /********************************************************************** * End of daily calculations, Prepare for next day and return. * **********************************************************************/ SWold = SWnew; - + #ifdef XPLOTS if ( xdisp ) flush_all_plots(); else #endif - if (quiet < 2) { + if (quiet < 2) { printf(" Running day %8d, %4.2f%% of days complete%c", jday, ntot*100./nDates, EOLN); fflush(stdout); } - write_diags(jday, calculate_lake_number()); - write_balance(jday); + } //# do while (ntot < nDates) if (quiet < 2) { printf("\n"); fflush(stdout); } /*----------########### End of main daily loop ################-----------*/ @@ -747,12 +776,14 @@ int do_subdaily_loop(int stepnum, int jday, int stoptime, int nsave, AED_REAL SW **************************************************************************/ iclock = startTOD; last_step = stepnum + ((stoptime - iclock) / noSecs); + if(stepnum == 0){ write_step = stepnum + nsave; - if ( iclock != 0 ) { - int t = mod(((stoptime - iclock) / noSecs), nsave); - if ( t != 0 ) write_step = stepnum + t; } - if ( write_step > last_step ) write_step = last_step; + //if ( iclock != 0 ) { + // int t = mod(((stoptime - iclock) / noSecs), nsave); + // if ( t != 0 ) write_step = stepnum + t; + //} + // printf("last step %d write_step %d\n", last_step, write_step); startTOD = 0; /* from now on start at the beginning of the day */ @@ -822,14 +853,13 @@ int do_subdaily_loop(int stepnum, int jday, int stoptime, int nsave, AED_REAL SW //# then do not output in the subdaily. Output writing is moved to the //# daily loop so it occurs after the inflow and output calculations. - if ( stepnum == write_step && (iclock + noSecs) != SecsPerDay) { + if ( stepnum == write_step && (iclock + noSecs) != SecsPerDay && stepnum != last_step) { #if PLOTS today = jday; #endif write_output(jday, iclock, nsave, stepnum); write_step += nsave; - if ( write_step > last_step ) write_step = last_step; #if PLOTS //if (++n_steps_done > END_STEPS) { int i; for (i = 0; i < NumLayers; i++) show_l_line(2, Lake[i].Height); flush_all_plots(); } plotstep++; @@ -846,6 +876,8 @@ int do_subdaily_loop(int stepnum, int jday, int stoptime, int nsave, AED_REAL SW iclock += noSecs; yearday += part_day_per_step; } //# do while (iclock < iSecsPerDay) + + //if ( write_step > last_step ) write_step = last_step; /************************************************************************** * End of sub-daily loop * **************************************************************************/