///////////////////////////////////////////////////////////////////////////////////////
/// \file framework.cpp
/// \brief Implementation of the framework() function
///
/// \author Ben Smith
/// $Date: 2014-05-14 12:43:55 +0200 (Wed, 14 May 2014) $
///
///////////////////////////////////////////////////////////////////////////////////////

#include "config.h"
#include "framework.h"
#include "commandlinearguments.h"
#include "guessserializer.h"
#include "parallel.h"
#include "inputmodule.h"
#include "driver.h"
#include "canexch.h"
#include "soilwater.h"
#include "somdynam.h"
#include "growth.h"
#include "vegdynam.h"
#include "landcover.h"
#include "bvoc.h"
#include "commonoutput.h"

// ecev3 - not in trunk. Savestate functions copied from RCA branch
#include "savestate.h"
#include "eceframework.h"
#include <string.h>
#include <sstream>
#include <netcdf.h>


using namespace std;
using namespace GuessParallel;

#include <memory>

// Declaration only. Definition below.
//void updatehighcoverage(Stand& stand);
//CLNvoid computeLAIforIFS(Gridcell& gridcell);
void computeIFSVegetationCover(Gridcell& gridcell);
void computeDailyLAIandFPCforIFS(Gridcell& gridcell);


/// Prints the date and time together with the name of this simulation
void print_logfile_heading() {
	xtring datetime;
	unixtime(datetime);

	xtring header = xtring("[LPJ-GUESS  ") + datetime + "]\n\n";
	dprintf((char*)header);

	// Print the title of this run
	std::string dashed_line(50, '-');
	dprintf("\n\n%s\n%s\n%s\n",
		dashed_line.c_str(), (char*)title, dashed_line.c_str());
}

// ecev3
void compute_gridcell_previous_cflux(Gridcell& gridcell) {

	// Compute previous years yearly C fluxes for gridcell
	// before stand fractions change due to landuse change

	if (date.day > 0)
		return;

	gridcell.previousyearsCfluxNat = 0.0;
	gridcell.previousyearsCfluxAnt = 0.0;

	Gridcell::iterator gc_itr = gridcell.begin();
	while (gc_itr != gridcell.end()) {

		// START OF LOOP THROUGH STANDS
		Stand& stand = *gc_itr;

		gridcell.previousyearsCfluxNat += stand.previousyearsCfluxNat * stand.get_gridcell_fraction();
		gridcell.previousyearsCfluxAnt += stand.previousyearsCfluxAnt * stand.get_gridcell_fraction();
	
		++gc_itr;
	}	// End of loop through stands

	gridcell.cLand = gridcell.ccont() + gridcell.previousyearsCfluxNat + gridcell.previousyearsCfluxAnt;
	gridcell.cFlux = 0.0;
}

// ecev3
void compute_gridcell_cflux(Gridcell& gridcell) {

	// Compute daily fluxes released to the atmosphere where previous years
	// yearly fluxes (Dec 31 fluxes) and annual fluxes created on Jan 1
	// have been spread out over this year

	// Gridcell-level C flux from harvest and harvest decomposition associated with landcover change. Updated only on Jan 1
	if (date.day == 0) {
		gridcell.previousyearsCfluxAnt += gridcell.landcover.dcflux_landuse_change + gridcell.landcover.dcflux_harvest_slow;
	}

	Gridcell::iterator gc_itr = gridcell.begin();
	while (gc_itr != gridcell.end()) {

		// START OF LOOP THROUGH STANDS
		Stand& stand = *gc_itr;

		double npp = 0.0;
		double rh = 0.0;

		for (unsigned int p = 0; p < stand.npatch(); p++) {
			Patch& patch = stand[p];

			// Individual C fluxes. All kg C/m2
			npp += patch.fluxes.get_daily_flux(Fluxes::NPP, date.day);
			rh += patch.fluxes.get_daily_flux(Fluxes::SOILC, date.day);

			patch.fluxes.report_flux(Fluxes::CO2NAT, patch.fluxes.get_daily_flux(Fluxes::SOILC, date.day) -
				patch.fluxes.get_daily_flux(Fluxes::NPP, date.day) +
				gridcell.previousyearsCfluxNat / (double)date.year_length());

			patch.fluxes.report_flux(Fluxes::CO2ANTT, gridcell.previousyearsCfluxAnt / (double)date.year_length());
		}

		stand.dcfluxnat_today = (rh / (float)stand.npatch() + gridcell.previousyearsCfluxNat / (double)date.year_length());
		stand.dcfluxant_today = gridcell.previousyearsCfluxAnt / (double)date.year_length();
		stand.dnpp_today = -npp / (float)stand.npatch();

		gridcell.cFlux += ((rh - npp) / (float)stand.npatch() + (gridcell.previousyearsCfluxNat + gridcell.previousyearsCfluxAnt) / (double)date.year_length()) * stand.get_gridcell_fraction();

		// On Dec 31 we must update this stand's previousyearsCflux{Nat,Ant} values to be spread over next year
		if (date.islastday && date.islastmonth) {

			// Reset yearly fluxes
			stand.previousyearsCfluxNat = 0.0;
			stand.previousyearsCfluxAnt = 0.0;

			double estc, firec, reprc, seedc, harvestc, debexcc, orgcleach;

			for (unsigned int p = 0; p < stand.npatch(); p++) {
				Patch& patch = stand[p];

				// Individual C fluxes. All kg C m-2 and all updated on Dec 31
				reprc = patch.fluxes.get_annual_flux(Fluxes::REPRC);
				firec = patch.fluxes.get_annual_flux(Fluxes::FIREC);
				estc = patch.fluxes.get_annual_flux(Fluxes::ESTC);
				seedc = patch.fluxes.get_annual_flux(Fluxes::SEEDC);
				harvestc = patch.fluxes.get_annual_flux(Fluxes::HARVESTC);
				debexcc = patch.fluxes.get_annual_flux(Fluxes::DEBEXCC);
				orgcleach = patch.fluxes.get_annual_flux(Fluxes::LEACHC);

				stand.previousyearsCfluxNat += (estc + firec + reprc + orgcleach - debexcc) / (double)stand.npatch();
				stand.previousyearsCfluxAnt += (seedc + harvestc) / (double)stand.npatch();

				//gridcell.cFlux += patch.soil.aorgCleach / (double)stand.npatch() * stand.get_gridcell_fraction();

				//patch.fluxes.report_flux(Fluxes::CO2NAT, patch.soil.aorgCleach); // Add to get the same result as cFlux_yearly (gridcell.cFlux)
			}
		}
		++gc_itr;
	}	// End of loop through stands
}

// ecev3 NOTE: 
// We keep the trunk versions of simulate_day and framework. This makes it easier to svn MERGE.

/// Simulate one day for a given Gridcell
/**
 * The climate object in the gridcell needs to be set up with
 * the day's forcing data before calling this function.
 *
 * \param gridcell            The gridcell to simulate
 * \param input_module        Used to get land cover fractions
 */

void simulate_day(Gridcell& gridcell, InputModule* input_module) { // ecev3 - TRUNK VERSION

	// Update daily climate drivers etc
	dailyaccounting_gridcell(gridcell, false); // ecev3 - added "false" to compile

	// Calculate daylength, insolation and potential evapotranspiration
	daylengthinsoleet(gridcell.climate);

	if (run_landcover) {
		if (run[CROPLAND]) {
			// Update crop sowing date calculation framework
			crop_sowing_gridcell(gridcell, false); // ecev3 - added false as firstsimulationday
		}
		if (date.day == 0) {
			// Update dynamic management options
			input_module->getmanagement(gridcell);

			// Dynamic landcover and crop fraction data during historical
			// period and create/kill stands.
			landcover_dynamics(gridcell, input_module);

			// Set forest management intensity for all stands this year
			cut_fractions(gridcell);
		}
	}


	Gridcell::iterator gc_itr = gridcell.begin();
	while (gc_itr != gridcell.end()) {

		// START OF LOOP THROUGH STANDS
		Stand& stand = *gc_itr;

		dailyaccounting_stand(stand);

		stand.firstobj();
		while (stand.isobj) {

			// START OF LOOP THROUGH PATCHES

			// Get reference to this patch
			Patch& patch = stand.getobj();
			// Update daily soil drivers including soil temperature
			dailyaccounting_patch(patch,false); // ecev3 - to get this to compile

			// Determine nitrogen fertilisation amount 
			if(run_landcover)
				nfert(patch);

			if (stand.landcover == CROPLAND) {
				// Calculate crop sowing dates
				crop_sowing_patch(patch);
				// Crop phenology
				crop_phenology(patch);
				// necessary updates after changing growingperiod status
				update_patch_fpc(patch);
			}

			// Leaf phenology for PFTs and individuals
			leaf_phenology(patch, gridcell.climate);
			// Interception
			interception(patch, gridcell.climate);
			initial_infiltration(patch, gridcell.climate);
			// Photosynthesis, respiration, evapotranspiration
			canopy_exchange(patch, gridcell.climate);
			// Sum total required irrigation
			irrigation(patch);
			// Soil water accounting, snow pack accounting
			soilwater(patch, gridcell.climate);
			// Daily C allocation (cropland)
			growth_daily(patch);
			// Soil organic matter and litter dynamics
			som_dynamics(patch);

			if (date.islastday && date.islastmonth) {

				// LAST DAY OF YEAR
				// Tissue turnover, allocation to new biomass and reproduction,
				// updated allometry
				growth(stand, patch);
			}
			stand.nextobj();
		}// End of loop through patches

		// Update crop rotation status
		crop_rotation(stand);

		if (date.islastday && date.islastmonth) {
			// LAST DAY OF YEAR
			stand.firstobj();
			while (stand.isobj) {

				// For each patch ...
				Patch& patch = stand.getobj();
				// Establishment, mortality and disturbance by fire
				vegetation_dynamics(stand, patch);
				stand.nextobj();
			}
		}

		++gc_itr;
	}	// End of loop through stands
}


int framework(const CommandLineArguments& args) {

	// The 'mission control' of the model, responsible for maintaining the
	// primary model data structures and containing all explicit loops through
	// space (grid cells/stands) and time (days and years).

	using std::auto_ptr;

	const char* input_module_name = args.get_input_module();

	auto_ptr<InputModule> input_module(InputModuleRegistry::get_instance().create_input_module(input_module_name));

	GuessOutput::OutputModuleContainer output_modules;
	GuessOutput::OutputModuleRegistry::get_instance().create_all_modules(output_modules);

	// Read the instruction file to obtain PFT static parameters and
	// simulation settings
	read_instruction_file(args.get_instruction_file());

	print_logfile_heading();

	// Initialise input/output
	input_module->init();
	output_modules.init();

	// Nitrogen limitation
	if (ifnlim && !ifcentury) {
		fail("\n\nIf nitrogen limitation is switched on then century soil module also needs to be switched on!");
	}

	// bvoc
	if (ifbvoc) {
		initbvoc();
	}

	// Create objects for (de)serializing grid cells
	auto_ptr<GuessSerializer> serializer;
	auto_ptr<GuessDeserializer> deserializer;

	if (save_state) {
		serializer = auto_ptr<GuessSerializer>(new GuessSerializer(state_path, GuessParallel::get_rank_specific(localcomm), GuessParallel::get_num_processes()));
	}

	if (restart) {
		deserializer = auto_ptr<GuessDeserializer>(new GuessDeserializer(state_path));
	}

	while (true) {

		// START OF LOOP THROUGH GRID CELLS

		// Initialise global variable date
		// (argument nyear not used in this implementation)
		date.init(1);

		// Create and initialise a new Gridcell object for each locality
		Gridcell gridcell;

		// Call input module to obtain latitude and driver data for this grid cell.
		if (!input_module->getgridcell(gridcell)) {
			break;
		}

		// Initialise certain climate and soil drivers
		gridcell.climate.initdrivers(gridcell.get_lat());

		if (run_landcover && !restart) {
			// Read landcover and cft fraction data from 
			// data files for the spinup period and create stands
			landcover_init(gridcell, input_module.get());
		}

		if (restart) {
			// Get the whole grid cell from file...
			deserializer->deserialize_gridcell(gridcell);
			// ...and jump to the restart year
			date.year = state_year;
		}

		// Call input/output to obtain climate, insolation and CO2 for this
		// day of the simulation. Function getclimate returns false if last year
		// has already been simulated for this grid cell

		while (input_module->getclimate(gridcell)) {

			// START OF LOOP THROUGH SIMULATION DAYS

			compute_gridcell_previous_cflux(gridcell);

			simulate_day(gridcell, input_module.get());

			// Compute daily fluxes released to the atmosphere where previous years
			// yearly fluxes (Dec 31 fluxes) and annual fluxes created on Jan 1
			// have been spread out over this year
			compute_gridcell_cflux(gridcell);

			output_modules.outdaily(gridcell);

			if (date.islastday && date.islastmonth) {
				// LAST DAY OF YEAR
				// Call output module to output results for end of year
				// or end of simulation for this grid cell
				output_modules.outannual(gridcell);

				gridcell.balance.check_year(gridcell);

				// Time to save state?
				if (date.year == state_year-1 && save_state) {
					serializer->serialize_gridcell(gridcell);
				}

				// Check whether to abort
				if (abort_request_received()) {
					return 99;
				}
			}

			// Advance timer to next simulation day
			date.next();

			// End of loop through simulation days
		}	//while (getclimate())

		gridcell.balance.check_period(gridcell);

	}		// End of loop through grid cells

	// END OF SIMULATION

	return 0;
}



///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
// EC-Earth Coupling Code - THE GUESS WORLD
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////

// The one and only world (vector of gridcells, one forest and one vegetated open land stand
// for each grid cell)
typedef std::vector<Gridcell*> World;
World world;

// Initialisation flag (whether GUESS has been initialised yet true/false)
bool ifinit=false;

int inst=0;

/// Number of years in the IFS output, spinup dataset
const int NYEARIFS = 10;

/// Number of years in the N deposition output, historical dataset (1850-01 - 2014-12)
const int NYEARNDEP = 165;

// ECE Parameters to be read from .ins file:

/// atmospheric nitrogen deposition (kgN/yr/ha)
double ndep;

// Determine RCP from guess.ins
double rcpforcing;

// Use IFS/H-TESSEl soil temp & moisture (true) or not (false)
//CLNbool ifssoilproperties;
bool useifssoilmoist;
bool useifssoiltemp;

xtring statefile;

// ecev3 - pointers
using std::auto_ptr;

// Objects for serializing and deserializing Gridcells
static std::auto_ptr<GuessSerializer> serializer_ece;
static std::auto_ptr<GuessDeserializer> deserializer_ece;

// For output
GuessOutput::OutputModuleContainer* output_modules_ece;

// For input
// Will give access to LULCC, crop types and management options for each cell. 
InputModule* input_module_ece;


void initlog(int inst) {

	printf("Attempting to open log file for instance number %d\n",inst);

	const char* file_log_prefix = "guess";

	xtring filename;
	filename.printf("%s%d.log", file_log_prefix, inst);

	set_shell(new CommandLineShell((char*)filename));
}



// EC-Earth VERSION of simulate_day. 

/// Simulate one day for a given Gridcell
/**
 * The climate object in the gridcell needs to be set up with
 * the day's forcing data before calling this function.
 *
 * \param gridcell            The gridcell to simulate
 * \param prescribe_ifs_soilmoist  New ecev3 option. Use IFS soil moisture if true, 
 *                            but standard LPJ-GUESS parameterisations if false 
 * \param prescribe_ifs_soiltemp  New ecev3 option. Use IFS soil temperature if true, 
 *                            but standard LPJ-GUESS parameterisations if false 
 * \param firstsimulationday  New ecev3 option sent to dailyaccounting_gridcell 
 *                             
 * \param printoutput         New ecev3 option. Update outputfiles if true. Set false at start of spinup 
 *                            
 */

void simulate_day_ece(Gridcell& gridcell, bool prescribe_ifs_soilmoist, bool prescribe_ifs_soiltemp, 
		      bool firstsimulationday, bool printoutput) {

	// Compute previous years yearly C fluxes for gridcell
	// before stand fractions change due to landuse change
	compute_gridcell_previous_cflux(gridcell);

	// Update daily climate drivers etc
	dailyaccounting_gridcell(gridcell, firstsimulationday);

	// Calculate daylength, insolation and potential evapotranspiration
	daylengthinsoleet(gridcell.climate);

	if (run_landcover) {
		if (run[CROPLAND]) {
			// Update crop sowing date calculation framework
			crop_sowing_gridcell(gridcell, firstsimulationday);
		}
		if (date.day == 0) {
			// Dynamic landcover and crop fraction data during historical
			// period and create/kill stands.
			landcover_dynamics(gridcell, input_module_ece);

			// Update dynamic management options
			input_module_ece->getmanagement(gridcell);
		}
	}


	Gridcell::iterator gc_itr = gridcell.begin();
	while (gc_itr != gridcell.end()) {

		// START OF LOOP THROUGH STANDS
		Stand& stand = *gc_itr;

		dailyaccounting_stand(stand);

		stand.firstobj();
		while (stand.isobj) {

			// START OF LOOP THROUGH PATCHES

			// Get reference to this patch
			Patch& patch = stand.getobj();
			// Update daily soil drivers including soil temperature
			dailyaccounting_patch(patch,prescribe_ifs_soiltemp); // ecev3 - override soil.temp?
															 // Determine nitrogen fertilisation amount 
			if (run_landcover)
				nfert(patch);

			if (stand.landcover == CROPLAND) {
				// Calculate crop sowing dates
				crop_sowing_patch(patch);
				// Crop phenology
				crop_phenology(patch);
				// necessary updates after changing growingperiod status
				update_patch_fpc(patch);
			}

			// Leaf phenology for PFTs and individuals
			leaf_phenology(patch, gridcell.climate);
			// Interception
			interception(patch, gridcell.climate);
			initial_infiltration(patch, gridcell.climate);
			// Photosynthesis, respiration, evapotranspiration
			canopy_exchange(patch, gridcell.climate);

			// Sum total required irrigation
			irrigation(patch);

			// Soil water accounting, snow pack accounting
			if (prescribe_ifs_soilmoist) {

				// Update wcont[0], wcont[0] and wcont_evap with IFS values
				patch.soil.wcont[0] = gridcell.climate.ifsw_surf;
				patch.soil.wcont[1] = gridcell.climate.ifsw_deep;
				patch.soil.wcont_evap = gridcell.climate.ifsw_surf; 

				// MUST call soil water, otherwise aaet and dperc are not updated.
				// These are needed for N limitation etc.
				soilwater(patch, gridcell.climate);

				// Update daily soil water in  both layers
				// Done in dailyaccounting_patch, but here we override with IFS values before 
				// these values are used in the fire module.
				patch.soil.wcont[0] = gridcell.climate.ifsw_surf;
				patch.soil.wcont[1] = gridcell.climate.ifsw_deep;
				patch.soil.dwcontupper[date.day] = patch.soil.wcont[0];
				patch.soil.dwcontlower[date.day] = patch.soil.wcont[1];
				patch.soil.dwcontlower_today = patch.soil.wcont[1];

			} 
			else {
				soilwater(patch, gridcell.climate);
			}

			// Daily C allocation (cropland)
			growth_daily(patch);
			// Soil organic matter and litter dynamics
			som_dynamics(patch);

			if (date.islastday && date.islastmonth) {

				// LAST DAY OF YEAR
				// Tissue turnover, allocation to new biomass and reproduction,
				// updated allometry
				growth(stand, patch);
			}
			stand.nextobj();
		} // End of loop through patches

		// Update crop rotation status
		crop_rotation(stand);

		if (date.islastday && date.islastmonth) {
			// LAST DAY OF YEAR
			stand.firstobj();
			while (stand.isobj) {

				// For each patch ...
				Patch& patch = stand.getobj();
				// Establishment, mortality and disturbance by fire
				vegetation_dynamics(stand, patch);
				stand.nextobj();
			}

		}

		// gridcell.nextobj();
		++gc_itr;
	}	// End of loop through stands

        // TODO check if outdaily should be after compute_gridcell_cflux() as in the standalone framework

	// Compute daily fluxes released to the atmosphere where previous years
	// yearly fluxes (Dec 31 fluxes) and annual fluxes created on Jan 1
	// have been spread out over this year
	compute_gridcell_cflux(gridcell);

	// ecev3 - not needed yet
	if (printoutput)
		(*output_modules_ece).outdaily(gridcell);

	if (date.islastday && date.islastmonth) {
		// LAST DAY OF YEAR

		// ecev3 - update the simulation year for this gridcell
		gridcell.simulationyear++;

		// Call output module to output results for end of year
		// or end of simulation for this grid cell
		if (printoutput) {

			(*output_modules_ece).outannual(gridcell);

			// ecev3 - compute vegetation type and cover fractions for this grid cell
			// outannual must be called first!
			computeIFSVegetationCover(gridcell);
		} 

		// for debugging:
		// dprintf("Stands: %4d... \n", gridcell.nbr_stands());

	} // last day?

}



///////////////////////////////////////////////////////////////////////////////////////
// INITGUESS
// Called at start of run to read settings from ins file
// For now ins file name is prescribed in the function body

void initguess(int myrank) {
	
	// Hard-wired ins file name
	xtring insfilename="guess.ins";
	
	// Get instance id
	inst = myrank;

	// Open log file
	initlog(inst);

	// Retrieve specified N value as read from ins file
	ndep=param["ndep"].num;

	// Determine RCP to follow
	rcpforcing=param["nrcp"].num;

	// Use IFS soil moisture and temp (1), or not (0) 
	useifssoilmoist=true;
	if (param["useifssoilmoist"].num == 0.0)
		useifssoilmoist=false;

	useifssoiltemp=true;
	if (param["useifssoiltemp"].num == 0.0)
		useifssoiltemp=false;

	// State file params
	restart=true;
	if (param["restart"].num == 0.0)
		restart=false;

	state_path=param["state_path"].str;
	state_name=param["state_name"].str;

	// Nitrogen limitation
	if (ifnlim && !ifcentury) {
		fail("\n\nIf nitrogen limitation is switched on then century soil module also needs to be switched on!");
	}

	// ecev3 - CMIP6
	if (!printcmip6 && printcmip6daily) {
		fail("\n\nOnly print daily CMIP6 output when printcmip6 is true too!");
	}

	// ecev3 - CRESCENDO
	if (!printcrescendo && printcrescendodaily) {
		fail("\n\nOnly print daily CRESCENDO output when printcrescendo is true too!");
	}

	// ecev3 - CRESCENDO FACE
	bool printCRESCENDO_FACE = false;
#ifdef CRESCENDO_FACE
	printCRESCENDO_FACE = true;
#endif
	if (printcrescendofacedaily && !printCRESCENDO_FACE) {
		fail("\n\nOnly print daily CRESCENDO FACE output when CRESCENDO_FACE is defined!");
	}
	if (!printcrescendofacedaily && printCRESCENDO_FACE) {
		fail("\n\nOnly define CRESCENDO_FACE if daily CRESCENDO FACE outputs are to be printed!");
	}

	// bvoc
	if (ifbvoc) {
	  initbvoc();
	}

	// Done it!
	ifinit=true;
}


// ecev3
void swConvert(Soiltype& stype,double &sw_surf,double &sw_deep) {

	if (sw_surf < stype.wilt_point)
		sw_surf = 0.0;
	else if (sw_surf > stype.field_cap)
		sw_surf = 1.0;
	else 
		sw_surf = (sw_surf - stype.wilt_point) / stype.whcap;

	if (sw_deep < stype.wilt_point)
		sw_deep = 0.0;
	else if (sw_deep > stype.field_cap)
		sw_deep = 1.0;
	else
		sw_deep = (sw_deep - stype.wilt_point) / stype.whcap;

}


// ecev3
/* 
Called daily to update gridcell.laiphen_high_today and gridcell.laiphen_high_today

computeIFSVegetationCover below updates these Gridcell values on Dec 31:

	gridcell.IFStypehigh
	gridcell.IFSfrachigh
	gridcell.IFStypelow
	gridcell.IFSfraclow
	gridcell.transferhightolow - use total LAI today if this is true

	These values are serialized, so they are available after a restart 
*/

void computeDailyLAIandFPCforIFS(Gridcell& gridcell) {

	const double MAX_DAILY_LAI = 10; // Max daily LAI - enforced at the end of this routine. Was 15.
	double ABSCUTOFF_FPC = 0.000; // Set to 0 to ensure we send nonzero values no matter how small. Was 0.005.

	double lai_tree_daily = 0.0;
	double lai_grass_daily = 0.0;
	double lai_crop_daily = 0.0;
	double lai_pasture_daily = 0.0;
	double lai_urban_daily = 0.0;
	double lai_peat_daily = 0.0;

	double fpc_tree_daily = 0.0;
	double fpc_grass_daily = 0.0;
	double fpc_crop_daily = 0.0;
	double fpc_pasture_daily = 0.0;
	double fpc_urban_daily = 0.0;
	double fpc_peat_daily = 0.0;

	// As updated in outannual
	double natural_frac = gridcell.landcover.frac[NATURAL];
	double cropland_frac = gridcell.landcover.frac[CROPLAND];
	double pasture_frac = gridcell.landcover.frac[PASTURE];
	double urban_frac = gridcell.landcover.frac[URBAN];
	double peatland_frac = gridcell.landcover.frac[PEATLAND];

	// Sum LAI and FPC across patches and PFTs
	Gridcell::iterator gc_itr = gridcell.begin();

	// Loop through Stands
	gc_itr = gridcell.begin();

	while (gc_itr != gridcell.end()) {
		Stand& stand = *gc_itr;
		stand.firstobj();

		double to_gridcell_average = stand.get_gridcell_fraction() / (double)stand.npatch();

		if (stand.landcover == NATURAL) { // Only look in NATURAL stands!

			// Loop through Patches for trees first 
			while (stand.isobj) {
				Patch& patch = stand.getobj();

				// Loop through Individuals
				Vegetation& vegetation = patch.vegetation;

				// Tree and grass fpc for this patch
				double fpc_tree_patch = 0.0;
				double fpc_grass_patch = 0.0;

				vegetation.firstobj();
				while (vegetation.isobj) {
					Individual& indiv = vegetation.getobj();

					if (indiv.id != -1 && indiv.alive) {

						double current_fpc = indiv.crownarea * indiv.densindiv * (1.0 - lambertbeer(indiv.lai_indiv_today()));

						if (indiv.pft.lifeform == TREE) {

							lai_tree_daily += indiv.lai_today() * to_gridcell_average;
							fpc_tree_patch += current_fpc * to_gridcell_average;
						}
						else {
							lai_grass_daily += indiv.lai_today() * to_gridcell_average;
							fpc_grass_patch += current_fpc * to_gridcell_average;
						}
					} // alive && id != -1?
					vegetation.nextobj();
				} // while/vegetation loop for trees

				// Tree FPC shouldn't be larger than total patch area
				if (fpc_tree_patch > to_gridcell_average)
					fpc_tree_patch = to_gridcell_average;

				// Tree and grass FPC shouldn't be larger than total patch area (trees shade grass)
				if (fpc_tree_patch + fpc_grass_patch > to_gridcell_average)
					fpc_grass_patch = to_gridcell_average - fpc_tree_patch;

				fpc_tree_daily += fpc_tree_patch;
				fpc_grass_daily += fpc_grass_patch;

				stand.nextobj();
			} // patch loop 
		}
		else if (stand.landcover == PASTURE || stand.landcover == URBAN || stand.landcover == PEATLAND) {

			// Look in pasture, peatland and urban stands for C3/C4

			// Loop through Patches
			while (stand.isobj) {
				Patch& patch = stand.getobj();
				
				// Loop through Individuals
				Vegetation& vegetation = patch.vegetation;

				vegetation.firstobj();
				while (vegetation.isobj) {
					Individual& indiv = vegetation.getobj();

					if (indiv.id != -1 && indiv.alive) {

						if (indiv.pft.lifeform == GRASS) {

							double current_fpc_daily = 1.0 - lambertbeer(indiv.lai_indiv_today());

							if (stand.landcover == PASTURE) {
								lai_pasture_daily += indiv.lai_today() * to_gridcell_average;
								fpc_pasture_daily += current_fpc_daily * to_gridcell_average;
							}
							else if (stand.landcover == URBAN) {
								lai_urban_daily += indiv.lai_today() * to_gridcell_average;
								fpc_urban_daily += current_fpc_daily * to_gridcell_average;
							}
							else {
								lai_peat_daily += indiv.lai_today() * to_gridcell_average;
								fpc_peat_daily += current_fpc_daily * to_gridcell_average;
							}
						} // GRASS?
						else {
							// Error!
							bool err = true;
						}
					} // alive && id != -1?
					vegetation.nextobj();
				} // while/vegetation loop
				stand.nextobj();
			} // patch loop
		}
		else if (stand.landcover == CROPLAND) {
			
			// Look in cropland to find irrigated and none irrigated / C3 or C4 crops

			// Loop through Patches
			while (stand.isobj) {
				Patch& patch = stand.getobj();

				// Loop through Individuals
				Vegetation& vegetation = patch.vegetation;
		
				vegetation.firstobj();
				while (vegetation.isobj) {
					Individual& indiv = vegetation.getobj();

					// Calculating real FCP as indiv.fpc for crops are always 1.0!
					double current_fpc_daily = 1.0 - lambertbeer(indiv.lai_indiv_today());

					lai_crop_daily += indiv.lai_today() * to_gridcell_average;
					fpc_crop_daily += current_fpc_daily * to_gridcell_average;

					vegetation.nextobj();
				} // while/vegetation loop

				stand.nextobj();
			} // patch loop
		}
		++gc_itr;
	} // stand loop

	// Rescale if FPC > Land cover fraction
	// If Natural tree bigger than cover fracton -> tree cover all fraction, grass zero
	if (fpc_tree_daily > natural_frac) {
		fpc_tree_daily = natural_frac;
		fpc_grass_daily = 0.0;
	}

	double fpc_natural_daily = fpc_grass_daily + fpc_tree_daily;
	if (fpc_natural_daily > natural_frac) {
		fpc_grass_daily = natural_frac - fpc_tree_daily; // Trees covers low vegetation
	}

	// Pasture
	if (fpc_pasture_daily > pasture_frac) {
		fpc_pasture_daily = pasture_frac;
	}

	// Urban
	if (fpc_urban_daily > urban_frac) {
		fpc_urban_daily = urban_frac;
	}

	// Peatland
	if (fpc_peat_daily > peatland_frac) {
		fpc_peat_daily = peatland_frac;
	}

	// Cropland
	if (fpc_crop_daily > cropland_frac) {
		fpc_crop_daily = cropland_frac;
	}

	double fpc_grass_total_daily = fpc_grass_daily + fpc_crop_daily + fpc_pasture_daily + fpc_urban_daily + fpc_peat_daily;

	// Lower limit - Low
	if (fpc_grass_total_daily < ABSCUTOFF_FPC) {
		fpc_grass_total_daily = 0.0; // Impose lower limit of 0.5%, 0 otherwise
	}

	// Lower limit - High
	if (fpc_tree_daily < ABSCUTOFF_FPC) {
		fpc_tree_daily = 0.0; // Impose lower limit of 0.5%, 0 otherwise
	}

	// Catch small rounding errors
	double coversum = fpc_tree_daily + fpc_grass_total_daily;
	if (coversum > 1.0 && coversum <= 1.001) {
		fpc_tree_daily /= coversum;
		fpc_grass_total_daily /= coversum;
		coversum = 1.0;
	}

	// Sanity check!!!
	if (coversum > 1.001)
		fail("\nInvalid gridcell.IFSfrachigh (%f) and/or gridcell.IFSfraclow (%f) value (sum = %f) in computeDailyLAIandFPCforIFS!",
			fpc_tree_daily, fpc_grass_total_daily, fpc_tree_daily + fpc_grass_total_daily);


	// Determine FPC based LAI
	double lai_fpc_tree_daily = 0.0;
	if (fpc_tree_daily > 0.0) {
		lai_fpc_tree_daily = lai_tree_daily / fpc_tree_daily;
	}

	double lai_fpc_grass_daily = 0.0;
	if (fpc_grass_total_daily > 0.0) {
		lai_fpc_grass_daily = (lai_grass_daily + lai_crop_daily + lai_pasture_daily + lai_urban_daily + lai_peat_daily) / fpc_grass_total_daily;
	}

	// Sanity checks - limit LAI to sensible values
	if (lai_fpc_tree_daily >= MAX_DAILY_LAI)
		lai_fpc_tree_daily = MAX_DAILY_LAI;

	if (lai_fpc_grass_daily >= MAX_DAILY_LAI)
		lai_fpc_grass_daily = MAX_DAILY_LAI;

	gridcell.laiphen_high_today = lai_fpc_tree_daily;
	gridcell.laiphen_low_today = lai_fpc_grass_daily;
	gridcell.IFSfrachigh = fpc_tree_daily;
	gridcell.IFSfraclow = fpc_grass_total_daily;
}


// ecev3 - called on Dec 31
void computeIFSVegetationCover(Gridcell& gridcell) {

	/*
	// LPJG vegetation mapped onto H-TESSEL vegetation types and fractions.
	// Types are one of (See IFS documentation, 36r1, Table 8.1):

	0 No high vegetation
	1 Crops, mixd farming
	2 Short grass
	3 Evergreen needleleaf trees
	4 Deciduous needleleaf trees
	5 Deciduous broadleaf trees
	6 Evergreen broadleaf trees
	7 Tall grass
	8 Desert
	9 Tundra
	10 Irrigated crops
	11 Semidesert
	12 Ice caps and glaciers
	13 Bogs and marshes
	14 Inland water
	15 Ocean
	16 Evergreen shrubs
	17 Deciduous shrubs
	18 Mixed forest/woodland
	19 Interrupted forest
	20 Water and land mixtures
	*/

	// -------------------------------------------------------------------------
	// 0. Initial constants and checks
	// -------------------------------------------------------------------------

	const double FORESTCUTOFF_FPC = 0.6;			// FPC of 60%
	const double DNFORESTCUTOFF_FPC = 0.4;			// FPC of 40%
	const double LOWTOHIGHCUTOFF_FPC = 0.05;		// FPC of 5%
	const double DOMINANCE_FRACTION = 0.6;			// 60% (IGBP)
	const double DOMINANCE_FRACTION_RAINGREEN = 0.3;// 30% (Otherwise raingreen dominance is not captured)
	const double TUNDRAGDD5LIMIT = 350.0;			// degree day
	const double DESERTGDD5LIMIT = 1500.0;			// degree day
	const double ABSCUTOFF_FPC = 0.000;				// MIN FPC (as in IFS files now)
	// const double DESERTCUTOFF_FPC = 0.05;			// FPC of 5%
	// const double SEMIDESERTCUTOFF_FPC = 0.4;		// FPC of 40%	
	const double DESERTCUTOFF_FPC_WCONT = 0.01;		// FPC * AWCONT of 0.01
	const double SEMIDESERTCUTOFF_FPC_WCONT = 0.15;	// FPC * AWCONT of 0.15	
	const double MAX_DAILY_LAI = 10;				// Max daily LAI - enforced at the end of this routine. Was 15.

	// Initialise these Gridcell values before updating them below
	gridcell.IFStypehigh = 0;
	double IFSfrachigh = 0.0;
	gridcell.IFStypelow = 0;
	double IFSfraclow = 0.0;
	gridcell.naturaltreeFPC = 0.0;
	gridcell.naturalgrassFPC = 0.0;

	// Whether we transfer possible tree characteristics from the NATURAL stand to the low veg information sent to IFS
	gridcell.transferhightolow = false;

	// Some quick checks
	if (gridcell.get_lat() < -65.0) {
		// Antarctica
		// gridcell.IFStypelow=12; // 12 Ice caps and glaciers
		return;
	}


	// As updated in outannual
	double natural_frac = gridcell.landcover.frac[NATURAL];
	double cropland_frac = gridcell.landcover.frac[CROPLAND];
	double pasture_frac = gridcell.landcover.frac[PASTURE];
	double urban_frac = gridcell.landcover.frac[URBAN];
	double peatland_frac = gridcell.landcover.frac[PEATLAND];

	// -------------------------------------------------------------------------
	// 1. Determine FPC for trees and grasses in Natural stands and real FPC for crops.
	// -------------------------------------------------------------------------

	double fpcEN = 0.0;
	double fpcDN = 0.0;
	double fpcDB = 0.0;
	double fpcRB = 0.0;
	double fpcEB = 0.0;
	double fpcgrass = 0.0;

	// To distinguish tall from short grass we need C3/C4 distinctions
	double fpcgrass_c3 = 0.0;
	double fpcgrass_c4 = 0.0;
	double fpcgrass_crop_c3 = 0.0;
	double fpcgrass_crop_c4 = 0.0;
	double fpcgrass_pasture_c3 = 0.0;
	double fpcgrass_pasture_c4 = 0.0;
	double fpcgrass_urban_c3 = 0.0;
	double fpcgrass_urban_c4 = 0.0;
	double fpcgrass_peatland_c3 = 0.0;
	double fpcgrass_peatland_c4 = 0.0;

	// To distinguish rainfed and irrigated crops
	double fpccrop_rain = 0.0;
	double fpccrop_irri = 0.0;
	double fpccrop = 0.0;

	Gridcell::iterator gc_itr = gridcell.begin();

	// Loop through Stands
	while (gc_itr != gridcell.end()) {
		Stand& stand = *gc_itr;
		stand.firstobj();

		double to_gridcell_average = stand.get_gridcell_fraction() / (double)stand.npatch();

		if (stand.landcover == NATURAL) { // Only look in NATURAL stands!

			// Loop through Patches
			while (stand.isobj) {
				Patch& patch = stand.getobj();

				// Loop through Individuals
				Vegetation& vegetation = patch.vegetation;

				vegetation.firstobj();
				while (vegetation.isobj) {
					Individual& indiv = vegetation.getobj();

					if (indiv.id != -1 && indiv.alive) {

						if (indiv.pft.lifeform == GRASS) {
							fpcgrass += indiv.fpc * to_gridcell_average;

							// C3 or C4?
							if (indiv.pft.pathway == C3) {
								fpcgrass_c3 += indiv.fpc * to_gridcell_average;
							}
							else {
								fpcgrass_c4 += indiv.fpc * to_gridcell_average;
							}
						}
						else if (indiv.pft.leafphysiognomy == NEEDLELEAF && indiv.pft.phenology == EVERGREEN) {
							fpcEN += indiv.fpc * to_gridcell_average;
						}
						else if (indiv.pft.leafphysiognomy == NEEDLELEAF && indiv.pft.phenology == SUMMERGREEN) {
							fpcDN += indiv.fpc * to_gridcell_average;
						}
						else if (indiv.pft.leafphysiognomy == BROADLEAF && indiv.pft.phenology == SUMMERGREEN) {
							fpcDB += indiv.fpc * to_gridcell_average;
						}
						else if (indiv.pft.leafphysiognomy == BROADLEAF && indiv.pft.phenology == EVERGREEN) {
							fpcEB += indiv.fpc * to_gridcell_average;
						}
						else if (indiv.pft.leafphysiognomy == BROADLEAF && indiv.pft.phenology == RAINGREEN) {
							fpcDB += indiv.fpc * to_gridcell_average;
							fpcRB += indiv.fpc * to_gridcell_average;
						}
						else {
							// Error!
							bool err = true;
						}
					} // alive && id != -1?
					vegetation.nextobj();
				} // while/vegetation loop
				stand.nextobj();
			} // patch loop
		}
		else if (stand.landcover == PASTURE || stand.landcover == URBAN || stand.landcover == PEATLAND) {

			// Look in pasture, peatland and urban stands for C3/C4

			// Loop through Patches
			while (stand.isobj) {
				Patch& patch = stand.getobj();

				// Loop through Individuals
				Vegetation& vegetation = patch.vegetation;

				vegetation.firstobj();
				while (vegetation.isobj) {
					Individual& indiv = vegetation.getobj();

					if (indiv.id != -1 && indiv.alive) {

						if (indiv.pft.lifeform == GRASS) {

							// C3 or C4?
							if (indiv.pft.pathway == C3) {

								// C3
								if (stand.landcover == PASTURE)
									fpcgrass_pasture_c3 += indiv.fpc * to_gridcell_average;
								else if (stand.landcover == URBAN)
									fpcgrass_urban_c3 += indiv.fpc * to_gridcell_average;
								else
									fpcgrass_peatland_c3 += indiv.fpc * to_gridcell_average;
							}
							else {

								// C4
								if (stand.landcover == PASTURE)
									fpcgrass_pasture_c4 += indiv.fpc * to_gridcell_average;
								else if (stand.landcover == URBAN)
									fpcgrass_urban_c4 += indiv.fpc * to_gridcell_average;
								else
									fpcgrass_peatland_c4 += indiv.fpc * to_gridcell_average;
							}
						} // GRASS?
						else {
							// Error!
							bool err = true;
						}
					} // alive && id != -1?
					vegetation.nextobj();
				} // while/vegetation loop
				stand.nextobj();
			} // patch loop
		}
		else if (stand.landcover == CROPLAND) {

			// Look in cropland to find irrigated and none irrigated / C3 or C4 crops

			// Loop through Patches
			while (stand.isobj) {
				Patch& patch = stand.getobj();

				// Loop through Individuals
				Vegetation& vegetation = patch.vegetation;

				vegetation.firstobj();
				while (vegetation.isobj) {
					Individual& indiv = vegetation.getobj();

					if (!indiv.pft.isintercropgrass) {

						// Calculating real FCP as indiv.fpc for crops are always 1.0!
						double indiv_fpc = 1.0 - lambertbeer(indiv.lai);

						if (patch.stand.isirrigated)
							fpccrop_irri += indiv_fpc * to_gridcell_average;
						else
							fpccrop_rain += indiv_fpc * to_gridcell_average;

						if (indiv.pft.pathway == C3)
							fpcgrass_crop_c3 += indiv_fpc * to_gridcell_average;
						else
							fpcgrass_crop_c4 += indiv_fpc * to_gridcell_average;

						fpccrop += indiv_fpc * to_gridcell_average;
					}
					vegetation.nextobj();
				} // while/vegetation loop

				stand.nextobj();
			} // patch loop
		}
		++gc_itr;

	} // stand loop

	  // Now determine the maximum type
	double fpcnatural = fpcEN + fpcDN + fpcDB + fpcEB + fpcgrass;
	double fpctree = fpcEN + fpcDN + fpcDB + fpcEB;
	double fpcpeatland = fpcgrass_peatland_c3 + fpcgrass_peatland_c4;
	double fpcpasture = fpcgrass_pasture_c3 + fpcgrass_pasture_c4;
	double fpcurban = fpcgrass_urban_c3 + fpcgrass_urban_c4;

	// Rescale if FPC > Land cover fraction
	// Natural
	if (fpcnatural > natural_frac) {
		fpcEN *= (natural_frac / fpcnatural);
		fpcDN *= (natural_frac / fpcnatural);
		fpcDB *= (natural_frac / fpcnatural);
		fpcRB *= (natural_frac / fpcnatural);
		fpcEB *= (natural_frac / fpcnatural);

		fpctree *= (natural_frac / fpcnatural);
		fpcgrass *= (natural_frac / fpcnatural);
		fpcgrass_c3 *= (natural_frac / fpcnatural);
		fpcgrass_c4 *= (natural_frac / fpcnatural);

		// Reset total to
		fpcnatural = natural_frac;
	}

	// Pasture
	if (fpcpasture > pasture_frac) {
		fpcgrass_pasture_c3 *= (pasture_frac / fpcpasture);
		fpcgrass_pasture_c4 *= (pasture_frac / fpcpasture);

		// Reset total to
		fpcpasture = pasture_frac;
	}

	// Urban
	if (fpcurban > urban_frac) {
		fpcgrass_urban_c3 *= (urban_frac / fpcurban);
		fpcgrass_urban_c4 *= (urban_frac / fpcurban);

		// Reset total to
		fpcurban = urban_frac;
	}

	// Peatland
	if (fpcpeatland > peatland_frac) {
		fpcgrass_peatland_c3 *= (peatland_frac / fpcpeatland);
		fpcgrass_peatland_c4 *= (peatland_frac / fpcpeatland);

		// Reset total to
		fpcpeatland = peatland_frac;
	}

	// Cropland
	if (fpccrop > cropland_frac) {
		fpccrop_rain *= (cropland_frac / fpccrop);
		fpccrop_irri *= (cropland_frac / fpccrop);
		fpcgrass_crop_c3 *= (cropland_frac / fpccrop);
		fpcgrass_crop_c4 *= (cropland_frac / fpccrop);

		// Reset total to
		fpccrop = cropland_frac;
 	}

	// -------------------------------------------------------------------------
	// 1. Lump together to totals
	// -------------------------------------------------------------------------

	// Note: we limit FPCs to be <= 1
	double fpcgridcell = fpccrop + fpcpasture + fpcnatural + fpcurban + fpcpeatland;

	// Gridcell C3/C4 grass distinctions - NATURAL, PASTURE, URBAN & PEATLAND - must be between 0 and 1
	double fpcgrasses_c3 = fpcgrass_c3 + fpcgrass_crop_c3 + fpcgrass_pasture_c3 + fpcgrass_urban_c3 + fpcgrass_peatland_c3;
	double fpcgrasses_c4 = fpcgrass_c4 + fpcgrass_crop_c4 + fpcgrass_pasture_c4 + fpcgrass_urban_c4 + fpcgrass_peatland_c4;
	double fpcgrasses = fpcgrasses_c3 + fpcgrasses_c4;

	// -------------------------------------------------------------------------
	// 2. Forest types
	// -------------------------------------------------------------------------

	int typehigh;

	// *** FOREST ***

	// 4 Deciduous needleleaf trees
	// 3 Evergreen needleleaf trees 
	// 5 Deciduous broadleaf trees
	// 6 Evergreen broadleaf trees
	// or
	// 18 Mixed forest/woodland
	
	if (fpcDN > DOMINANCE_FRACTION * fpctree)
		typehigh = 4;
	else if (fpcEN > DOMINANCE_FRACTION * fpctree)
		typehigh = 3;
	else if (fpcDB > DOMINANCE_FRACTION * fpctree || fpcRB > DOMINANCE_FRACTION_RAINGREEN * fpctree)
		typehigh = 5;
	else if (fpcEB > DOMINANCE_FRACTION * fpctree)
		typehigh = 6;
	else if (!negligible(fpctree))
		typehigh = 18; // Mixed if no type exceeds 60% of the total
	else
		typehigh = 0;

	gridcell.IFStypehigh = typehigh;

	// -------------------------------------------------------------------------
	// 3. Low types and fractions
	// -------------------------------------------------------------------------

	//  Quick, initial check for desert and glacier-type conditions
	if (negligible(fpcgrasses)) {

		gridcell.IFStypelow = 0; // No Vegetation
	}
	// Peatland check - and this must depend on the hard-coded wetland fraction
	else if (fpcpeatland > (fpcgrasses - fpcpeatland))
		gridcell.IFStypelow = 13; // Bogs and marshes
	else {

		// Not Bogs and marshes ....

		// Crops, mixed farming dominate
		if (fpccrop > fpcgrass && fpccrop > fpcpasture) {

			if (fpccrop_irri > fpccrop_rain)
				gridcell.IFStypelow = 10; // Irrigated crops 10
			else
				gridcell.IFStypelow = 1; // Crops, mixed farming
		}
		else {

			// Natural GRASSES (and maybe trees) dominate the low cover types
			if (gridcell.climate.agdd5_5.mean() < TUNDRAGDD5LIMIT) {
				gridcell.IFStypelow = 9; // Tundra (Hickler et al. 2006)
			}
			else {
				if (fpcgrasses_c3 > fpcgrasses_c4)
					gridcell.IFStypelow = 2; // Short grass 2
				else
					gridcell.IFStypelow = 7; // Tall grass 7
			}
		}
	}
} // computeIFSVegetationCover



std::string lpjg_to_string(int d) {
	std::ostringstream os;
	os << d;
	return os.str();
}



bool getdailyIFSforcing(const char* varname, const int& startyear, const int numyear, double data[NYEARIFS][366], int daysinyear, int ifsindex, bool depth) {

	string filenametail = "_dayavg.nc";

	// var39_1971_dayavg.nc

	xtring ifs_spinup_dir=param["ifs_spinup_dir"].str;

	int yrix = 0;

	double yearavg = 0.0; // Test of annual averages

	// New code for reading AMIP files
	const char* filevar = varname;

	if (startyear < 1900) {

		// ECE v3.2 AMIP forcing only
		if (!strcmp(varname,"var39"))
			filevar = "var039";
		if (!strcmp(varname, "var40"))
			filevar = "var040";	
		if (!strcmp(varname, "var41"))
			filevar = "var041";	
		if (!strcmp(varname, "var42"))
			filevar = "var042";
	
	}

	// ecev3 - AMIP fix. No AMIP data before 1870, so use that data instead
	int datastartyear = startyear;
	if (startyear < 1870)
		datastartyear = 1870;

	for (int yr = datastartyear; yr < datastartyear+numyear; yr++) {

		// First create NC file name for this year of data 

		string yearstr;
		std::string full_path((char*)ifs_spinup_dir);
		
		// AMIP
		full_path += filevar + std::string("_") + lpjg_to_string(yr) + filenametail;

		// Open NC file
		int status, ncid;
		status = nc_open(full_path.c_str(), NC_NOWRITE, &ncid);
		if (status != NC_NOERR)
			return false;

		// Get the ID for varname
		int var_id;
		status = nc_inq_varid (ncid, varname, &var_id);
		if (status != NC_NOERR)
			return false;

		// check for leap years
		if (!(yr % 4) && (yr % 100 | !(yr % 400)))
			daysinyear = 366;
		else
			daysinyear = 365;

		// Read data
		// double dataarray[365*88834];
		//status = nc_get_var_double(ncid, var_id, dataarray);
		
		for (int dd = 0; dd < daysinyear; dd++) {
			size_t index[] = {dd,0,ifsindex-1}; // NB! because ifsindex has base 1
			size_t depthindex[] = {dd,0,0,ifsindex-1}; // E.g. soil temp & moisture
			
			double val;
			
			if (!depth)
				status = nc_get_var1_double(ncid, var_id, index, &val);
			else 
				status = nc_get_var1_double(ncid, var_id, depthindex, &val);
			
			if (status != NC_NOERR)
				return false;
			data[yrix][dd] = val;

			yearavg+=data[yrix][dd];
		}

		nc_close(ncid);

		// Get the data for variable longitude
		//status = nc_get_var_double(ncid, lon_id, lon);
		//if (status != NC_NOERR) handle_error(status);

		// Read one year of data for this point (ifsindex)

		yrix++;
	}

	yearavg/=(numyear*365);

	return true;

} // getdailyIFSforcing


// Faster version of the above function.
// Read in the data for all NYEARIFS (10) years at once for this point.
// We created one file per variable as follows, for var143:
// cdo mergetime *var143* var143_1870-1879.nc
// ncpdq -a x,time,y var143_1870-1879.nc var143_1870-1879_fast.nc

bool getdailyIFSforcing_fast(const char* varname, const int& startyear, const int numyear, double data[NYEARIFS][366], int daysinyear, int ifsindex, bool depth) {

	string filenametail = "_1870-1879_fast.nc";

	xtring ifs_spinup_dir = param["ifs_spinup_dir"].str;

	// New code for reading AMIP files
	const char* filevar = varname;

	if (startyear < 1900) {

		// ECE v3.2 AMIP forcing only
		if (!strcmp(varname, "var39"))
			filevar = "var039";
		if (!strcmp(varname, "var40"))
			filevar = "var040";
		if (!strcmp(varname, "var41"))
			filevar = "var041";
		if (!strcmp(varname, "var42"))
			filevar = "var042";

	}

	// ecev3 - AMIP fix. No AMIP data before 1870, so use that data instead
	int datastartyear = startyear;
	if (startyear < 1870)
		datastartyear = 1870;

	// First create NC file name for this year of data 

	string yearstr;
	std::string full_path((char*)ifs_spinup_dir);

	// AMIP - merged and reordered
	full_path += filevar + filenametail;
	dprintf("full_path %s \n",full_path.c_str());
	// Open NC file
	int status, ncid;
	status = nc_open(full_path.c_str(), NC_NOWRITE, &ncid);
	if (status != NC_NOERR) {
		fprintf(stderr, "Error opening %s: %s \n",full_path.c_str(),nc_strerror(status));
		return false;
	}
	dprintf("file read \n");
	// Get the ID for varname
	int var_id;
	status = nc_inq_varid(ncid, varname, &var_id);
	if (status != NC_NOERR) {
		fprintf(stderr, "Error inquiring var %s: %s \n",varname,nc_strerror(status));
		return false;
	}

	// Read data
	double data_arr[3652]; // 1870-79 daily data, incl. 2 leap years

	size_t start[] = { ifsindex - 1,0,0 };			// NB! because ifsindex has base 1
	size_t start_depth[] = { ifsindex - 1,0,0,0 };	// NB! because ifsindex has base 1

	size_t count[] = {1,3652,1};					// NB! because ifsindex has base 1
	size_t count_depth[] = {1,3652,1,1};			// NB! because ifsindex has base 1

	if (!depth)
		status = nc_get_vara_double(ncid, var_id, start, count, data_arr);
	else
		status = nc_get_vara_double(ncid, var_id, start_depth, count_depth, data_arr);

	nc_close(ncid);


	int yrix = 0;
	int dayix = 0;
	for (int yr = datastartyear; yr < datastartyear + NYEARIFS; yr++) {

		// check for leap years
		if (!(yr % 4) && (yr % 100 | !(yr % 400)))
			daysinyear = 366;
		else
			daysinyear = 365;

		for (int dd = 0; dd < daysinyear; dd++) {
			data[yrix][dd] = data_arr[dayix];
			dayix++;
		}

		yrix++;
	}

	return true;

} // getdailyIFSforcing_fast


  // Read in ERA20C or GSWP3 data for NYEARIFS (10) years at once for this point.
  // We created one file per variable as follows, for var143:
  // cdo mergetime *var143* var143_1870-1879.nc
  // ncpdq -a x,time,y var143_1901-1910.nc var143_1901-1910_fast.nc

bool getdailyERA20Cforcing_fast(const char* varname, const int& startyear, const int numyear, double data[NYEARIFS][366], int daysinyear, int ifsindex, bool depth) {

	int datastartyear = datastartyear = startyear;
	xtring filenametail;
	filenametail.printf("_%4d-%4d_fast.nc", startyear, startyear+NYEARIFS-1);

	xtring ifs_spinup_dir = param["ifs_spinup_dir"].str;

	const char* filevar = varname;

	// First create NC file name for this year of data 

	string yearstr;
	xtring full_path = ifs_spinup_dir + filevar + filenametail;

	//dprintf("Reading fast forcing file %s for era20c spinup \n", (const char*)full_path);

	// Open NC file
	int status, ncid;
	status = nc_open((const char*)full_path, NC_NOWRITE, &ncid);
	if (status != NC_NOERR)
		return false;

	// Get the ID for varname
	int var_id;
	status = nc_inq_varid(ncid, varname, &var_id);
	if (status != NC_NOERR)
		return false;


	// Read data
	double data_arr[3652]; // 1901-1910 daily data, incl. 2 leap years

	size_t start[] = { ifsindex - 1,0,0 };			// NB! because ifsindex has base 1
	size_t start_depth[] = { ifsindex - 1,0,0,0 };	// NB! because ifsindex has base 1

	size_t count[] = { 1,3652,1 };					// NB! because ifsindex has base 1
	size_t count_depth[] = { 1,3652,1,1 };			// NB! because ifsindex has base 1

	if (!depth)
		status = nc_get_vara_double(ncid, var_id, start, count, data_arr);
	else
		status = nc_get_vara_double(ncid, var_id, start_depth, count_depth, data_arr);

	nc_close(ncid);


	int yrix = 0;
	int dayix = 0;
	for (int yr = datastartyear; yr < datastartyear + NYEARIFS; yr++) {

		// check for leap years
		if (!(yr % 4) && (yr % 100 | !(yr % 400)))
			daysinyear = 366;
		else
			daysinyear = 365;

		for (int dd = 0; dd < daysinyear; dd++) {
			data[yrix][dd] = data_arr[dayix];
			dayix++;
		}

		yrix++;
	}

	return true;

} // getdailyERA20Cforcing_fast


// Used for nearest-neighbour N dep
void getndep_nearest_index(double lon, double lat, int ncid, int& lon_index, int& lat_index) {

	// Search all coordinates in a square around (lon, lat), but first go down to
	// multiple of 0.5 - ecev3 - add 0.25
	double center_lon = floor(lon * 2) / 2 + 0.25;
	double center_lat = floor(lat * 2) / 2 + 0.25;

	int nrlat = 96;
	int nrlon = 144;

	// Get the ID for longitude
	int status, lon_id;
	status = nc_inq_varid(ncid, "lon", &lon_id);
	if (status != NC_NOERR)
		fail("Can't find N deposition longitude");

	double val;
	double diff = 360.0;

	// Loop through latitudes to find the nearest
	for (int lo = 0; lo < nrlon; lo++) {
		size_t index[] = { lo, 0, 1 };
		status = nc_get_var1_double(ncid, lon_id, index, &val);
		
		// same longitude scale 
		if (val >= 180.0) 
			val = val - 360.0;
		
		// determine if this latitude is nearer
		if (abs((val + 1.25) - center_lon) < diff) {
			lon_index = lo;
			diff = abs((val + 1.25) - center_lon);
		}
	}

	// Get the ID for latitude
	int lat_id;
	status = nc_inq_varid(ncid, "lat", &lat_id);
	if (status != NC_NOERR)
		fail("Can't find N deposition latitude");

	double val1, val1_c, val2;
	diff = 180.0;
	lat_index = -999;

	// Loop through latitudes to find the nearest
	size_t index[] = { 0, 0, 1 };
	status = nc_get_var1_double(ncid, lat_id, index, &val1);
	for (int la = 1; la < nrlat; la++) {
		size_t index[] = { la, 0, 1 };
		status = nc_get_var1_double(ncid, lat_id, index, &val2);
		
		// calculate center for previous latitude cell
		val1_c = val1 - (val1 - val2) / 2.0;

		// determine if the previous latitude is nearer
		if (abs(val1_c - center_lat) < diff) {
			lat_index = la - 1;
			diff = abs(val1_c - center_lat);
		}
		val1 = val2;

		// if the end if reached without updating the latitude index, 
		// then the last latitude cell is the correct one.
		if (la == 95 && lat_index < 0)
			lat_index = la;
	}
}


// Used for nearest-neighbour N dep
bool getmonthlyNdepforcing(const char* varname, const int numyear, double data[NYEARNDEP][12], double lon, double lat, int& ndeplonindex, int& ndeplatindex) {

	string filenametail = "_input4MIPs_surfaceFluxes_CMIP_NCAR-CCMI-1-0_gr_185001-201412.nc";

	xtring ndep_cmip6_dir = param["ndep_cmip6_dir"].str;

	const char* filevar = varname;

	if (!strcmp(varname, "drynhx"))
		filevar = "drynhx";
	if (!strcmp(varname, "drynoy"))
		filevar = "drynoy";
	if (!strcmp(varname, "wetnhx"))
		filevar = "wetnhx";
	if (!strcmp(varname, "wetnoy"))
		filevar = "wetnoy";

	// First create NC file name for this data 
	std::string full_path((char*)ndep_cmip6_dir);

	// AMIP
	full_path += filevar + filenametail;

	// Open NC file
	int status, ncid;
	status = nc_open(full_path.c_str(), NC_NOWRITE, &ncid);
	if (status != NC_NOERR)
		return false;

	if (ndeplonindex < 0)
		getndep_nearest_index(lon, lat, ncid, ndeplonindex, ndeplatindex);

	// Get the ID for varname
	int var_id;
	status = nc_inq_varid(ncid, varname, &var_id);
	if (status != NC_NOERR)
		return false;

	// Read data
	for (int yy = 0; yy < numyear; yy++) {
		for (int mm = 0; mm < 12; mm++) {
			size_t index[] = { yy*12+mm, ndeplatindex - 1, ndeplonindex - 1 }; // NB! because ndep index has base 1

			double val;

			status = nc_get_var1_double(ncid, var_id, index, &val);

			if (status != NC_NOERR)
				return false;
			data[yy][mm] = val * 86400.0; // convert from kg N m-2 s-1 to kg N m-2 day-1
		}
	}

	nc_close(ncid);

	return true;

} // getmonthlyNdepforcing



bool getmonthlyNdepforcingConserve(const char* varname, const int numyear, double data[NYEARNDEP][12], int ndepindex) {

	// 1st (".nc") or 2nd order remapping file?
	string filenametail = "2.nc";

	xtring ndep_cmip6_dir = param["ndep_cmip6_dir"].str;

	const char* filevar = varname;

	if (!strcmp(varname, "drynhx"))
		filevar = "drynhx";
	if (!strcmp(varname, "drynoy"))
		filevar = "drynoy";
	if (!strcmp(varname, "wetnhx"))
		filevar = "wetnhx";
	if (!strcmp(varname, "wetnoy"))
		filevar = "wetnoy";

	// First create NC file name for this data 
	std::string full_path((char*)ndep_cmip6_dir);

	// AMIP
	full_path += filevar + filenametail;

	// Open NC file
	int status, ncid;
	status = nc_open(full_path.c_str(), NC_NOWRITE, &ncid);
	if (status != NC_NOERR)
		return false;

	// Get the ID for varname
	int var_id;
	status = nc_inq_varid(ncid, varname, &var_id);
	if (status != NC_NOERR)
		return false;

	// Read data
	for (int yy = 0; yy < numyear; yy++) {
		for (int mm = 0; mm < 12; mm++) {
			size_t index[] = {yy * 12 + mm, ndepindex};

			double val;

			status = nc_get_var1_double(ncid, var_id, index, &val);

			// < 0?
			val = max(val, 0.0);

			if (status != NC_NOERR)
				return false;

			data[yy][mm] = val * 86400.0; // convert from kg N m-2 s-1 to kg N m-2 day-1
		}
	}

	nc_close(ncid);

	return true;

} // getmonthlyNdepforcingConserve




/// Spinup LPJ-GUESS using IFS output data for EC-Earth restarts
/**
 *
 * \param gridcell				The gridcell to simulate
 * \param year					The year to save the state for. e.g. 1850
 * \param prescribe_ifs_soilmoist	Use IFS soil moisture if true, 
 *								but standard LPJ-GUESS parameterisations if false 
 * \param prescribe_ifs_soiltemp	Use IFS soil temperature if true, 
 *								but standard LPJ-GUESS parameterisations if false 
 */
bool IFSspinup(Gridcell& gridcell, int year, double co2spinup, bool prescribe_ifs_soilmoist, 
	       bool prescribe_ifs_soiltemp) {

	bool spinon=true;

	// Read netcdf files with daily ECE output. 
	dprintf("Reading data before LPJ-GUESS IFS spinup \n");

	int repetitions = (int) nyear_spinup / NYEARIFS;
	int repix = 0;
	int calyear = year;
	int yrix = 0;
	int spinupyear = 0;

	// Keep a track of the simulation year for this cell.  
	gridcell.simulationyear = spinupyear;

	// ecev3 - isspinup true will ensure fixed 1850 land use during spinup
	gridcell.isspinup = true; // ONLY true during spinup!


	// Arrays to be populated by data from IFS output
	double ifstemp[NYEARIFS][366];
	double ifsprecip[NYEARIFS][366];
	double ifsprecipstrat[NYEARIFS][366];
	double ifsprecipconv[NYEARIFS][366];
	double ifsswradnet[NYEARIFS][366];	// net SW at the surface
	double ifslwradnet[NYEARIFS][366];	// net LW at the surface
	double ifstempsoil2[NYEARIFS][366];
	double ifstempsoil3[NYEARIFS][366];
	double ifssw1[NYEARIFS][366];
	double ifssw2[NYEARIFS][366];
	double ifssw3[NYEARIFS][366];
	double ifssw4[NYEARIFS][366];

	// Reset date
	gridcell.date.set(year,0,0,0,0,0);
	date = gridcell.date;

	// Populate arrays one by one using gridcell.ifs_index for this gridcell

	// T2M (K)
	const char* vartoread = "var167";
	bool dataOK = getdailyIFSforcing_fast(vartoread, year, NYEARIFS, ifstemp, gridcell.date.year_length(), gridcell.ifs_index, false);
	if (!dataOK) {
		dprintf("Error reading T2M data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SURFACE SOLAR RAD - SW (W m-2 s - average of 4 6h periods)
	vartoread = "var176";
	dataOK = getdailyIFSforcing_fast(vartoread, year, NYEARIFS, ifsswradnet, gridcell.date.year_length(), gridcell.ifs_index, false);
	if (!dataOK) {
		dprintf("Error reading NETSW data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SURFACE THERMAL RAD - LW (W m-2 s - average of 4 6h periods)
	vartoread = "var177";
	dataOK = getdailyIFSforcing_fast(vartoread, year, NYEARIFS, ifslwradnet, gridcell.date.year_length(), gridcell.ifs_index, false);
	if (!dataOK) {
		dprintf("Error reading NETLW data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL TEMP 2 (K)
	vartoread = "var170";
	dataOK = getdailyIFSforcing_fast(vartoread, year, NYEARIFS, ifstempsoil2, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading SOIL TEMP 2 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL TEMP 3 (K)
	vartoread = "var183";
	dataOK = getdailyIFSforcing_fast(vartoread, year, NYEARIFS, ifstempsoil3, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading SOIL TEMP 3 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL WATER 1 (m3 m-3)
	vartoread = "var39";
	dataOK = getdailyIFSforcing_fast(vartoread, year, NYEARIFS, ifssw1, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading SOIL WATER 1 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL WATER 2 (m3 m-3)
	vartoread = "var40";
	dataOK = getdailyIFSforcing_fast(vartoread, year, NYEARIFS, ifssw2, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading SOIL WATER 2 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL WATER 3 (m3 m-3)
	vartoread = "var41";
	dataOK = getdailyIFSforcing_fast(vartoread, year, NYEARIFS, ifssw3, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading SOIL WATER 3 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL WATER 4 (m3 m-3)
	vartoread = "var42";
	dataOK = getdailyIFSforcing_fast(vartoread, year, NYEARIFS, ifssw4, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading SOIL WATER 4 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// STRATIFORM PRECIP (m - average of 4 6h periods). NB! This includes both rain and snow
	vartoread = "var142";
	dataOK = getdailyIFSforcing_fast(vartoread, year, NYEARIFS, ifsprecipstrat, gridcell.date.year_length(), gridcell.ifs_index, false);
	if (!dataOK) {
		dprintf("Error reading STRATIFORM PRECIP data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// CONVECTIVE PRECIP (m - average of 4 6h periods) NB! This includes both rain and snow
	vartoread = "var143";
	dataOK = getdailyIFSforcing_fast(vartoread, year, NYEARIFS, ifsprecipconv, gridcell.date.year_length(), gridcell.ifs_index, false);
	if (!dataOK) {
		dprintf("Error reading CONVECTIVE PRECIP data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}


	// CONVERSIONS and averages:
	// **************************

	float lpjgsw1[NYEARIFS][366];
	float lpjgsw2[NYEARIFS][366];

	// Precip: aggregation and conversion to mm/day (*4 for full day acc precip & * 1000 m to mm) 
	// SW: Conversion to W m-2. (* 4 as to get daily J m-2 (accumulated), / secs_in_day to get W m-2 
	double secs_in_day = 60.0 * 60.0 * 24.0;

	// Annual averages - Forget about leap years
	double t2mavg = 0.0;
	double soiltavg = 0.0;
	double precavg = 0.0;
	double SWavg = 0.0;
	double LWavg = 0.0;
	double lpjg_wcont = 0.0;
	int ndays = 0;


	for (int yy = 0; yy < NYEARIFS; yy++) {
		for (int dd = 0; dd < 366; dd++ ) {
			// Precip:
			ifsprecip[yy][dd] = (ifsprecipstrat[yy][dd] + ifsprecipconv[yy][dd]) * 4 * 1000; // ecev3 bugfix - removed snow
			ifsprecip[yy][dd] = max(ifsprecip[yy][dd],0.0); // remove small -ve values
			if (ifsprecip[yy][dd] < 0.000001) ifsprecip[yy][dd] = 0.0; // remove very small or negative precip amounts (< 0.000001 mm/day)

			// SW:
			ifsswradnet[yy][dd] *= 4.0 / secs_in_day;
			if (ifsswradnet[yy][dd] < 0.000001) ifsswradnet[yy][dd] = 0.0; // remove very small or negative SSR amounts (< 0.000001 W m-2)
			
			// LW:
			ifslwradnet[yy][dd] *= 4.0 / secs_in_day; // can be positive or negative

			// Soil water
			ifs_to_lpjg_soilwater(ifssw1[yy][dd], ifssw2[yy][dd], ifssw3[yy][dd],
				ifssw4[yy][dd], lpjgsw1[yy][dd], lpjgsw2[yy][dd]);

			if (dd < 365) {
				ndays++;
				t2mavg += ifstemp[yy][dd];
				SWavg += ifsswradnet[yy][dd];
				LWavg += ifslwradnet[yy][dd];
				precavg += ifsprecip[yy][dd];
				soiltavg += ifstempsoil2[yy][dd];
				lpjg_wcont += lpjgsw1[yy][dd];
			}
		}
	}

	t2mavg /= ndays;
	precavg /= ndays;
	SWavg /= ndays;
	LWavg /= ndays;
	soiltavg /= ndays;
	lpjg_wcont /= ndays;


	// Now loop through the X years until Y years of spinup have been performed
	dprintf("Data loaded. Now performing LPJ-GUESS spinup \n");
	dprintf("Spinup mean T2M, Prec, SW, LW, SoilT, SoilW for cell with lon=%f, lat=%f are: %f, %f, %f, %f, %f, %f \n", gridcell.get_lon(), gridcell.get_lat(), 
		t2mavg, precavg, SWavg, LWavg, soiltavg, lpjg_wcont);

	// CO2 now passed into this routine
	// double co2spinup = 284.725; // ppm - fixed for spinup period - preIndustrial, 1850
	// double co2spinup = 353.855; // ppm - fixed for spinup period - modern, 1990
	// double co2spinup = 324.985; // ppm - fixed for spinup period - modern, 1970

	//////////////// ERROR fix N dep index, nearest cell
	gridcell.ndep_lon_index = -999;
	gridcell.ndep_lat_index = -999;

	// Monthly dry NH4 deposition (kg m-2 s-1)
	vartoread = "drynhx";
	double drynhx_data[NYEARNDEP][12];
	// dataOK = getmonthlyNdepforcing(vartoread, NYEARNDEP, drynhx_data, gridcell.get_lon(), gridcell.get_lat(), gridcell.ndep_lon_index, gridcell.ndep_lat_index);
	dataOK = getmonthlyNdepforcingConserve(vartoread, NYEARNDEP, drynhx_data, gridcell.ndep_index);
	if (!dataOK) {
		dprintf("Error reading dry NHx deposition data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// Monthly wet NH4 deposition (kg m-2 s-1)
	vartoread = "wetnhx";
	double wetnhx_data[NYEARNDEP][12];
	//dataOK = getmonthlyNdepforcing(vartoread, NYEARNDEP, wetnhx_data, gridcell.get_lon(), gridcell.get_lat(), gridcell.ndep_lon_index, gridcell.ndep_lat_index);
	dataOK = getmonthlyNdepforcingConserve(vartoread, NYEARNDEP, wetnhx_data, gridcell.ndep_index);
	if (!dataOK) {
		dprintf("Error reading wet NHx deposition data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// Monthly dry NO3 deposition (kg m-2 s-1)
	vartoread = "drynoy";
	double drynoy_data[NYEARNDEP][12];
	//dataOK = getmonthlyNdepforcing(vartoread, NYEARNDEP, drynoy_data, gridcell.get_lon(), gridcell.get_lat(), gridcell.ndep_lon_index, gridcell.ndep_lat_index);
	dataOK = getmonthlyNdepforcingConserve(vartoread, NYEARNDEP, drynoy_data, gridcell.ndep_index);
	if (!dataOK) {
		dprintf("Error reading dry NOy deposition data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// Monthly wet NO3 deposition (kg m-2 s-1)
	vartoread = "wetnoy";
	double wetnoy_data[NYEARNDEP][12];
	//dataOK = getmonthlyNdepforcing(vartoread, NYEARNDEP, wetnoy_data, gridcell.get_lon(), gridcell.get_lat(), gridcell.ndep_lon_index, gridcell.ndep_lat_index);
	dataOK = getmonthlyNdepforcingConserve(vartoread, NYEARNDEP, wetnoy_data, gridcell.ndep_index);
	if (!dataOK) {
		dprintf("Error reading wet NOy deposition data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// Get monthly ndep values and convert to daily
	double mndrydep[12];
	double mnwetdep[12];

	/// Daily N deposition for current year
	double dndep[Date::MAX_YEAR_LENGTH];

	bool veryfirstday = true; // needed to initialise below
	bool printoutput = false; // Will be set to true for the final 50 years of the spinup simulation
	bool printdlai = false;	// print daily LAI on the last day of the spinup.
	FILE *ofp; // Daily LAI

	while (spinon) {

		// START OF LOOP THROUGH SIMULATION DAYS

		// First get the climate and soil drivers for this day of the simulation:
		// temp, precip, netswrad, temp_soil,soilw_surf,soilw_deep
		double temp = 0.0;
		double precip = 0.0;
		double netswrad = 0.0;
		double netlwrad = 0.0;
		double temp_soil = 0.0;
		double temp_soil_2 = 0.0;
		double temp_soil_3 = 0.0;
		double soilw_surf = 0.0;
		double soilw_deep = 0.0;

		// T2M
		temp = ifstemp[yrix][gridcell.date.day];
		temp-=K2degC; // Convert to degC

		// PRECIP, NET SW & LW
		precip = ifsprecip[yrix][gridcell.date.day];
		netswrad = ifsswradnet[yrix][gridcell.date.day];
		netlwrad = ifslwradnet[yrix][gridcell.date.day];

		if (veryfirstday) {

			// ecev3 - Only need to populate the dndep array ONCE during the spinup

			// When using .bin files:
			// gridcell.ndepts.get_one_calendar_year(gridcell.date.year, mndrydep, mnwetdep);

			// Use CMIP6 N dep historical data
			for (int mm = 0; mm < 12; mm++) {
				mndrydep[mm] = drynhx_data[0][mm] + drynoy_data[0][mm]; // values for year 1850 during spinup
				mnwetdep[mm] = wetnhx_data[0][mm] + wetnoy_data[0][mm]; // values for year 1850 during spinup
			}

			// Distribute N deposition
			distribute_ndep(mndrydep, mnwetdep, ifsprecip[yrix], gridcell.date.ndaymonth, dndep);
		}

		// SOIL TEMP
		temp_soil_2 = ifstempsoil2[yrix][gridcell.date.day];
		temp_soil_3 = ifstempsoil3[yrix][gridcell.date.day];
		// Interpolate soil temp to 25cm
		temp_soil = ifs_to_lpjg_soiltemp(temp_soil_2,temp_soil_3);
		temp_soil-=K2degC; // Convert to degC

		// SOIL WATER
		soilw_surf = (double)lpjgsw1[yrix][gridcell.date.day];
		soilw_deep = (double)lpjgsw2[yrix][gridcell.date.day];

		// Convert from m3 m-3 to AWC [0-1] depending on soil type
		swConvert(gridcell.soiltype,soilw_surf,soilw_deep);

		// Daily N deposition (Read from file)
		// double dndep  = ndep / (365.0 * 10000.0); // units kgN/ha/year to kgN/m2/day
		// double dndep_today  = ndep / (365.0 * 10000.0); // units kgN/ha/year to kgN/m2/day
		double dndep_today  = dndep[gridcell.date.day]; // units: kgN/m2/day
		double dnfert = 0.0;

		// NB! ecev3 - Could call getgridcell here if LU cover changes during spinup...

		// Transfer today's climate data to the Climate object in Gridcell.
		gridcell.climate.setdrivers(temp, precip, netswrad, netlwrad, co2spinup, temp_soil, soilw_surf, soilw_deep, dndep_today, dnfert);

		// Simulate vegetation today!
		simulate_day_ece(gridcell, prescribe_ifs_soilmoist, prescribe_ifs_soiltemp, 
				 veryfirstday, printoutput);
		veryfirstday= false; // NB! gridcell.simulationyear is updated on Dec 31 in simulate_day_ece

		// Daily LAI?
		if (spinupyear == nyear_spinup - 1 && gridcell.date.day == 0) {
			printdlai = true;
			// Print out dLAI?
			ofp = fopen("dailyLAIandFPC.txt", "a");
		}

		if (printdlai) {
			// Now calculate the LAI for high and low fractions
			computeDailyLAIandFPCforIFS(gridcell);
			fprintf(ofp, "%8.6f\t %8.6f\t %d\t %d\t %8.5f\t %8.5f\t %8.5f\t %8.5f\n", gridcell.get_lon(), gridcell.get_lat(), date.get_calendar_year(), date.day + 1,
				gridcell.laiphen_low_today, gridcell.laiphen_high_today, gridcell.IFSfraclow,gridcell.IFSfrachigh);
		}


		// Advance day and check for new years, number of repetitions etc.

		int oldyear = gridcell.date.year;
		// Next day of the simulation
		gridcell.date.next();
		date=gridcell.date;
		int newyear = gridcell.date.year;

		if (newyear != oldyear) {
			// If Jan 1...
			calyear++;
			yrix++;
			spinupyear++;

			// if (spinupyear >= 0) // for testing & debugging
			if (spinupyear >= nyear_spinup-5)
				printoutput = true; // Only print output if we're in the final 10 years of the spinup

			// test with: 
			// dprintf("Repetition %4d... \n", spinupyear);
			// dprintf("Nr of stands %4d\n", gridcell.nbr_stands());

			if (calyear == year + NYEARIFS) {
				repix++;		// One more repetition completed...
				calyear = year; // so reset the clock 
				gridcell.date.set(year,0,0,0,0,0);
				date=gridcell.date;
				yrix=0;			// Refer to correct array elements above
			}
		}

		if (repix==repetitions) // End of spinup if we've run all repetitions
			spinon=false;
	}

	// Close LAI file
	if (printdlai) {
		fclose(ofp);
	}


	dprintf("LPJ-GUESS IFS spinup complete\n");

	gridcell.isspinup = false; // set to false 
	
	return true;

} // IFSspinup




  /// Spinup LPJ-GUESS using ERA20C or GSWP3 output data for EC-Earth restarts
  /**
  *
  * \param gridcell				The gridcell to simulate
  * \param year					The year to save the state for. e.g. 1850
  * \param prescribe_ifs_soilmoist	Use IFS soil moisture if true,
  *								but standard LPJ-GUESS parameterisations if false
  * \param prescribe_ifs_soiltemp	Use IFS soil temperature if true,
  *								but standard LPJ-GUESS parameterisations if false
  */
bool ERA20Cspinup(Gridcell& gridcell, int year, double co2spinup, bool prescribe_ifs_soilmoist,
	bool prescribe_ifs_soiltemp) {

	bool spinon = true;

	// Read netcdf files with daily ECE output. 
	dprintf("Reading data before LPJ-GUESS ERA20C spinup \n");

	int repetitions = (int)nyear_spinup / NYEARIFS;
	int repix = 0;
	int calyear = year;
	int yrix = 0;
	int spinupyear = 0;

	// Keep a track of the simulation year for this cell.  
	gridcell.simulationyear = spinupyear;

	// ecev3 - isspinup true will ensure fixed 1850 land use during spinup
	gridcell.isspinup = true; // ONLY true during spinup!


							  // Arrays to be populated by data from IFS output
	double ifstemp[NYEARIFS][366];
	double ifsprecip[NYEARIFS][366];
	double ifsswradnet[NYEARIFS][366];	// net SW at the surface
	double ifslwradnet[NYEARIFS][366];	// net LW at the surface
	double ifstempsoil2[NYEARIFS][366];
	double ifstempsoil3[NYEARIFS][366];
	double ifssw1[NYEARIFS][366];
	double ifssw2[NYEARIFS][366];
	double ifssw3[NYEARIFS][366];
	double ifssw4[NYEARIFS][366];

	// Reset date
	gridcell.date.set(year, 0, 0, 0, 0, 0);
	date = gridcell.date;

	// Populate arrays one by one using gridcell.ifs_index for this gridcell

	// T2M (K)
	const char* vartoread = "var167";
	bool dataOK = getdailyERA20Cforcing_fast(vartoread, year, NYEARIFS, ifstemp, gridcell.date.year_length(), gridcell.ifs_index, false);
	if (!dataOK) {
		dprintf("Error reading ERA20C T2M data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SURFACE SOLAR RAD - SW (W m-2)
	vartoread = "var176";
	dataOK = getdailyERA20Cforcing_fast(vartoread, year, NYEARIFS, ifsswradnet, gridcell.date.year_length(), gridcell.ifs_index, false);
	if (!dataOK) {
		dprintf("Error reading ERA20C NETSW spinup data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SURFACE THERMAL RAD - LW (W m-2)
	vartoread = "var177";
	dataOK = getdailyERA20Cforcing_fast(vartoread, year, NYEARIFS, ifslwradnet, gridcell.date.year_length(), gridcell.ifs_index, false);
	if (!dataOK) {
		dprintf("Error reading ERA20C NETLW data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL TEMP 2 (K)
	vartoread = "var170";
	dataOK = getdailyERA20Cforcing_fast(vartoread, year, NYEARIFS, ifstempsoil2, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading ERA20C SOIL TEMP 2 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL TEMP 3 (K)
	vartoread = "var183";
	dataOK = getdailyERA20Cforcing_fast(vartoread, year, NYEARIFS, ifstempsoil3, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading ERA20C SOIL TEMP 3 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL WATER 1 (m3 m-3)
	vartoread = "var39";
	dataOK = getdailyERA20Cforcing_fast(vartoread, year, NYEARIFS, ifssw1, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading ERA20C SOIL WATER 1 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL WATER 2 (m3 m-3)
	vartoread = "var40";
	dataOK = getdailyERA20Cforcing_fast(vartoread, year, NYEARIFS, ifssw2, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading ERA20C SOIL WATER 2 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL WATER 3 (m3 m-3)
	vartoread = "var41";
	dataOK = getdailyERA20Cforcing_fast(vartoread, year, NYEARIFS, ifssw3, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading ERA20C SOIL WATER 3 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// SOIL WATER 4 (m3 m-3)
	vartoread = "var42";
	dataOK = getdailyERA20Cforcing_fast(vartoread, year, NYEARIFS, ifssw4, gridcell.date.year_length(), gridcell.ifs_index, true);
	if (!dataOK) {
		dprintf("Error reading ERA20C SOIL WATER 4 data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// PRECIP (mm). NB! This includes both rain and snow
	vartoread = "var228";
	dataOK = getdailyERA20Cforcing_fast(vartoread, year, NYEARIFS, ifsprecip, gridcell.date.year_length(), gridcell.ifs_index, false);
	if (!dataOK) {
		dprintf("Error reading ERA20C PRECIP data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}


	// CONVERSIONS and averages:
	// **************************

	float lpjgsw1[NYEARIFS][366];
	float lpjgsw2[NYEARIFS][366];

	// Precip: aggregation and conversion to mm/day (*4 for full day acc precip & * 1000 m to mm) 
	// SW: Conversion to W m-2. (* 4 as to get daily J m-2 (accumulated), / secs_in_day to get W m-2 
	double secs_in_day = 60.0 * 60.0 * 24.0;

	// Annual averages - Forget about leap years
	double t2mavg = 0.0;
	double soiltavg = 0.0;
	double precavg = 0.0;
	double SWavg = 0.0;
	double LWavg = 0.0;
	double lpjg_wcont = 0.0;
	int ndays = 0;


	for (int yy = 0; yy < NYEARIFS; yy++) {
		for (int dd = 0; dd < 366; dd++) {
			// Precip:
			ifsprecip[yy][dd] = max(ifsprecip[yy][dd], 0.0); // remove small -ve values
			if (ifsprecip[yy][dd] < 0.000001) ifsprecip[yy][dd] = 0.0; // remove very small precip amounts (< 0.000001 mm/day)

			// SW
			if (ifsswradnet[yy][dd] < 0.000001) ifsswradnet[yy][dd] = 0.0; // remove very small or negative SSR amounts (< 0.000001 W m-2)

			// Soil water
			ifs_to_lpjg_soilwater(ifssw1[yy][dd], ifssw2[yy][dd], ifssw3[yy][dd], ifssw4[yy][dd], lpjgsw1[yy][dd], lpjgsw2[yy][dd]);

			if (dd < 365) {
				ndays++;
				t2mavg += ifstemp[yy][dd];
				SWavg += ifsswradnet[yy][dd];
				LWavg += ifslwradnet[yy][dd];
				precavg += ifsprecip[yy][dd];
				soiltavg += ifstempsoil2[yy][dd];
				lpjg_wcont += lpjgsw1[yy][dd];
			}
		}
	}

	t2mavg /= ndays;
	precavg /= ndays;
	SWavg /= ndays;
	LWavg /= ndays;
	soiltavg /= ndays;
	lpjg_wcont /= ndays;


	// Now loop through the X years until Y years of spinup have been performed
	dprintf("Data loaded. Now performing LPJ-GUESS ERA20C spinup \n");
	dprintf("Spinup mean T2M, Prec, SW, LW, SoilT, SoilW for cell with lon=%f, lat=%f are: %f, %f, %f, %f, %f, %f \n", gridcell.get_lon(), gridcell.get_lat(),
		t2mavg, precavg, SWavg, LWavg, soiltavg, lpjg_wcont);

	// CO2 now passed into this routine
	// double co2spinup = 284.725; // ppm - fixed for spinup period - preIndustrial, 1850
	// double co2spinup = 353.855; // ppm - fixed for spinup period - modern, 1990
	// double co2spinup = 324.985; // ppm - fixed for spinup period - modern, 1970

	//////////////// ERROR fix N dep index, nearest cell
	gridcell.ndep_lon_index = -999;
	gridcell.ndep_lat_index = -999;

	// Monthly dry NH4 deposition (kg m-2 s-1)
	vartoread = "drynhx";
	double drynhx_data[NYEARNDEP][12];
	dataOK = getmonthlyNdepforcingConserve(vartoread, NYEARNDEP, drynhx_data, gridcell.ndep_index);
	if (!dataOK) {
		dprintf("Error reading dry NHx deposition data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// Monthly wet NH4 deposition (kg m-2 s-1)
	vartoread = "wetnhx";
	double wetnhx_data[NYEARNDEP][12];
	dataOK = getmonthlyNdepforcingConserve(vartoread, NYEARNDEP, wetnhx_data, gridcell.ndep_index);
	if (!dataOK) {
		dprintf("Error reading wet NHx deposition data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// Monthly dry NO3 deposition (kg m-2 s-1)
	vartoread = "drynoy";
	double drynoy_data[NYEARNDEP][12];
	dataOK = getmonthlyNdepforcingConserve(vartoread, NYEARNDEP, drynoy_data, gridcell.ndep_index);
	if (!dataOK) {
		dprintf("Error reading dry NOy deposition data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// Monthly wet NO3 deposition (kg m-2 s-1)
	vartoread = "wetnoy";
	double wetnoy_data[NYEARNDEP][12];
	dataOK = getmonthlyNdepforcingConserve(vartoread, NYEARNDEP, wetnoy_data, gridcell.ndep_index);
	if (!dataOK) {
		dprintf("Error reading wet NOy deposition data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// Get monthly ndep values and convert to daily
	double mndrydep[12];
	double mnwetdep[12];

	/// Daily N deposition for current year
	double dndep[Date::MAX_YEAR_LENGTH];

	bool veryfirstday = true; // needed to initialise below
	bool printoutput = false; // Will be set to true for the final 50 years of the spinup simulation
	bool printdlai = false;	// print daily LAI on the last day of the spinup.
	FILE *ofp; // Daily LAI

	while (spinon) {

		// START OF LOOP THROUGH SIMULATION DAYS

		// First get the climate and soil drivers for this day of the simulation:
		// temp, precip, netswrad, temp_soil,soilw_surf,soilw_deep
		double temp = 0.0;
		double precip = 0.0;
		double netswrad = 0.0;
		double netlwrad = 0.0;
		double temp_soil = 0.0;
		double temp_soil_2 = 0.0;
		double temp_soil_3 = 0.0;
		double soilw_surf = 0.0;
		double soilw_deep = 0.0;

		// T2M
		temp = ifstemp[yrix][gridcell.date.day];
		temp -= K2degC; // Convert to degC

		// PRECIP, NET SW & LW
		precip = ifsprecip[yrix][gridcell.date.day];
		precip = precip*1.e3*secs_in_day/1000; // Convert from kg m-2 s-1 to mm
		netswrad = ifsswradnet[yrix][gridcell.date.day];
		netlwrad = ifslwradnet[yrix][gridcell.date.day];

		if (veryfirstday) {

			// ecev3 - Only need to populate the dndep array ONCE during the spinup

			// When using .bin files:
			// gridcell.ndepts.get_one_calendar_year(gridcell.date.year, mndrydep, mnwetdep);

			// Use CMIP6 N dep historical data
			for (int mm = 0; mm < 12; mm++) {
				mndrydep[mm] = drynhx_data[0][mm] + drynoy_data[0][mm]; // values for year 1850 during spinup
				mnwetdep[mm] = wetnhx_data[0][mm] + wetnoy_data[0][mm]; // values for year 1850 during spinup
			}

			// Distribute N deposition
			distribute_ndep(mndrydep, mnwetdep, ifsprecip[yrix], gridcell.date.ndaymonth, dndep);
		}

		// SOIL TEMP
		temp_soil_2 = ifstempsoil2[yrix][gridcell.date.day];
		temp_soil_3 = ifstempsoil3[yrix][gridcell.date.day];
		// Interpolate soil temp to 25cm
		temp_soil = ifs_to_lpjg_soiltemp(temp_soil_2, temp_soil_3);
		temp_soil -= K2degC; // Convert to degC

		// SOIL WATER
		soilw_surf = (double)lpjgsw1[yrix][gridcell.date.day];
		soilw_deep = (double)lpjgsw2[yrix][gridcell.date.day];

		// Convert from m3 m-3 to AWC [0-1] depending on soil type
		swConvert(gridcell.soiltype, soilw_surf, soilw_deep);

		// Daily N deposition (Read from file)
		// double dndep  = ndep / (365.0 * 10000.0); // units kgN/ha/year to kgN/m2/day
		// double dndep_today  = ndep / (365.0 * 10000.0); // units kgN/ha/year to kgN/m2/day
		double dndep_today = dndep[gridcell.date.day]; // units: kgN/m2/day
		double dnfert = 0.0;

		// NB! ecev3 - Could call getgridcell here if LU cover changes during spinup...

		// Transfer today's climate data to the Climate object in Gridcell.
		gridcell.climate.setdrivers(temp, precip, netswrad, netlwrad, co2spinup, temp_soil, soilw_surf, soilw_deep, dndep_today, dnfert);

		// Simulate vegetation today!
		simulate_day_ece(gridcell, prescribe_ifs_soilmoist, prescribe_ifs_soiltemp,
			veryfirstday, printoutput);
		veryfirstday = false; // NB! gridcell.simulationyear is updated on Dec 31 in simulate_day_ece

		// Daily LAI?
		if (spinupyear == nyear_spinup - 1 && gridcell.date.day == 0) {
			printdlai = true;
			// Print out dLAI?
			ofp = fopen("dailyLAIandFPC.txt", "a");
		}

		if (printdlai) {
			// Now calculate the LAI for high and low fractions
			computeDailyLAIandFPCforIFS(gridcell);
			fprintf(ofp, "%8.6f\t %8.6f\t %d\t %d\t %8.5f\t %8.5f\t %8.5f\t %8.5f\n", gridcell.get_lon(), gridcell.get_lat(), date.get_calendar_year(), date.day + 1,
				gridcell.laiphen_low_today, gridcell.laiphen_high_today, gridcell.IFSfraclow, gridcell.IFSfrachigh);
		}


		// Advance day and check for new years, number of repetitions etc.

		int oldyear = gridcell.date.year;
		// Next day of the simulation
		gridcell.date.next();
		date = gridcell.date;
		int newyear = gridcell.date.year;

		if (newyear != oldyear) {
			// If Jan 1...
			calyear++;
			yrix++;
			spinupyear++;

			// if (spinupyear >= 0) // for testing & debugging
			if (spinupyear >= nyear_spinup - 5)
				printoutput = true; // Only print output if we're in the final 10 years of the spinup

			if (calyear == year + NYEARIFS) {
				repix++;		// One more repetition completed...
				calyear = year; // so reset the clock 
				gridcell.date.set(year, 0, 0, 0, 0, 0);
				date = gridcell.date;
				yrix = 0;			// Refer to correct array elements above
			}
		}

		if (repix == repetitions) // End of spinup if we've run all repetitions
			spinon = false;
	}

	// Close LAI file
	if (printdlai) {
		fclose(ofp);
	}


	dprintf("LPJ-GUESS ERA20C spinup complete\n");

	gridcell.isspinup = false; // set to false 

	return true;

} // ERA20Cspinup



///////////////////////////////////////////////////////////////////////////////////////
// GUESS_COUPLED
// 

void guess_coupled(int& id,
					int localrank,
					int numlpjgprocesses,
					int isfinal,
					double lon,
					double lat,
					int ifs_soilcd,
					int ifs_index,
					int ndep_index,
					int year,
					int sim_year,
					int month,
					int monthday,
					int hour, 
					double temp, 
					double precip, 
					double netswrad,
					double netlwrad,
					double co2std,
					double temp_soil,
					double soilw_surf,
					double soilw_deep, // Arg18
					double vegl, 
					int vegl_type,
					float& cfluxnattoday,
					float& cfluxanttoday,
					float& npptoday,
					float& laiphen_high,
					float& laiphen_low,
					float& ifsvegfrachigh,
					int& ifsvegtypehigh,				  
					float& ifsvegfraclow,
					int& ifsvegtypelow,
					bool islpjgspinup,
					int fixedNDepafter, 
					int fixedLUafter)
{
	// ARGUMENTS:
	// id					= id from Gridcell object, 0 if initialising
	// is_final				= 0: not the last timestep. 1: the final timestep - so save state!
	// ifs_soilcd			= soil code used in IFS (1-8)
	// ifs_index			= IFS index, i.e. cell 1 to 35718 (T159), or XXXXX (T255)
	// tile					= forest or open land
	// year					= calender! year (e.g. 1964)
	// sim_year				= simulation year (0 to ...)
	// month				= month (0-11)
	// dayofmonth			= what it says (0-30)
	// hour					= hour (0-23)
	// temp					= instantaneous temperature 2m above ground (deg C)
	// precip				= precip (total mm this time step)		
	// netswrad				= net downward shortwave radiation (J/m2/s)
	// netlwrad				= net longwave radiation (J/m2/s)
	// co2std				= ambient CO2 concentration (ppmv)
	// temp_soil			= soil temperature (deg C)
	// soilw_surf			= m3 m-3 from IFS - will convert below to soil water content AWC (0-1) - surface layer 0-50 cm
	// soilw_deep			= m3 m-3 from IFS - will convert below to soil water content AWC (0-1) - deep layer 50-150 cm
	// vegl					= fraction of IFS gridcell covered by the low type
	// vegl_type			= IFS/H-TESSEL low vegetation type - needed to distinguish crops/agriculture, and tundra
	// cfluxnattoday		= natural c clux for the gridcell 
	// cfluxanttoday	       	= anthropogenic c flux for the gridcell 
	// npptoday      	        = npp for the gridcell 
	// laiphen_high			= LAI from high vegetation today
	// laiphen_low			= LAI from low vegetation today
	// ifsvegfrachigh		= fraction of high vegetation in this cell
	//						  updated after spinup and on Dec 31
	// ifsvegtypehigh		= Dominant IFS/H-TESSEL type of high vegetation in this cell
	//						  Values allowed: 3 (EN), 4 (DN), 5 (DB), 6 (EB), 18 (Mixed forest), 19 (Interrupted forest)
	//						  updated after spinup on Dec 31, for FOREST/HIGH stands only. 0 for GRASS/LOW stands
	// ifsvegfraclow		= fraction of low vegetation in this cell
	//						  updated after spinup and on Dec 31
	// ifsvegtypelow		= Dominant IFS/H-TESSEL type of low vegetation in this cell
	//						  Values allowed: 
	//						  updated after spinup on Dec 31
	// islpjgspinup			= Is this a spinup?
	// fixedNDepafter		= year at which N deposition is fixed. Needed for DECK runs. Value -1 if not needed
	// fixedLUafter			= year at which land use is fixed. Needed for DECK runs. Value -1 if not needed

	bool initialising=false;

	// Initialise variables to be returned to IFS
	cfluxnattoday=0.0;
	cfluxanttoday=0.0;
	npptoday=0.0;
	laiphen_high=0.0;
	laiphen_low=0.0;
	ifsvegfrachigh=0.0;
	ifsvegtypehigh=0;
	ifsvegfraclow=0.0;
	ifsvegtypelow=0;


	// Set to TRUE if we want debugguing print statements
	bool LPJG_debug = false;


	// Initialise GUESS if not yet done
	if (!ifinit) {

		// Input and output modules
		const char* input_module_name = "ece"; // ecev3 - use new EC-Earth InputModule
		InputModuleRegistry& theinstance = InputModuleRegistry::get_instance();
		input_module_ece = theinstance.create_input_module(input_module_name);
		
		output_modules_ece = new GuessOutput::OutputModuleContainer();
		GuessOutput::OutputModuleRegistry::get_instance().create_all_modules(*output_modules_ece);

		// Read the instruction file to obtain PFT static parameters and
		// simulation settings
		read_instruction_file("guess.ins");

		print_logfile_heading();

		// Initialise input/output - MUST be before the output module is initialized
 		input_module_ece->init();

		// Initialise output
		(*output_modules_ece).init();

		// Initialise LPJ-GUESS
		initguess(localrank);

		// Create objects for (de)serializing grid cells
		if (!islpjgspinup && restart==1) { // Do not restart from file if we're spinning up.
			dprintf("LPJG: deserialization preparation for state path %s and name %s \n", (char*)state_path, (char*)state_name);
			deserializer_ece.reset(create_deserializer(state_path, state_name, year));
		}

	}

	if (LPJG_debug) 
		dprintf("LPJG DEBUG 1 - Date is YEAR: %04d, MONTH: %02d, DOY:%03d, HOUR:%02d\n",year,month+1,monthday+1,hour);


	if (!id) {
		
		// Create new Gridcell object and initialise
		
	 	if (LPJG_debug) dprintf("LPJG DEBUG 2 - Date is YEAR: %04d, MONTH: %02d, DOY:%03d, HOUR:%02d\n",year,month+1,monthday+1,hour);
 
		initialising=true;

		world.push_back(new Gridcell());
		Gridcell& gridcell = *world.back();
		gridcell.id = world.size();
		gridcell.date.init(1);
		gridcell.ifs_index = ifs_index;		// Record this gridcell's IFS index - used to get IFS spinup data
		gridcell.ndep_index = ndep_index;	// Record this gridcell's N deposition index - used to access N deposition data
		date.init(1);						// date.year needs to be zero when landcover_init/getlandcover is called
			
		// Store id
		id = gridcell.id;

		// Tell framework the coordinates of this grid cell
		gridcell.set_coordinates(lon, lat);

		// The insolation data is sent in as net SW over the whole day, not just daylight hours
		gridcell.climate.instype=NETSWRAD_TS;

		// Tell framework the soil type of this grid cell
		// Set soilcode to the value used by IFS.
		ifssoilparameters(gridcell.soiltype,ifs_soilcd);

		// Initialise dates AFTER call to landcover_init. This ensures that the stands are created.
		// NB - ecev3 - this makes sure that stand.first_year is set in landcover_init!!!
		date.set(year,month,monthday,hour,0,0);
		gridcell.date = date;
		gridcell.set_first_year(date.year); // Now initialise the first year. This ensures that anetps_ff_est_initial is initialised

		// Initialise certain climate and soil drivers
		gridcell.climate.initdrivers(gridcell.get_lat());
		
		// RESTART from state file?
		if (deserializer_ece.get() && !islpjgspinup) {
			deserializer_ece->deserialize_gridcell(gridcell);
			dprintf("LPJG: Loaded state from file\n");

			if (calc_phen_after_restart) {
				// Update phenology - not saved in earlier versions of the code
				// Now it is saved, so there is no need for this code anymore, but it is kept to 
				// ensure for consistency with non-C4MIP CMIP6 runs. Set in guess.ins
				for(unsigned int i = 0; i < gridcell.nbr_stands(); i++) {
				
					// For each stand ...
					Stand& stand = gridcell[i];

					stand.firstobj();
					while (stand.isobj) {

						// For each patch ...
						Patch& patch = stand.getobj();
						leaf_phenology(patch,gridcell.climate);

						stand.nextobj();
					}
				}
			}

		} else
			 dprintf("LPJG: Did not load state from file\n");

		// ecev3 - set a year from which to freeze Land Use Change
		// Moved here, after deserialization because fixedLUafter is serialized. 
		//TODO: remove from serialization
		gridcell.fixedLUafter = fixedLUafter;
		int thisyear = date.get_calendar_year();
		if (gridcell.fixedLUafter >= 0 && thisyear >= gridcell.fixedLUafter) {
			gross_land_transfer = 0;
			if (date.day == 0)
				dprintf("Setting gross_land_transfer = 0. fixedLUfter = %d, year = %d \n", gridcell.fixedLUafter, thisyear);
		}


		// Call input module to obtain updated LU driver data for this grid cell.
		if (!input_module_ece->getgridcell(gridcell)) {
			dprintf("LPJG: error in getgridcell !\n");
			return;
		}

		if(run_landcover && islpjgspinup) {
			// ecev3 - moved from above. We should call this when we're spinning up, and on Jan 1 each year
			// Read static landcover and cft fraction data from ins-file and/or from data files for the spinup period and create stands
			// make sure that landcover is initialized with isspinup = true
			gridcell.isspinup = true;
			landcover_init(gridcell, input_module_ece); // ecev3
		}


		// SPINUP?
		if (islpjgspinup) {

			// Always use IFS data to spinup
			bool spinupOK = false;

			if (!ERA20CSPINUP)
				spinupOK = IFSspinup(gridcell, year, co2std, useifssoilmoist, useifssoiltemp);
			else // ERA20C or GSWP3 forcing
				spinupOK = ERA20Cspinup(gridcell, year, co2std, useifssoilmoist, useifssoiltemp);
			
			if (!spinupOK) {
				 dprintf("Error in spinup \n\n\n");
				 fail();
			}

			// Save the state so we don't need a spinup next time
			dprintf("Rank %i is saving state after spinup for cell with lon=%f, lat=%f\n", localrank, lon, lat);

			if (!serializer_ece.get()) {
				serializer_ece.reset(create_serializer(state_path, state_name, 
					                                   year,
					                                   inst,
					                                   numlpjgprocesses));
			}
			
			serializer_ece->serialize_gridcell(gridcell);		

			// Save the state so we don't need a spinup next time
			dprintf("Rank %i has saved the state after spinup for cell with lon=%f, lat=%f\n", localrank, lon, lat);

		}

		if (LPJG_debug) dprintf("LPJG DEBUG - initialisation OK \n\n\n");

		// End of initialisation
	}


	// No need for return values when spinning up
	if (islpjgspinup) {
		cfluxnattoday=0.0;
		cfluxanttoday=0.0;
		npptoday=0.0;
		laiphen_high=0.0;
		laiphen_low=0.0;
		ifsvegfrachigh=0.0;
		ifsvegtypehigh=0;
		ifsvegfraclow=0.0;
		ifsvegtypelow=0;
		world.clear(); // reduce memory usage during spinup
		return;
	}

	// Retrieve gridcell object
	Gridcell& gridcell = *world[id-1];

	// Reset dates 
	date.set(year,month,monthday,hour,0,0);
	gridcell.date = date;

	// Convert from m3 m-3 to AWC [0-1] depending on soil type
	swConvert(gridcell.soiltype,soilw_surf,soilw_deep);

	/// Daily N deposition for current year
	double dndep[Date::MAX_YEAR_LENGTH];

	if (gridcell.date.day == 0) {
		// Jan 1 only
		// We use the following call to instead read conservatively-mapped N deposition data into 
		// gridcell.mndrydep[12] and gridcell.mnwetdep[12] in getndepCMIP6.
		// Saves the data in gridcell.mndrydep[12] and gridcell.mnwetdep[12], so only need to call this once per year, per gridcell

		bool Ndepdatafound = gridcell.getndepCMIP6(param["ndep_cmip6_dir"].str, fixedNDepafter);
		if (!Ndepdatafound) {
			dprintf("Rank %i could not read N deposition files for cell with lon=%f, lat=%f\n", localrank, lon, lat);
			fail();
		}
	}

	double drizzle[Date::MAX_YEAR_LENGTH];
	for (int dd=0; dd<Date::MAX_YEAR_LENGTH; dd++) drizzle[dd]=1.0; // Assume 1 mm/day

	// Distribute N deposition, assuming it rains a bit every day
	distribute_ndep(gridcell.mndrydep, gridcell.mnwetdep, drizzle, date.ndaymonth, dndep);
	
	double dndep_today = dndep[gridcell.date.day]; // NB! "day" is actually month day
	double dnfert = 0.0;

	// Call input module to obtain LU driver data for this grid cell.
	if (gridcell.date.day == 0) {
		if (!input_module_ece->getgridcell(gridcell)) {
			dprintf("LPJG: error in getgridcell !\n");
			return;
		}
	}


	// Transfer today's climate data to the Climate object in Gridcell.
	gridcell.climate.setdrivers(temp, precip, netswrad, netlwrad, co2std,temp_soil,soilw_surf,soilw_deep, dndep_today, dnfert);

	// simulate vegetation today (and assume that Gridcell properties have already been initialised)
	simulate_day_ece(gridcell, useifssoilmoist, useifssoiltemp, false, true);

	// Don't update date. This is driven outside the call to guess_coupled and updated in Reset dates above
	/*
	gridcell.date.next();
	date=gridcell.date;
	*/


	// Calculate return values
	// -------------

	// 1. C FLUX TODAY

	// NB! a positive C flux (cfluxtoday or npptoday) implies a net C RELEASE by this gridcell, 
	// and a negative value implies a net C UPTAKE

	// Update C flux variable of interest
	unsigned int nrst = gridcell.nbr_stands();
	for(unsigned int i = 0; i < nrst; i++) {
		// For each stand ...
		cfluxnattoday += gridcell[i].dcfluxnat_today * gridcell[i].get_gridcell_fraction();
		cfluxanttoday += gridcell[i].dcfluxant_today * gridcell[i].get_gridcell_fraction();
		npptoday += gridcell[i].dnpp_today * gridcell[i].get_gridcell_fraction();
	}

	// 2. VEGETATION STATE

	// Dominant vegetation type and cover fractions for this grid cell are updated on Dec 31
	// each year, when computeIFSVegetationCover(gridcell) is called from simulate_day_ece() above, 
	// just after outannual

	// Now calculate the LAI for high and low fractions
	computeDailyLAIandFPCforIFS(gridcell);

	// Send the vegetation state information back to IFS
	// Values updated daily:
	laiphen_high		= gridcell.laiphen_high_today;
	laiphen_low			= gridcell.laiphen_low_today;
	ifsvegfrachigh		= gridcell.IFSfrachigh;
	ifsvegfraclow		= gridcell.IFSfraclow;
	// Values updated on Dec 31:
	ifsvegtypehigh		= gridcell.IFStypehigh;
	ifsvegtypelow		= gridcell.IFStypelow;

	/*
	if (date.month==11 && date.dayofmonth==30)
 		dprintf("Year, number of stands: %04d, %04d\n", date.year, gridcell.nbr_stands());
	*/

	// SAVESTATE???
	// -------------

	// Finally, we determine if the state file should be created

	if (isfinal && !initialising) {

        dprintf("Date info for isfinal is: %04d-%02d-%02d, %02d\n",date.year,date.month+1,date.dayofmonth+1);
				
		// Serialize by increasing the year by 1, i.e. so we can restart on Jan 1 the following year
		if (!serializer_ece.get()) {
			serializer_ece.reset(create_serializer(state_path, state_name, year+1, inst, numlpjgprocesses));
		}
		
		serializer_ece->serialize_gridcell(gridcell);
		
	} else {

		// When it's no longer time to save state we'll get rid of the serializer
		// object so that the files get flushed. The next time it's time to
		// save state, we'll create a new serializer object with a new date.
		serializer_ece.reset();
	}

}


// cleans up after run has finished
void guess_coupled_finish()
{
	delete output_modules_ece;
}