ecearth3/sources/lpjg/framework/eceframework.cpp

// Paul Miller, 15 February 2015
// 18 August svn test 
//
#ifdef HAVE_MPI
#include <mpi.h>
#endif

#include "config.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "gutil.h"
#include <math.h>
#include <netcdf.h>
#include "parallel.h"
#include <iostream>

#include <vector>

using namespace GuessParallel;
using namespace std;

#include "shell.h"
#include "oasis-cpp-interface.h"
#include "OasisCoupler.h"
#include <iostream>
#include "commandlinearguments.h"
#include "eceframework.h"
#include "framework.h"


//////////////////////////////////////////////////////////////////////////////////////
// GLOBAL PARAMETERS

// CMIP6 CO2 FORCING

static char file_co2[] = "mole_fraction_of_carbon_dioxide_in_air_input4MIPs_lpjg.nc"; // concatenated CO2 file provided by run_script
// static char file_co2[] = "co2_histo_cmip6_lpjg.nc"; // version converted from the above using CDO as follows:
// cdo -f nc copy -selname,mole_fraction_of_carbon_dioxide_in_air CMIP6_histo_mole_fraction_of_carbon_dioxide_in_air_input4MIPs_gr1-GMNHSH.nc co2_histo_cmip6_lpjg.nc


// NETCDF files with associated variable names

#ifdef GRID_T159 
// ecev3 - T159
// *** Coordinates, masks and gridcell areas

static const int MAXGRID = 10407; //  10407 for 6k and PI onwards, 12245 for LGM
	// number of LAND grid cells at T159 resolution
static const int NX_ATMO=35718;
	// number of GLOBAL grid cells at T159 resolution
static char parn_lon[]="L080.lon";
	// longitude parameter name in filen_grid
static char parn_lat[]="L080.lat";
	// latitude parameter name in filen_grid
static char parn_mask[]="L080.msk";	
	// mask parameter name in filen_mask
static char parn_areas[] = "L080.srf";
	// mask parameter name in filen_areas
static char filen_soilcd[]="slt_T159.nc"; // "slt_T159_LGM.nc" for T159 LGM
	// path to soil data file 
static char parn_lon_cor[] = "L080.clo";
	// corners
static char parn_lat_cor[] = "L080.cla";
	// corners
#endif

#ifdef GRID_T255  
// ecev3 - T255
// *** Coordinates, masks and gridcell areas
static const int MAXGRID=25799; 
	// number of LAND grid cells at T255 resolution
static const int NX_ATMO=88838;
	// number of GLOBAL grid cells at T255 resolution
static char parn_lon[]="L128.lon";
	// longitude parameter name in filen_grid
static char parn_lat[]="L128.lat";
	// latitude parameter name in filen_grid
static char parn_mask[]="L128.msk";
	// mask parameter name in filen_mask
static char parn_areas[]="L128.srf";
	// mask parameter name in filen_areas
static char filen_soilcd[]="slt_T255.nc";
	// path to soil data file
static char parn_lon_cor[] = "L128.clo";
	// corners
static char parn_lat_cor[] = "L128.cla";
	// corners
#endif

// *** Grid definition files (common to all resolutions)
static char filen_grid[] = "grids.nc";
	// path to grid data file
static char filen_areas[]="areas.nc";
	// path to areas data file
static char filen_mask[] = "masks.nc";
	// path to mask data file

// T159 LGM
// static char filen_mask[] = "masks_LGM.nc";
	// path to mask data file


// *** Paths to and variables in netcdf files extracted from ICMGG*INIT files 
static char filen_cvl[]="cvl.nc";
	// path to LOW cover fraction file 
static char parn_cvl[]="var27";
	// soil code parameter name in filen_cvl
static char filen_cvh[]="cvh.nc";
	// path to HIGH cover fraction file 
static char parn_cvh[]="var28";
	// soil code parameter name in filen_cvh
static char filen_tvl[]="tvl.nc";
	// path to LOW cover TYPE file 
static char parn_tvl[]="var29";
	// soil code parameter name in file filen_tvl
static char filen_tvh[]="tvh.nc";
	// path to HIGH cover TYPE file 
static char parn_tvh[]="var30";
	// soil code parameter name in file filen_tvh
static char parn_soilcd[] = "var43";
	// soil code parameter name in filen_soilcd


// *** Initialisation from .rc files
static char file_timesteps[]=".//lpjg_steps.rc";
	// path to IFS file that determines the LPJ-GUESS timesteps, resolution and coupling mode
int resolution; 
	// 255 or 159
bool activateTM5coupling;
	// communicate with TM5 if true
const int NYEAR_ECEARTH=271; 
	// i.e. 2110 - 1840 + 1
	// number of years of historical climate in CRU and files (see below)
const int FIRSTHISTYEAR=1840; 
        // Number of years CO2 is availale
const int NYEAR_CO2=2501;
const int FIRSTYEAR_CO2=0;

bool ifs_A4xCO2 = false;
	// abrupt 4xCO2 after 1850, for DECK Experiments
bool ifs_1PCTCO2 = false;
	// 1% / year (exponential) increase from 1850 on, for DECK Experiments
int ifs_FIXEDYEARCO2 = -1;
	// set CO2 values to chosen year for DECK Experiments
int fixedNDepafter = -1;
	// hold NDep constant from this year on
int fixedLUafter = -1;
	// suppress LU from this year on


int TIMESTEP = 86400;
struct{        
	int year;
	int month;
	int day;
} STARTDATE={1850,1,1}; // initialise with default values
int NTIMESTEP=1;


// *** Date format specification
const bool IFLEAPYEARS=true;
static int NDAYMONTH[12]={31,28,31,30,31,30,31,31,30,31,30,31};
	// if IFLEAPYEARS=true, February will be assigned an extra day when necessary



//////////////////////////////////////////////////////////////////////////////////////
// GLOBAL DATA

static int ngridcell; 
	// number of grid cells

static double co2[NYEAR_CO2];
	// CO2 data for each year of historical data set

static int soilcode[MAXGRID];
	// Soil code for each grid cell
	
// Fields to pass to/from OASIS and i/o fields
static const int NY_ATMO=1;
static float coup_lsmask[NX_ATMO][NY_ATMO];

// Received FROM IFS
static double coup_temper[NX_ATMO*NY_ATMO];
static double coup_precip[NX_ATMO*NY_ATMO];
static double coup_vegl[NX_ATMO*NY_ATMO];
static double coup_vegl_type[NX_ATMO*NY_ATMO];
static double coup_vegh[NX_ATMO*NY_ATMO];
static double coup_snowc[NX_ATMO*NY_ATMO];
static double coup_snowd[NX_ATMO*NY_ATMO];
const int NHTESSELSOILLAYERS = 4; // Number of soil layers in HTESSEL
static double coup_st1l[NX_ATMO*NY_ATMO];
static double coup_st2l[NX_ATMO*NY_ATMO];
static double coup_st3l[NX_ATMO*NY_ATMO];
static double coup_st4l[NX_ATMO*NY_ATMO];
static double coup_sm1l[NX_ATMO*NY_ATMO];
static double coup_sm2l[NX_ATMO*NY_ATMO];
static double coup_sm3l[NX_ATMO*NY_ATMO];
static double coup_sm4l[NX_ATMO*NY_ATMO];
static double coup_sw_radiat[NX_ATMO*NY_ATMO];
static double coup_lw_radiat[NX_ATMO*NY_ATMO];

// Sent TO IFS
static double coup_lowlai[NX_ATMO*NY_ATMO];
static double coup_highlai[NX_ATMO*NY_ATMO];
static double coup_lpjg_typeh[NX_ATMO*NY_ATMO];
static double coup_lpjg_frach[NX_ATMO*NY_ATMO];
static double coup_lpjg_typel[NX_ATMO*NY_ATMO];
static double coup_lpjg_fracl[NX_ATMO*NY_ATMO];

// Received FROM TM5 (when coupled)
static double coup_co2_field[NX_ATMO*NY_ATMO];
// Sent TO TM5 (when coupled)
static double coup_dcflux[NX_ATMO*NY_ATMO];
static double coup_dnpp[NX_ATMO*NY_ATMO];

// Store this year's CO2, either from file or TM5
static double coup_co2;



//////////////////////////////////////////////////////////////////////////////////////
// GRIDLIST

struct EceCoord {

	// Type for storing grid cell properties
	// longitude and latitude
	// ID for the Gridcell
	// soilcode, area, type and cover fractions

	int id;
	float lon;
	float lat;
	int soilcode;	// soilcode
	int IFSindex;	// record the IFS index
    float area;		// gridcell area in m2
	float cvl, cvh;	// Cover fractions from GCM ICMGG*INIT files
	int tvl, tvh;	// Cover types from GCM ICMGG*INIT files
	int ndep_index;	// 0 to NX_ATMO-1 - for accessing N deposition data in CMIP6 files 

	// Only used when used when printing grids
	/*
	double cornerla[4];
	double cornerlo[4];
	*/

	EceCoord() {
		id=0;
		lon=0.0;
		lat=0.0;
		IFSindex=0;
		soilcode=0;
		area=0.0;
		cvl=0.0;
		cvh=0.0;
		tvl=0;
		tvh=0;
		ndep_index=-1;
	};

	~EceCoord() {
	}
};


void handle_error(int status) {
	if (status != NC_NOERR) {
		fprintf(stderr, "%s\n", nc_strerror(status));
		//CLNexit(-1);
    }
}


int opennetcdf(int& ncid, const char* path) {

	try {

		int status;

		status = nc_open(path, NC_NOWRITE, &ncid);
		if (status != NC_NOERR) { 
			handle_error(status);
			return 0;
		}
		else {
			return 1;
		}
	
	} catch (...) {
		cout << "Unknown error in opennetcdf()!" << endl;
		return 0;	
	}
}



int readnetcdf_latlon(const int& ncid, double* lat, double* lon) {

	try{

		int status;

		int lat_id, lon_id;
		char* lonstr = "lon";
		char* latstr = "lat";

		// Get the ID for longitude
		status = nc_inq_varid (ncid, lonstr, &lon_id);
		if (status != NC_NOERR) handle_error(status);

		// Get the ID for latitude
		status = nc_inq_varid (ncid, latstr, &lat_id);
		if (status != NC_NOERR) handle_error(status);
			
		// Get the data for variable longitude
		status = nc_get_var_double(ncid, lon_id, lon);
		if (status != NC_NOERR) handle_error(status);

		// Get the data for variable latitude
		status = nc_get_var_double(ncid, lat_id, lat);
		if (status != NC_NOERR) handle_error(status);

		return 1;

	} catch(...) {
		cout << "Unknown error in readnetcdf_latlon()!" << endl;
		return 0;
	}
}



int getglobaldata_double(const int& ncid, const char* varname, double* vardata) {

	try{

		int status;
		int var_id;

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

		// Get the data for varname
		status = nc_get_var_double(ncid, var_id, vardata);
		if (status != NC_NOERR) handle_error(status);

		return 1;

	} catch(...) {
		cout << "Unknown error in getglobaldata_double()!" << endl;
		return 0;
	}
}


int getglobaldataarray_double(const int& ncid, const char* varname, int corners, int lons, int lats, double* vardata) {

	try{

		int status;
		int var_id;

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

		static size_t start[] = {0, 0, 0}; /* start at first value */
		static size_t count[] = {corners, lons, lats};

		// Get the data for varname
		status = nc_get_vara_double(ncid, var_id, start, count, vardata);
		if (status != NC_NOERR) handle_error(status);

		return 1;

	} catch(...) {
		cout << "Unknown error in getglobaldata_double()!" << endl;
		return 0;
	}
}

int getglobaldata_int(const int& ncid, const char* varname, int* vardata) {

	try{

		int status;
		int var_id;

		// Get the ID for varname
		status = nc_inq_varid (ncid, varname, &var_id);
		if (status != NC_NOERR) handle_error(status);
			
		// Get the data for varname
		status = nc_get_var_int(ncid, var_id, vardata);
		if (status != NC_NOERR) handle_error(status);

		return 1;

	} catch(...) {
		cout << "Unknown error in getglobaldata_int()!" << endl;
		return 0;
	}
}

// ecev3
static ListArray_id<EceCoord> gridlist;
	// Will maintain a list of EceCoord objects containing coordinates
	// and other properties of the grid cells to simulate


// Reads lpjg_steps.rc
bool read_ifs_timesteps() {

	dprintf("Experiment Setup =========== \n");
        FILE* in_ts=fopen(file_timesteps,"rt");
        if (!in_ts) {
			dprintf("Could not open %s for input\n",file_timesteps);
			return false;
        }

        bool eof=false;
		bool is_err=false;
        // Read TIMESTEP
        eof=!readfor(in_ts,"i",&TIMESTEP);

		// Validate TIMESTEP
		if (TIMESTEP>24*3600) {
			dprintf("Maximum allowable timestep is %d seconds (24 hours)\n",24*3600);
			fclose(in_ts);
			return false;
		}

		if (!eof) {

			// Read STARTDATE
			readfor(in_ts,"i",&STARTDATE.year);
			readfor(in_ts,"i",&STARTDATE.month);
			readfor(in_ts,"i",&STARTDATE.day);

			// Validate STARTDATE
			bool valid=true;
			if (STARTDATE.month<1 || STARTDATE.month>12 || STARTDATE.day<1) 
				valid=false;
			else if (IFLEAPYEARS && STARTDATE.month==2 && !(STARTDATE.year%4) && (STARTDATE.year%400)) { // leap year?
				if (STARTDATE.day>NDAYMONTH[1]+1) 
					valid=false;
			}
			else if (STARTDATE.day>NDAYMONTH[STARTDATE.month-1]) 
				valid=false;

			if (!valid) {
				dprintf("Invalid start date %d-%02d-%02d\n", STARTDATE.year,STARTDATE.month,STARTDATE.day);
				fclose(in_ts);
				return false;
			}

			if (STARTDATE.year<0 /* || STARTDATE.year>=FIRSTHISTYEAR+NYEAR_ECEARTH */) { // Larger years can be used in fixed-year experiments
				dprintf("Specified start date %d-%02d-%02d is outside range of available data\n",
					STARTDATE.year,STARTDATE.month,STARTDATE.day);
				fclose(in_ts);
				return false;
			}

			// Read NTIMESTEP
			eof=!readfor(in_ts,"i",&NTIMESTEP);

			// Read resolution
			readfor(in_ts,"i",&resolution);
			if (resolution != 159 && resolution != 255) {
				dprintf("Invalid resolution %d\n", resolution);
				fclose(in_ts);
				return false;
			}

			// Read TM5 coupling option
			int TM5couple = -1;
			readfor(in_ts,"i",&TM5couple);
			if (TM5couple==1)
				activateTM5coupling = true;
			else if (TM5couple==0)
				activateTM5coupling = false;
			else {
				dprintf("Invalid TM5 coupling option %d\n", TM5couple);
				fclose(in_ts);
				return false;
			}

			// Read DECK experiment settings
			// Abrupt 4 x CO2
			xtring sdum[1];
			string is_true_l = "t";
			string is_true_u = "T";
			eof = !readfor(in_ts,"a1",&sdum);
			if ( eof ) {
				dprintf("Entry for A4xCO2 not found in %s \n", file_timesteps);
				is_err = true;
			} 
			else {
				if ( !is_true_l.compare((char *)sdum[0]) || 
				     !is_true_u.compare((char *)sdum[0]) ) {				
					ifs_A4xCO2 = 1;
				}
				else {
					ifs_A4xCO2 = 0;
				}
				dprintf("A4xCO2       set to %i \n",ifs_A4xCO2);
			} 
				
			// 1 percent CO2 per year
			eof = !readfor(in_ts,"a1",&sdum);
			if ( eof ) {
				dprintf("Entry for 1PCTCO2 not found in %s \n", file_timesteps);
				is_err = true;
			}
			else {
				if ( !is_true_l.compare((char *)sdum[0]) || 
				     !is_true_u.compare((char *)sdum[0]) ) {				
					ifs_1PCTCO2 = 1;
				}
				else {
					ifs_1PCTCO2 = 0;
				}			    
				dprintf("1PCTCO2      set to %i \n",ifs_1PCTCO2);
			}
			// fixed year CO2
			eof = !readfor(in_ts,"i",&ifs_FIXEDYEARCO2);
			if ( eof ) {
				dprintf("Entry for FIXEDYEARCO2 not found in %s \n", file_timesteps);
				is_err = true;
			} 
			else {
				dprintf("CO2 is       set fix at %d \n",ifs_FIXEDYEARCO2);
			}
			// fixed NDep from this year on
			eof = !readfor(in_ts,"i",&fixedNDepafter);
			if ( eof ) {
				dprintf("Entry for fixedNDepafter year not found in %s \n", file_timesteps);
				is_err = true;
			} 
			// fixed Land Use from this year on
			eof = !readfor(in_ts,"i",&fixedLUafter);
			if ( eof ) {
				dprintf("Entry for fixedLUafter year not found in %s \n", file_timesteps);
				is_err = true;
			} 
			// Override settings for fixed experiments
			if ( ifs_A4xCO2 || ifs_1PCTCO2 ) {
				fixedNDepafter = CMIP6STARTYEAR;
				fixedLUafter   = CMIP6STARTYEAR;
			}
			dprintf("fixedNDepafter set to %d \n",fixedNDepafter);
			dprintf("fixedLUafter   set to %d \n",fixedLUafter);

		}
		dprintf("End experiment Setup ======= \n");
		fclose(in_ts);
		if ( is_err ) {
			return false;
		}
		else {
			return true;
		}
}


// Reading grid information from files provided by the atmosphere model
// LPJ-GUESS will run on the same grid
// This is where to read grid information from IFS-provided masks.nc, grids.nc and areas.nc

bool read_grid_info(bool readcorners){

	int ix;
	
	double* nlon =new double[NX_ATMO];
	double* nlat =new double[NX_ATMO];
	double* nmask =new double[NX_ATMO];
	double* nareas =new double[NX_ATMO];

	double* nlon_cor = NULL;
	double* nlat_cor = NULL;

	bool ok_flag = true;

	if (readcorners) {
		nlon_cor = new double[4*NX_ATMO];
		nlat_cor = new double[4*NX_ATMO];
	}

	int ncid = 0;
	int gridOK = opennetcdf(ncid, filen_grid);
	int lonsOK = getglobaldata_double(ncid, parn_lon, nlon);
	int latsOK = getglobaldata_double(ncid, parn_lat, nlat);

	if (readcorners) {
		int cloOK = getglobaldataarray_double(ncid, parn_lon_cor, 4, 1, NX_ATMO, nlon_cor);
		int claOK = getglobaldataarray_double(ncid, parn_lat_cor, 4, 1, NX_ATMO, nlat_cor);
		if ( ! cloOK || ! claOK ) {
			dprintf("LPJ-GUESS: Error reading corners from file:%s\n",filen_grid); 
			ok_flag = false;
		}
	}

	nc_close(ncid);

	int ncid2 = 0;
	int maskfileOK = opennetcdf(ncid2, filen_mask);
	int maskOK = getglobaldata_double(ncid2, parn_mask, nmask);
	nc_close(ncid2);

	int ncid3 = 0;
	int areasfileOK = opennetcdf(ncid3, filen_areas);
	int areasOK = getglobaldata_double(ncid3, parn_areas, nareas);
	nc_close(ncid3);

	ngridcell=0;
	int ifsindex = 0;
	
	for(ix=0;ix<NX_ATMO;ix++){
		
		// Record the IFS index (e.g. 1 to 35718 (T159)) in EceCoord now too.	    
		ifsindex++;

		if(nmask[ix]){

			ngridcell++;
			coup_lsmask[0][ix] = nmask[ix]; 
			EceCoord& c=gridlist.createobj(); // add new coordinate to grid list
			c.id=0; // ecev3 - always 0 to begin with, before the LPJG Gridcell is created 
			c.lon=nlon[ix];
			if(c.lon>=180.)c.lon=c.lon-360.;
			c.lat=nlat[ix];

			if (readcorners) {

				/* UNCOMMENT if corners are needed!

				c.cornerla[0]=nlat_cor[ix];
				c.cornerla[1]=nlat_cor[ix+1*NX_ATMO];
				c.cornerla[2]=nlat_cor[ix+2*NX_ATMO];
				c.cornerla[3]=nlat_cor[ix+3*NX_ATMO];

				c.cornerlo[0]=nlon_cor[ix];
				c.cornerlo[1]=nlon_cor[ix+1*NX_ATMO];
				c.cornerlo[2]=nlon_cor[ix+2*NX_ATMO];
				c.cornerlo[3]=nlon_cor[ix+3*NX_ATMO];
				if(c.cornerlo[0]>=180.)c.cornerlo[0]=c.cornerlo[0]-360.;
				if(c.cornerlo[1]>=180.)c.cornerlo[1]=c.cornerlo[1]-360.;
				if(c.cornerlo[2]>=180.)c.cornerlo[2]=c.cornerlo[2]-360.;
				if(c.cornerlo[3]>=180.)c.cornerlo[3]=c.cornerlo[3]-360.;
				*/
			}

			c.IFSindex=ifsindex; // 1 to 35718 (T159) or 88838 (T255)
			c.area = nareas[ix]; // area in m2
			c.ndep_index = ngridcell - 1;
		}
	}


	if (ngridcell>MAXGRID) {
		dprintf("Too many grid cells in gridlist file %s\nMaximum allowed is %d\n", filen_mask,MAXGRID);
		delete[] nlon;
		delete[] nlat;
		delete[] nmask;
		delete[] nareas;

		if (readcorners) {
			delete[] nlon_cor;
			delete[] nlat_cor;
		}
		ok_flag = false;
		//CLNexit(99);
	}
	
	delete[] nlon;
	delete[] nlat;
	delete[] nmask;
	delete[] nareas;
	
	if (readcorners) {
		delete[] nlon_cor;
		delete[] nlat_cor;
	}

	// Error handling cheap way
	if ( ! gridOK ) {
		dprintf("LPJ-GUESS: Error opening file:%s\n",filen_grid); 
	}
	else {
		if ( ! lonsOK ) 
			dprintf("LPJ-GUESS: Error reading lons from file:%s\n",filen_grid); 
		if ( ! latsOK ) 
			dprintf("LPJ-GUESS: Error reading lats from file:%s\n",filen_grid);
	} 

	if ( ! maskfileOK ) 
		dprintf("LPJ-GUESS: Error opening file:%s\n",filen_mask); 
	else if ( ! maskOK ) 
		dprintf("LPJ-GUESS: Error reading mask from file:%s\n",filen_mask); 
		
	if ( ! areasfileOK ) 
		dprintf("LPJ-GUESS: Error opening file:%s\n",filen_areas); 
	else if ( ! areasOK ) 
		dprintf("LPJ-GUESS: Error reading areas from file:%s\n",filen_areas); 
		
	if ( ! gridOK || ! lonsOK || ! latsOK || ! maskfileOK || ! maskOK || ! areasfileOK || ! areasOK ) 
		ok_flag = false;


	return ok_flag;
}


bool read_soil_veg_info(bool printgridlist) {

	// Reading soil and vegetation information from NetCDF input files
	// Will also print the gridlist if the outcommented code is used.

	int ix;
	double* dcoup_soilcd =new double[NX_ATMO];
	double* dcvl =new double[NX_ATMO];
	double* dcvh =new double[NX_ATMO];
	double* dtvl =new double[NX_ATMO];
	double* dtvh =new double[NX_ATMO];

	int ncid = 0;

	bool ok_flag = true;

	// soilcode
	int soilfileOK = opennetcdf(ncid, filen_soilcd);
	int soilOK = getglobaldata_double(ncid, parn_soilcd, dcoup_soilcd);
	nc_close(ncid);
	if ( ! soilfileOK ) 
		dprintf("LPJ-GUESS: Error opening file:%s\n",filen_soilcd); 
	else if ( ! soilOK )
		dprintf("LPJ-GUESS: Error reading from file:%s\n",filen_soilcd); 


	if (printgridlist) {

		// cvl
		int cvlfileOK = opennetcdf(ncid, filen_cvl);
		int cvlOK = getglobaldata_double(ncid, parn_cvl, dcvl);
		nc_close(ncid);
		if ( ! cvlfileOK )
			dprintf("LPJ-GUESS: Error opening file:%s\n",filen_cvl); 
		else if ( ! cvlOK )
			dprintf("LPJ-GUESS: Error reading from file:%s\n",filen_cvl); 

		// cvh
		int cvhfileOK = opennetcdf(ncid, filen_cvh);
		int cvhOK = getglobaldata_double(ncid, parn_cvh, dcvh);
		nc_close(ncid);
		if ( ! cvhfileOK )
			dprintf("LPJ-GUESS: Error opening file:%s\n",filen_cvh); 
		else if ( ! cvhOK )
			dprintf("LPJ-GUESS: Error reading from file:%s\n",filen_cvh); 

		// tvl
		int tvlfileOK = opennetcdf(ncid, filen_tvl);
		int tvlOK = getglobaldata_double(ncid, parn_tvl, dtvl);
		nc_close(ncid);
		if ( ! tvlfileOK )
			dprintf("LPJ-GUESS: Error opening file:%s\n",filen_tvl); 
		else if ( ! tvlOK )
			dprintf("LPJ-GUESS: Error reading from file:%s\n",filen_tvl); 

		// tvh
		int tvhfileOK = opennetcdf(ncid, filen_tvh);
		int tvhOK = getglobaldata_double(ncid, parn_tvh, dtvh);
		nc_close(ncid);
		if ( ! tvhfileOK )
			dprintf("LPJ-GUESS: Error opening file:%s\n",filen_tvh); 
		else if ( ! tvhOK )
			dprintf("LPJ-GUESS: Error reading from file:%s\n",filen_tvh); 

	}

	// Print out gridlist?
	
	FILE *ofp;

	if (printgridlist)
		ofp = fopen("ece_gridlist_T255.txt", "w");

	// ofp = fopen("ece_gridlist_T159.txt", "w"); // T159
	// ofp = fopen("ece_gridlist_LGM.txt", "w"); // T159 - LGM


	// storing soil code in Coord c
	gridlist.firstobj();

	int mask_index = 0;

	for(ix=0;ix<NX_ATMO;ix++){
		if(coup_lsmask[0][ix]==1.0) {

			EceCoord& c=gridlist.getobj();			
			int cd = (int)dcoup_soilcd[ix];
			if(cd <= 0 || cd > 7){
				dprintf("Bad soil code! \n");
				return false;
			}
			c.soilcode=cd;

			if (printgridlist) {
				c.cvl = dcvl[ix];
				c.cvh = dcvh[ix];
				c.tvl = (int)dtvl[ix];
				c.tvh = (int)dtvh[ix];
			}
			else {
				c.cvl = 0.0;
				c.cvh = 0.0;
				c.tvl = 0;
				c.tvh = 0;
			}

			// print gridlist?
			if (printgridlist)
				fprintf(ofp, "%8.5f\t %8.5f\t %8.5f\t %d\t %d\t %d\t %8.5f\t %8.5f\t %d\t %d\t \n", c.lon, c.lat, c.area, c.soilcode, c.IFSindex, mask_index, c.cvl, c.cvh, c.tvl, c.tvh);
			
			gridlist.nextobj();
			mask_index++;
		}
	}

	delete[] dcoup_soilcd;
	delete[] dcvl;
	delete[] dcvh;
	delete[] dtvl;
	delete[] dtvh;

	// print gridlist?
	if (printgridlist)
		fclose(ofp);

	// error handling
	//if ( ! soilfileOK || ! soilOK || ! cvlfileOK || ! cvlOK || ! cvhfileOK || ! cvhOK || 	
	//     ! tvlfileOK || ! tvlOK || ! tvhfileOK || ! tvhOK )

	if (!soilfileOK || !soilOK)
		ok_flag = false;
		
	return ok_flag;
}


static void readco2_cmip6(bool& error_flag) {

	// Reads in atmospheric CO2 concentrations for EC_Earth historical from 
	// Jan 1 0 (1 BC) - to 2014 (see http://www.climate-energy-college.net/cmip6)
	// from netCDF CMIP6 CO2 file. 

	// netCDF related variables
	int status;
	int file_id = 0;

	status = nc_open(file_co2, NC_NOWRITE, &file_id);
	//if (status != NC_NOERR) handle_error(status);
	if (file_id <= 0) {
		dprintf("!!!ERROR!!! readco_cmip6: could not open CO2 file %s for input\n", file_co2);
		error_flag = true;
	}

	// Get dimension time
	int time_id;
	status = nc_inq_unlimdim(file_id, &time_id);
	if (status != NC_NOERR) {
		handle_error(status);
		error_flag = true;
	}
	//time_id = ncdimid(file_id, 'time');
	size_t time_size;
	status = nc_inq_dimlen(file_id, time_id, &time_size);
	if (status != NC_NOERR) {
		handle_error(status);
		error_flag = true;
	}

	// get variable containing co2
	int co2_id;
	char* dsetname = "mole_fraction_of_carbon_dioxide_in_air";
	status = nc_inq_varid(file_id, dsetname, &co2_id);
	if (status != NC_NOERR) {
		handle_error(status);
		error_flag = true;
	}

	// define hyperslab to read from netCDF file
	static size_t start[] = { 0, 0 };
	static size_t count[] = { time_size, 1 };
	float co2_in[NYEAR_CO2];

	// read co2 and assign to global array 
	status = nc_get_vara_float(file_id, co2_id, start, count, co2_in);
	if (status != NC_NOERR) {
		handle_error(status);
		error_flag = true;
	}
	for (int i = 0; i < NYEAR_CO2; i++) {
		co2[i] = (double)co2_in[i];
	}

	// close file
	status = nc_close(file_id);
	if (status != NC_NOERR) {
		handle_error(status);
		error_flag = true;
	}
}


/* if needed when printing cover data too
bool read_soil_veg_cover_info() {

	// Reading soil and vegetation information from NetCDF input files

	int ix;
	double* dcoup_soilcd = new double[NX_ATMO];
	double* dcvl = new double[NX_ATMO];
	double* dcvh = new double[NX_ATMO];
	double* dtvl = new double[NX_ATMO];
	double* dtvh = new double[NX_ATMO];

	int ncid = 0;

	// soilcode
	int soilfileOK = opennetcdf(ncid, filen_soilcd);
	int soilOK = getglobaldata_double(ncid, parn_soilcd, dcoup_soilcd);
	nc_close(ncid);

	// cvl
	int cvfileOK = opennetcdf(ncid, filen_cvl);
	int cvOK = getglobaldata_double(ncid, parn_cvl, dcvl);
	nc_close(ncid);

	// cvh
	cvfileOK = opennetcdf(ncid, filen_cvh);
	cvOK = getglobaldata_double(ncid, parn_cvh, dcvh);
	nc_close(ncid);

	// tvl
	int tvfileOK = opennetcdf(ncid, filen_tvl);
	int tvOK = getglobaldata_double(ncid, parn_tvl, dtvl);
	nc_close(ncid);

	// tvh
	tvfileOK = opennetcdf(ncid, filen_tvh);
	tvOK = getglobaldata_double(ncid, parn_tvh, dtvh);
	nc_close(ncid);

	// vegcover
	// static char filen_vcf[] = "veginfo_T255_con2.nc";
	static char filen_vcf[] = "veginfo_T159_con2.nc";
	// static char filen_vcf[] = "veginfo_T159_bilinear.nc"; // bilinear

	// Categories: bare_ground,Broadleaf,Needleleaf,Evergreen,Decidous,Tree_cover,Herb
	tvfileOK = opennetcdf(ncid, filen_vcf);
	double* dvcf1 = new double[NX_ATMO];
	double* dvcf2 = new double[NX_ATMO];
	double* dvcf3 = new double[NX_ATMO];
	double* dvcf4 = new double[NX_ATMO];
	double* dvcf5 = new double[NX_ATMO];
	double* dvcf6 = new double[NX_ATMO];
	double* dvcf7 = new double[NX_ATMO];
	tvOK = getglobaldata_double(ncid, "bare_ground", dvcf1);
	tvOK = getglobaldata_double(ncid, "Broadleaf", dvcf2);
	tvOK = getglobaldata_double(ncid, "Needleleaf", dvcf3);
	tvOK = getglobaldata_double(ncid, "Evergreen", dvcf4);
	tvOK = getglobaldata_double(ncid, "Decidous", dvcf5);
	tvOK = getglobaldata_double(ncid, "Tree_cover", dvcf6);
	tvOK = getglobaldata_double(ncid, "Herb", dvcf7);
	nc_close(ncid);

	// Print out veg cover
	FILE *ofp2;
	ofp2 = fopen("ece_vegcover_T159_con2.txt", "w");

	// storing soil code in Coord c
	gridlist.firstobj();

	int mask_index = 0;

	for (ix = 0; ix<NX_ATMO; ix++) {
		if (coup_lsmask[0][ix] == 1.0) {

			EceCoord& c = gridlist.getobj();
			int cd = (int)dcoup_soilcd[ix];
			if (cd <= 0 || cd > 7) {
				dprintf("Bad soil code! \n");
				return false;
			}
			c.soilcode = cd;

			c.cvl = dcvl[ix];
			c.cvh = dcvh[ix];
			c.tvl = (int)dtvl[ix];
			c.tvh = (int)dtvh[ix];

			fprintf(ofp2, "%8.5f\t %8.5f\t %8.5f\t %d\t %d\t %d\t %8.5f\t %8.5f\t %d\t %d\t %8.5f\t %8.5f\t %8.5f\t %8.5f\t %8.5f\t %8.5f\t %8.5f \n", 
				c.lon, c.lat, c.area, c.soilcode, c.IFSindex, mask_index, c.cvl, c.cvh, c.tvl, c.tvh,
				dvcf1[ix], dvcf2[ix], dvcf3[ix], dvcf4[ix], dvcf5[ix], dvcf6[ix], dvcf7[ix]);

			gridlist.nextobj();
			mask_index++;
		}
	}

	delete[] dcoup_soilcd;
	delete[] dcvl;
	delete[] dcvh;
	delete[] dtvl;
	delete[] dtvh;
	delete[] dvcf1;
	delete[] dvcf2;
	delete[] dvcf3;
	delete[] dvcf4;
	delete[] dvcf5;
	delete[] dvcf6;
	delete[] dvcf7;

	// clode vcf file
	fclose(ofp2);

	dprintf("Soil and static vegetation and cover data read in and printed\n");
	return true;
}

*/

/*
// Used for creating files to be used with cdo, remapcon
bool print_grid() {

	// Print the grid if needed for cdo operations!

	const long int maxpoints = 100000; // For testing. Otherwise set to 100000

	// Print out gridlist?
	FILE *ofp;
	//ofp = fopen("grids_T255.txt", "w");
	ofp = fopen("grids_T159_all.txt", "w");

	//fprintf(ofp, "%s\n%s\n%s\n", "gridtype = cell", "gridsize = 25799", "nvertex=4"); // grids_T255
	//fprintf(ofp, "%s\n%s\n%s\n", "gridtype = cell", "gridsize = 88838", "nvertex=4"); // grids_T255_all
	//fprintf(ofp, "%s\n%s\n%s\n", "gridtype = cell", "gridsize = 10407", "nvertex=4"); // grids_T159
	fprintf(ofp, "%s\n%s\n%s\n", "gridtype = cell", "gridsize = 35718", "nvertex=4"); // grids_T159_all
	
	// storing soil code in Coord c
	gridlist.firstobj();

	fprintf(ofp, "%s\n","xvals =");

	int test = 0;
	while(gridlist.isobj && test < maxpoints) {
		EceCoord& c=gridlist.getobj();

		fprintf(ofp, "%8.6f\n",c.lon);

		gridlist.nextobj();
		test++;
	}

	test = 0;
	fprintf(ofp, "%s\n","xbounds =");

	gridlist.firstobj();

	while(gridlist.isobj && test < maxpoints) {
		EceCoord& c=gridlist.getobj();

		// Uncomment when needed!
		fprintf(ofp, "%8.6f %8.6f %8.6f %8.6f\n",c.cornerlo[0],c.cornerlo[1],c.cornerlo[2],c.cornerlo[3]);

		gridlist.nextobj();
		test++;
	}


	test = 0;
	gridlist.firstobj();

	fprintf(ofp, "%s\n","yvals =");

	while(gridlist.isobj && test < maxpoints) {
		EceCoord& c=gridlist.getobj();

		fprintf(ofp, "%8.6f\n",c.lat);

		gridlist.nextobj();
		test++;
	}

	test = 0;
	fprintf(ofp, "%s\n","ybounds =");

	gridlist.firstobj();

	while(gridlist.isobj && test < maxpoints) {
		EceCoord& c=gridlist.getobj();

		// Uncomment when needed!
		fprintf(ofp, "%8.6f %8.6f %8.6f %8.6f\n",c.cornerla[0],c.cornerla[1],c.cornerla[2],c.cornerla[3]);

		gridlist.nextobj();
		test++;
	}

	// print gridlist?
	fclose(ofp);

	dprintf("Grid printed \n");
	return true;
}

*/


/*
static void readco2() {

    // Reads in atmospheric CO2 concentrations for EC_Earth historical period, 1840-2010
    // from ascii text file with records in format: <year> [CO2] [CH4] [N2O] [CFC11] [CFC12]

    int year,calender_year;

    FILE* in=fopen(file_co2,"rt");
    if (!in) {
		dprintf("readco2 in guessio: could not open CO2 file %s for input\n",file_co2);
		exit(99);
    } else {
		dprintf("readco2 in guessio - CO2 file opened OK\n");
	}
	
	// I deleted the first read the two header lines
	readfor(in,"");
    readfor(in,"");

	for (year=0;year<NYEAR_ECEARTH-100;year++) { 
		// was NYEAR_HIST. Changed because ghg_histo.txt has 171 years of data
		double ch4,n2o,cfc11,cfc12;
		readfor(in,"i,f,f,f,f,f",&calender_year,&co2[year],&ch4,&n2o,&cfc11,&cfc12);
		if (calender_year!=FIRSTHISTYEAR+year) {
			dprintf("readco2: %s, line %d, cal_year %d, FIRSTHISTYEAR %d - incorrect year specified",file_co2,year,calender_year, FIRSTHISTYEAR);
			exit(99);
		}
	}

	// TEST?
	dprintf("readco2 in guessio : CO2 in 1840: %12.6f \n",co2[0]);
	dprintf("readco2 in guessio : CO2 in 2010: %12.6f \n",co2[NYEAR_ECEARTH-101]);

	fclose(in);
}
*/


/*
// Use unix2dos instead:

static void readandwrite_lu() {

	long int line;

	FILE* in = fopen("lu_1850_2015_gridece.txt", "rt");
	if (!in) {
		dprintf("could not open LU file %s for input\n");
		exit(99);
	}

	FILE* out = fopen("lu_1850_2015_gridece_win.txt", "wt");
	if (!out) {
		dprintf("could not open LU file %s for output\n");
		exit(99);
	}

	bool eof = false;
	xtring headervals[9];

	// header line
	readfor(in, "9a", headervals); // URBAN, PASTURE, CROPLAND, NATURAL, PEATLAND, BARREN;
	// fprintf(out, "9%8s\n", headervals[0], headervals[1], headervals[2], headervals[3], headervals[4], headervals[5], headervals[6], headervals[7], headervals[8]);
	for (int i = 0; i<9; i++)
		fprintf(out, "%8s ", (char *)headervals[i]);
	fprintf(out, "\n");

	double dlon, dlat;
	long double luvals[6];
	int year; 

	line = 0;

	while (!eof) {

		// Read next record in file
		eof = !readfor(in, "f,f,i,6f.14", &dlon, &dlat, &year, luvals);

		if (!eof && !(dlon == 0.0 && dlat == 0.0)) { // ignore blank lines at end (if any)
		
			fprintf(out, "%8.6f %8.6f %d %16.14f %16.14f %16.14f %16.14f %16.14f %16.14f\n",dlon, dlat, year, luvals[0], luvals[1], luvals[2], luvals[3], luvals[4], luvals[5]);

			line++;
		}
	}

	fclose(in);
	fclose(out);

	dprintf("Finished LU!!!\n");
}


static void readandwrite_crop() {

	long int line;

	FILE* in = fopen("crop_1850_2015_gridece.txt", "rt");
	if (!in) {
		dprintf("could not open crop file %s for input\n");
		exit(99);
	}

	FILE* out = fopen("crop_1850_2015_gridece_win.txt", "wt");
	if (!out) {
		dprintf("could not open crop file %s for output\n");
		exit(99);
	}

	bool eof = false;
	xtring headervals[8];

	// header line
	readfor(in, "8a", headervals); // CC3ann  CC3per  CC3nfx  CC4ann  CC4per
	for (int i = 0; i<8; i++)
		fprintf(out, "%8s ", (char *)headervals[i]);
	fprintf(out, "\n");

	double dlon, dlat;
	long double luvals[5];
	int year;

	line = 0;

	while (!eof) {

		// Read next record in file
		eof = !readfor(in, "f,f,i,5f.6", &dlon, &dlat, &year, luvals);

		if (!eof && !(dlon == 0.0 && dlat == 0.0)) { // ignore blank lines at end (if any)

			fprintf(out, "%8.6f %8.6f %d %8.6f %8.6f %8.6f %8.6f %8.6f\n", dlon, dlat, year, luvals[0], luvals[1], luvals[2], luvals[3], luvals[4]);

			line++;
		}
	}

	fclose(in);
	fclose(out);

	dprintf("Finished CROP!!!\n");
}


static void readandwrite_nfert() {

	long int line;

	FILE* in = fopen("nfert_1850_2015_gridece.txt", "rt");
	if (!in) {
		dprintf("could not open nfert file %s for input\n");
		exit(99);
	}

	FILE* out = fopen("nfert_1850_2015_gridece_win.txt", "wt");
	if (!out) {
		dprintf("could not open nfert file %s for output\n");
		exit(99);
	}

	bool eof = false;
	xtring headervals[8];

	// header line
	readfor(in, "8a", headervals); // CC3ann  CC3per  CC3nfx  CC4ann  CC4per
	for (int i = 0; i<8; i++)
		fprintf(out, "%8s ", (char *)headervals[i]);
	fprintf(out, "\n");

	double dlon, dlat;
	long double luvals[5];
	int year;

	line = 0;

	while (!eof) {

		// Read next record in file
		eof = !readfor(in, "f,f,i,5f.6", &dlon, &dlat, &year, luvals);

		if (!eof && !(dlon == 0.0 && dlat == 0.0)) { // ignore blank lines at end (if any)

			fprintf(out, "%8.6f %8.6f %d %8.6f %8.6f %8.6f %8.6f %8.6f\n", dlon, dlat, year, luvals[0], luvals[1], luvals[2], luvals[3], luvals[4]);

			line++;
		}
	}

	fclose(in);
	fclose(out);

	dp rintf("Finished N FERT!!!\n");
}

*/

//////////////////////////////////////////////////////////////////////////////////////
// READ ALL INPUT DATA

void read_input_data(bool& error_flag) { 
	
	// ecev3 - MUCH simplified

	bool gridOK = read_grid_info(false);
	if (!gridOK) {
		dprintf("LPJ-GUESS: Problem with atmosphere model's mask and grid information.  \n");
		error_flag = true;
		//exit(99);
	}

	// Read IFS timesteps
	bool rcfileOK = read_ifs_timesteps();
	if (!rcfileOK) {
		dprintf("Problem with EC-Earth .rc file. Quitting. \n");
		error_flag = true;
		//exit(99);
	}

	// LU reformatting (but best to use unix2dos):
	// readandwrite_lu();
	// readandwrite_crop();
	// readandwrite_nfert();

	// Read CO2 from CMIP6
	readco2_cmip6(error_flag);
	// readco2(); // old, CMIP5 method

	// Read soil and static vegetation information
	bool soilOK = read_soil_veg_info(false);
	if (!soilOK)  {
		dprintf("Problem with EC-Earth Soil info. \n");
		error_flag = true;
	}
	// Print SCRIP grid description file?
	// print_grid();
 
}


//////////////////////////////////////////////////////////////////////////////////////
// CALLS TO OASIS

void call_oasis_get(int timestep){
  
	// Obtaining/"get" atmosphere data from OASIS

	// LAI check when timestep < 20
	if (timestep < 20) {
		dprintf("call_oasis_get on timestep %i - LPJ-GUESS: OASIS communicating\n",timestep);
	}

	// Exchange fields
	// ecev3 - use activateTM5coupling, set from .rc file. Added coup_vegl_type.
	OasisCoupler::couple_get(timestep*TIMESTEP,NX_ATMO,NY_ATMO, activateTM5coupling, coup_temper, 
				 coup_precip, coup_snowc, coup_snowd,coup_st1l, coup_st2l, coup_st3l, 
				 coup_st4l, coup_sm1l, coup_sm2l, coup_sm3l, coup_sm4l, coup_sw_radiat, coup_lw_radiat, 
				 coup_co2_field);

	if (timestep < 20) {
		dprintf("call_oasis_get - LPJ-GUESS: OASIS communicated!\n");
	}
}



void call_oasis_put(int timestep){
  
	// Send/"put" vegetation data to OASIS
  
	int ix;

	// LAI check when timestep < 5
	if (timestep < 5) {

		dprintf("call_oasis_put on timestep %i - LPJ-GUESS: OASIS communicating\n",timestep);

		float max_llai = 0.0;
		float max_hlai = 0.0;
		float mean_llai = 0.0;
		float mean_hlai = 0.0;
		float mean_cflux = 0.0;
		float mean_npp = 0.0;

		for(ix=0;ix<NX_ATMO;ix++){
			if(coup_lsmask[ix][0]){
				
				if (coup_lowlai[ix] > max_llai) max_llai = (float)coup_lowlai[ix]; 
				if (coup_highlai[ix] > max_hlai) max_hlai = (float)coup_highlai[ix];
			
				mean_hlai += (float)coup_highlai[ix]/MAXGRID;	
				mean_llai += (float)coup_lowlai[ix]/MAXGRID;
				mean_cflux += (float)coup_dcflux[ix]/MAXGRID; 
				mean_npp += (float)coup_dnpp[ix]/MAXGRID; 
			}
		}

		dprintf("call_oasis_put - Values on timestep %i - max LOW, max HIGH, mean LOW, mean HIGH LAI, mean CFLUX, mean NPP are %12.6f, %12.6f, %12.6f, %12.6f, %12.6f, %12.6f\n",
			timestep,max_llai,max_hlai,mean_llai,mean_hlai,mean_cflux, mean_npp);
	
	}

	// ecev3 - use activateTM5coupling
	OasisCoupler::couple_put(timestep*TIMESTEP,NX_ATMO,NY_ATMO,activateTM5coupling,coup_lowlai, coup_highlai, coup_lpjg_typeh, 
		coup_lpjg_frach, coup_lpjg_typel, coup_lpjg_fracl, coup_dcflux, coup_dnpp);

	// dprintf("call_oasis_put - LPJ-GUESS: OASIS communicated!\n");
        
}


void call_coupler_get(int timestep,double temp[MAXGRID],double prec[MAXGRID],
	double vegl[MAXGRID], int vegl_type[MAXGRID], double vegh[MAXGRID], double snowc[MAXGRID], double snowd[MAXGRID],
	double stl[MAXGRID][NHTESSELSOILLAYERS], double sml[MAXGRID][NHTESSELSOILLAYERS], 
	double swrad[MAXGRID], double lwrad[MAXGRID],double co2_field[MAXGRID]) {
 
	int g;
	int ix;
	int iy;
	
	const double TEMP_KTOC=-273.15;

	// This is the point at which to get the fields from OASIS
	call_oasis_get(timestep);

	// timestep = timestep number from 0 for start date
	// temp = temperature (degC instantaneous)
	// prec = precipitation (mm since last time step)
	// swrad = net downward shortwave radiation (W/m2 instantaneous)
	// lwrad = net downward longwave radiation (W/m2 instantaneous)
	// co2 = atmospheric CO2 concentration (ppmv)
	// temp, prec and rad are 1-dimensional arrays containing data for all grid cells

	// reduce complete grid arrays received from OASIS to one-dimensional arrays within GUESS
	g=0;

	// Doubles to hold global averages for this timestep. Use in debugging!
	double temp_avg = 0.0;
	double prec_avg = 0.0;
	double swrad_avg = 0.0;
	double lwrad_avg = 0.0;

	double vegl_avg = 0.0;
	double vegh_avg = 0.0;

	double snowc_avg = 0.0;
	double snowd_avg = 0.0;

	double st1l_avg = 0.0;
	double st2l_avg = 0.0;
	double st3l_avg = 0.0;
	double st4l_avg = 0.0;

	double sm1l_avg = 0.0;
	double sm2l_avg = 0.0;
	double sm3l_avg = 0.0;
	double sm4l_avg = 0.0;

	double co2_avg = 0.0;

	for(iy=0;iy<NY_ATMO;iy++){
		for(ix=0;ix<NX_ATMO;ix++){

			// global averages
			temp_avg += coup_temper[ix]/(double)NX_ATMO;
			prec_avg += coup_precip[ix]/(double)NX_ATMO;	

			swrad_avg += coup_sw_radiat[ix] / (double)NX_ATMO;
			lwrad_avg += coup_lw_radiat[ix] / (double)NX_ATMO;

			snowc_avg += coup_snowc[ix]/(double)NX_ATMO;
			snowd_avg += coup_snowd[ix]/(double)NX_ATMO;

			st1l_avg += coup_st1l[ix]/(double)NX_ATMO;
			st2l_avg += coup_st2l[ix]/(double)NX_ATMO;
			st3l_avg += coup_st3l[ix]/(double)NX_ATMO;
			st4l_avg += coup_st4l[ix]/(double)NX_ATMO;

			sm1l_avg += coup_sm1l[ix]/(double)NX_ATMO;
			sm2l_avg += coup_sm2l[ix]/(double)NX_ATMO;
			sm3l_avg += coup_sm3l[ix]/(double)NX_ATMO;
			sm4l_avg += coup_sm4l[ix]/(double)NX_ATMO;

			co2_avg += coup_co2_field[ix]/(double)NX_ATMO;

			if(coup_lsmask[ix][iy]){
	      
				// Unit conversions
				// temperature (K) --> (oC)
				// surface solar radiation (W m-2 s) --> (W m-2)
				// precipitation (m) --> (mm)
				temp[g]=coup_temper[ix]+TEMP_KTOC;
     
				// Convert from kg m-2 s-1 to mm
				prec[g]=coup_precip[ix]*1.e3*TIMESTEP/1000; 
				if (prec[g] < 0.000001) prec[g] = 0.0; // For negative or very low values

				snowc[g] = coup_snowc[ix]; // kg m-2
				snowd[g] = coup_snowd[ix]; // kg m-3

				stl[g][0] = coup_st1l[ix]+TEMP_KTOC; // all K to degC
				stl[g][1] = coup_st2l[ix]+TEMP_KTOC;
				stl[g][2] = coup_st3l[ix]+TEMP_KTOC;
				stl[g][3] = coup_st4l[ix]+TEMP_KTOC;
	
				sml[g][0] = coup_sm1l[ix]; // all m3 m-3
				sml[g][1] = coup_sm2l[ix];
				sml[g][2] = coup_sm3l[ix];
				sml[g][3] = coup_sm4l[ix];

				// ECE sends W m-2
				swrad[g] = coup_sw_radiat[ix];
				lwrad[g] = coup_lw_radiat[ix];
				if (swrad[g] < 0.000001) swrad[g] = 0.0; // For negative or very low values
	     
				co2_field[g] = coup_co2_field[ix];

				g++;
			}
		}
	}


	if (timestep < 20) {

		// Print checksums
		dprintf("ECEFRAMEWORK: (checksum TEMP) %12.6f\n",temp_avg);
		dprintf("ECEFRAMEWORK: (checksum PREC) %12.6f\n",prec_avg);
        dprintf("ECEFRAMEWORK: (checksum SNOWC) %12.6f\n",snowc_avg);
        dprintf("ECEFRAMEWORK: (checksum SNOWD) %12.6f\n",snowd_avg);
        dprintf("ECEFRAMEWORK: (checksum ST1L) %12.6f\n",st1l_avg);
        dprintf("ECEFRAMEWORK: (checksum ST2L) %12.6f\n",st2l_avg);
        dprintf("ECEFRAMEWORK: (checksum ST3L) %12.6f\n",st3l_avg);
        dprintf("ECEFRAMEWORK: (checksum ST4L) %12.6f\n",st4l_avg);
        dprintf("ECEFRAMEWORK: (checksum SM1L) %12.6f\n",sm1l_avg);
        dprintf("ECEFRAMEWORK: (checksum SM2L) %12.6f\n",sm2l_avg);
        dprintf("ECEFRAMEWORK: (checksum SM3L) %12.6f\n",sm3l_avg);
        dprintf("ECEFRAMEWORK: (checksum SM4L) %12.6f\n",sm4l_avg);
		dprintf("ECEFRAMEWORK: (checksum SW RAD) %12.6f\n", swrad_avg);
		dprintf("ECEFRAMEWORK: (checksum LR RAD) %12.6f\n", lwrad_avg);
		dprintf("ECEFRAMEWORK: (checksum CO2) %12.6f\n",co2_avg);
	}
   
	return;
}



void call_coupler_put(int timestep,double lailow[MAXGRID], double laihigh[MAXGRID],	// LAI
	int lpjg_typeh[MAXGRID], double lpjg_frach[MAXGRID],							// VEG H
	int lpjg_typel[MAXGRID], double lpjg_fracl[MAXGRID],							// VEG L
	double dcflux[MAXGRID], double dnpp[MAXGRID]) {									// C FLUX
 
	int g;
	int ix;
	int iy;
	
	// Expand one-dimensional arrays within GUESS to complete grid for sending to OASIS
	g=0;
	for(iy=0;iy<NY_ATMO;iy++) {
		for(ix=0;ix<NX_ATMO;ix++) {

			// Initialise
			coup_lowlai[ix]=0.;
			coup_highlai[ix]=0.;
	        coup_lpjg_typeh[ix]=0.;
		    coup_lpjg_frach[ix]=0.;
	        coup_lpjg_typel[ix]=0.;
		    coup_lpjg_fracl[ix]=0.;
			coup_dcflux[ix]=0.;
			coup_dnpp[ix]=0.;

			// Fill land cells
			if(coup_lsmask[ix][iy]){
				// LAI
				coup_lowlai[ix]=lailow[g];
				coup_highlai[ix]=laihigh[g];
				// HIGH VEG
				coup_lpjg_typeh[ix]=(double)lpjg_typeh[g]; // convert int to double
				coup_lpjg_frach[ix]=lpjg_frach[g];
				// LOW VEG
				coup_lpjg_typel[ix]=(double)lpjg_typel[g]; // convert int to double
				coup_lpjg_fracl[ix]=lpjg_fracl[g];
				// C FLUX
				coup_dcflux[ix]=dcflux[g];
				coup_dnpp[ix]=dnpp[g];
				g++;
			}
		}
	}

	// This is the point at which to send the fields to OASIS 
	call_oasis_put(timestep);

	return;
}


//////////////////////////////////////////////////////////////////////////////////////
// HELPER FUNCTIONS

// IFS layers have thicknesses: 7, 21, 72 and 189cm, respectively. 
// Layer centres are given (in IFS output files) as 3, 17, 64 and 177 cm.

// Interpolate linearly between soil temperatures in the 2nd and 3rd IFS soil layers to get 
// T25. Layer 2 centre: 17cm. Layer 3 centre: 64cm. 
float ifs_to_lpjg_soiltemp(double temp_soil_2, double temp_soil_3) {

	float slope = ((float)temp_soil_3 - (float)temp_soil_2) / (64.0 - 17.0);
	float t25 = slope * (25.0 - 17.0) + (float)temp_soil_2;

	return t25;
}


void ifs_to_lpjg_soilwater(double sm1l, double sm2l, double sm3l, double sm4l, float& lpjg_sw_upper, float& lpjg_sw_lower) {

	// Calculate weighted averages of IFS layer soil water contents (m3 m-3) to get values for
	// LPJ-GUESS upper and lower soil layers
	lpjg_sw_upper = (float) (7.0 * sm1l + 21.0 * sm2l + 22.0 * sm3l) / 50.0;
    lpjg_sw_lower = (float) (50.0 * sm3l + 50.0 * sm4l) / 100.0;

	return;
} 


// ecev3 
// Run LPJ-GUESS for ONE day, for ALL gridcells
// called from runlpjguess
void runlpjguess_today(
	
	// Spinup?
	bool islpjgspinup,
	// Date and time information
	int timestep,int isfinal,int yearc,int sim_yr,int monthc,int dayc,int time,
	// Fields received FROM IFS (but only the master process has the data when this routine is called) 
	double temp[MAXGRID],double prec[MAXGRID],double swrad[MAXGRID], double lwrad[MAXGRID], 
	double snowc[MAXGRID], double snowd[MAXGRID],
	double stl[MAXGRID][NHTESSELSOILLAYERS], double sml[MAXGRID][NHTESSELSOILLAYERS],
	// Fields sent TO IFS 
	double lailow[MAXGRID], double laihigh[MAXGRID], 
	int lpjg_typeh[MAXGRID], double lpjg_frach[MAXGRID],
	int lpjg_typel[MAXGRID], double lpjg_fracl[MAXGRID], // ecev3 - new
	// For TM5 communication
	double co2_ppm[MAXGRID], double dcflux[MAXGRID], double dnpp[MAXGRID]) {

	/*
	Main loop through all the land cells TODAY, and calls LPJG for each cell
	Climate inputs from runlpjguess.
	Returns LAI, carbon fluxes etc. for ALL cells today
	*/

	int ifs_soilcd = -1;
	int ifs_index = -1;
	int ndep_index = -1;

	int g=0;
	int hourc=0;;
	float alat=0.0;
	float alon=0.0;

	// Inputs
	float temp_2m=0.0;
	float temp_soil=0.0;
	float soilw_surf=0.0;
	float soilw_deep=0.0;
	float swrad_net_inst = 0.0;
	float lwrad_net_inst = 0.0;
	float temp_soil_1=0.0;
	float temp_soil_2=0.0;
	float temp_soil_3=0.0;
	float temp_soil_4=0.0;
	float precip=0.0;
	float co2=0.0;
	float vegl_cell=0.0;
	int vegl_type_cell=0;
	float vegh_cell=0.0;
	int vegh_type_cell=0;

	// Outputs
	float cfluxtoday=0.0;
	float npptoday=0.0;
	float laiphen_high=0.0;
	float laiphen_low=0.0;
	float ifsvegfraclow=0.0; 
	int ifsvegtypelow=0;
	float ifsvegfrachigh=0.0; 
	int ifsvegtypehigh=0;

	// MPI variables
	#ifdef HAVE_MPI
		MPI_Status status;
	#endif

	int rank; 
	int num_procs;


	// ecev3 - will need when we go from DAILY timesteps to using the diurnal code
	hourc=time/3600; // 0, 6, 12 or 18
	
	// Determine the gridcells to simulate on this processor
	rank = get_rank_specific(localcomm);
	num_procs = get_num_local_processes(localcomm);
	

	// Report the cells this node will simulate
	if (timestep == 0) dprintf("Simulating %i cells on node %i: starting at %i in steps of %i \n", (int)((MAXGRID-rank)/num_procs), rank, rank, num_procs);

	// Initialise arrays to be returned to IFS and TM5
	// Shuffling of gridcells now done along the latidues such 
	// that every proc gets hi and lo lat gridcells to ensure equal memory
	// usage. LN 12/2017
	for (g = rank; g<MAXGRID; g+=num_procs) {
	  
		dcflux[g]=0.0;
		dnpp[g]=0.0;
		lailow[g]=0.0;
		laihigh[g]=0.0;
		lpjg_typeh[g]=0;
		lpjg_frach[g]=0.0;
		lpjg_typel[g]=0;
		lpjg_fracl[g]=0.0;

		// ECEtest - outcomment these lines to return nonzero values to IFS
		/*
		lailow[g]=2;		// Short grass
		laihigh[g]=5;		// EN trees
		lpjg_typeh[g]=3;	// EN trees
		lpjg_frach[g]=0.5;	
		lpjg_typel[g]=2;	// Short grass
		lpjg_fracl[g]=0.5;
		dcflux[g]=0.001;	// kgC/m2 
		dnpp[g]=0.001;		// kgC/m2 
		*/
	}

	// ECEtest - outcomment this line if you want to test the exchange of fields only. 
    // return; 


	// Transfer soil temp and moisture information to 1D arrays before MPI broadcasting
	double stl1[MAXGRID],stl2[MAXGRID],stl3[MAXGRID],stl4[MAXGRID];
	double sml1[MAXGRID],sml2[MAXGRID],sml3[MAXGRID],sml4[MAXGRID];

	if (rank==0) {
		// Only process 0 has the data from OASIS/IFS
		for (int ii = 0; ii < MAXGRID; ii++) {
			stl1[ii] = stl[ii][0];
			stl2[ii] = stl[ii][1];
			stl3[ii] = stl[ii][2];
			stl4[ii] = stl[ii][3];

			sml1[ii] = sml[ii][0];
			sml2[ii] = sml[ii][1];
			sml3[ii] = sml[ii][2];
			sml4[ii] = sml[ii][3];
		}
	}

	// FULL simulation

	// Report the cells this node will simulate
	if (timestep<2) dprintf("2. Simulating %i cells on node %i: starting at %i in steps of %i \n", (int)((MAXGRID-rank)/num_procs), rank, rank, num_procs);
                
	#ifdef HAVE_MPI

    // ECEtest - outcomment this line if you want to test the exchange of fields only.
    // return;
 
	if (!islpjgspinup) {

		if (timestep<2) dprintf("Broadcasting on node %i\n", rank);

		// Broadcast forcing information to all LOCAL processes
		MPI_Bcast(temp,MAXGRID,MPI_DOUBLE,0,localcomm);

		if (timestep<2) dprintf("Broadcasted temp on node %i\n", rank);

		MPI_Bcast(prec   ,MAXGRID,MPI_DOUBLE,0,localcomm);
		MPI_Bcast(swrad  ,MAXGRID,MPI_DOUBLE,0,localcomm);
		MPI_Bcast(lwrad  ,MAXGRID,MPI_DOUBLE,0,localcomm);
		MPI_Bcast(stl1   ,MAXGRID,MPI_DOUBLE,0,localcomm);
		MPI_Bcast(stl2   ,MAXGRID,MPI_DOUBLE,0,localcomm);
		MPI_Bcast(stl3   ,MAXGRID,MPI_DOUBLE,0,localcomm);
		MPI_Bcast(stl4   ,MAXGRID,MPI_DOUBLE,0,localcomm);
		MPI_Bcast(sml1   ,MAXGRID,MPI_DOUBLE,0,localcomm);
		MPI_Bcast(sml2   ,MAXGRID,MPI_DOUBLE,0,localcomm);
		MPI_Bcast(sml3   ,MAXGRID,MPI_DOUBLE,0,localcomm);
		MPI_Bcast(sml4   ,MAXGRID,MPI_DOUBLE,0,localcomm);
		MPI_Bcast(co2_ppm,MAXGRID,MPI_DOUBLE,0,localcomm);

		if (timestep<2) dprintf("Finished broadcasting on node %i\n", rank);

	}

	#endif

	// Some typical T255 cells, gg = 
	// 305			- Tundra
	// 976			- Larch (4) and evergreen shrubs (16)
	// 1707			- Tundra 
	// 1335			- Canadian bog 
	// 2326			- Russian bog
	// 2363			- Spassky Pad (Larch 62.8150°N, 129.8370°E; 220 m)
	// 2708			- DNF (Larch)
	// 5526			- Central Germany
	// 6654			- N Mongolia (short grass, TVL = 2, TVH = 0)
	// 9311			- Las Vegas (99% Semidesert in IFS)
	// 10423		- Semi desert
	// 13097		- Bangladesh (98% Irrigated crops in IFS)
	// 13203		- Sahara point
	// 14354		- Arabian peninsula (semidesert, TVL = 11, TVH = 0)
	// 16100		- S America, vegcodes (crop 1, mixed 19)
	// 16446		- Interrupted forest
	// 18211		- Amazon point
	// 19248		- Africa (codes 7 & 5)
	// 19437		- Another S America point (7, 19)
	// 21107		- Africa, vegcodes (tall grass 7, mixed 19)
	// 20579-20620	- transect across S. America at 20 deg south.
	// 20504-20534	- transect across S. Africa at 20 deg south.

	// Some typical T159 cells, gg = 
	// 207			- Tundra, on Arctic coast
	// 5526			- E India


	// An array of cells for testing purposes
	const int numtestcells = 20;
	int testcells[numtestcells] = { 305, 976, 1707, 1335, 2326, 2363, 2708, 5526, 6654, 9311, 10423, 13097, 13203, 14354, 16100, 16446, 18211, 19248, 19437, 21107 };

	// int testcells[2] = {207, 5526}; // T159

	// Loop through the selected cells
	//CLN ORIGINAL :for (int gg = firstcell; gg<lastcell; gg++) {
	for (int gg = rank; gg<MAXGRID; gg+=num_procs) {
	// for (int ii = 0; ii<numtestcells; ii++) {
	// for (int ii = 14; ii<15; ii++) {
		// int gg = testcells[ii]; // test cells only. Remove before checking in code.

		if (timestep==0) 
			dprintf("Simulating %i cells on node %i: starting at %i in steps of %i \n", (int)((MAXGRID-rank)/num_procs), rank, rank, num_procs);
		
		// PREPARE INPUT DATA TO LPJ-GUESS

		// Reset for each cell, every timestep
		cfluxtoday=0.0;
		npptoday=0.0;
		laiphen_high=0.0;
		laiphen_low=0.0;

		// Get reference to this grid cell
		EceCoord& c=gridlist[gg];
		alat=(float)c.lat;
		alon=(float)c.lon;
		ifs_soilcd = c.soilcode; // soil code (1-8)
		ifs_index = c.IFSindex;
		ndep_index = c.ndep_index;

		temp_2m=(float)temp[gg]; // K
		temp_soil_1=(float)stl1[gg];
		temp_soil_2=(float)stl2[gg];
		temp_soil_3=(float)stl3[gg];
		temp_soil_4=(float)stl4[gg];

		// Interpolate soil temp to 25cm. Already in deg C
		temp_soil = ifs_to_lpjg_soiltemp(temp_soil_2,temp_soil_3); // deg C

		// Now send precip to LPJ-GUESS
		precip = (float)prec[gg];

		// CO2 received from TM5, or file.
		co2 = co2_ppm[gg];

		// Use IFS soil moisture to update soilw_surf and soilw_deep
		ifs_to_lpjg_soilwater(sml1[gg], sml2[gg], sml3[gg], sml4[gg], soilw_surf, soilw_deep);

		// SW radiation
		swrad_net_inst = (float)swrad[gg];

		// LW radiation
		lwrad_net_inst = (float)lwrad[gg];

		// Veg. cover and type from netcdf input files
		vegl_cell = c.cvl;
		vegl_type_cell = c.tvl;
		vegh_cell = c.cvh;
		vegh_type_cell = c.tvh;

		
		// ecev3 - TEST DATA for one tropical gridcell. MUST outcomment!!!!
		/*
		soilw_surf = 0.6;
		soilw_deep = 0.6;
		if (monthc < 0 || monthc > 13) {
			precip = 1.0;
			temp_soil = 3.0;
			swrad_net_inst = 100.0;
			temp_2m = 22.0;
		}
		else { // Amazon point
			precip = 6.0;
			temp_soil = 26.0;
			swrad_net_inst = 170.0;
			temp_2m = 24.0;
		}
		*/


		/*
		// ecev3 - TEST DATA for one German gridcell. MUST outcomment!!!!
		soilw_surf = 0.33;
		soilw_deep = 0.33;
		if (monthc < 4 || monthc > 8) {
		precip = 2.4;
		temp_soil = 5.0;
		swrad_net_inst = 90.0;
		temp_2m = 5.0;
		}
		else {
		precip = 2.4;
		temp_soil = 15.0;
		swrad_net_inst = 150.0;
		temp_2m = 15.0;
		}
		*/

		/*
		// ecev3 - TEST DATA for one German gridcell. MUST outcomment!!!!
		soilw_surf = 0.0;
		soilw_deep = 0.0;
		if (monthc < 4 || monthc > 8) {
		precip = 0.1;
		temp_soil = -1.0;
		swrad_net_inst = 30.0;
		temp_2m = -10.0;
		}
		else {
		precip = 0.5;
		temp_soil = 1.0;
		swrad_net_inst = 150.0;
		temp_2m = 1.0;
		}
		*/

		// Call LPJ-GUESS for THIS cell TODAY !!!!! 
		// SEND: Dates, coordinates, climate and CO2 for this cell, today
		// RECEIVE: LAI for high and low Stands in the Gridcell, fraction and dominant type of high vegetation, as well as carbon fluxes

		// HERE is where we call the LPJG trunk framework!!!!!!! 
		// Something VERY SIMILAR to a direct call to guess_coupled - see RCA-GUESS for inspiration
        if (false && timestep < 2) {
			// ECEtest
			dprintf("GUESS: before guess_coupled on timestep %i for cell %i with lat %f, lon %f\n",timestep,ifs_index,alat,alon);
			dprintf("GUESS: temp %f, precip %f, CO2 %f\n",temp_2m, precip, co2);
			dprintf("GUESS: swrad_net_inst %f, lwrad_net_inst %f, temp_soil %f, soil code %i\n",swrad_net_inst, lwrad_net_inst, temp_soil, ifs_soilcd);
			dprintf("GUESS: soilw_surf %f, soilw_deep %f, vegl_cell %f, vegl_type_cell %i, vegh_cell %f, vegh_type_cell %i\n",soilw_surf, soilw_deep, 
				vegl_cell, vegl_type_cell, vegh_cell, vegh_type_cell);
		}

		guess_coupled(c.id, rank, num_procs, isfinal, alon, alat, ifs_soilcd, ifs_index, ndep_index, 
			yearc, sim_yr, monthc, dayc, hourc, 
			temp_2m, precip, swrad_net_inst, lwrad_net_inst, co2, temp_soil, soilw_surf, soilw_deep, // Arg19
			vegl_cell, vegl_type_cell, // Inputs
			// Updated fields follow...
			cfluxtoday, npptoday,
			laiphen_high, laiphen_low, 
			ifsvegfrachigh, ifsvegtypehigh, ifsvegfraclow, ifsvegtypelow, 
			islpjgspinup, fixedNDepafter, fixedLUafter);


		// Store lailow and laihigh 
		lailow[gg] = laiphen_low;
		laihigh[gg] = laiphen_high;

		if (ifsvegtypehigh == 0)
			laihigh[gg] = 0.0;
		if (ifsvegtypelow == 0)
			lailow[gg] = 0.0;


		// TYPE is actually only updated once a year
		lpjg_typeh[gg] = (int)ifsvegtypehigh;
		lpjg_frach[gg] = (double)ifsvegfrachigh;
		lpjg_typel[gg] = (int)ifsvegtypelow;
		lpjg_fracl[gg] = (double)ifsvegfraclow;

		// Check for errors
		if (ifsvegtypehigh < 0 || ifsvegtypehigh > 19 || ifsvegfrachigh < 0.0 || ifsvegfrachigh > 1.0 ||
			ifsvegtypelow < 0 || ifsvegtypelow > 19 || ifsvegfraclow < 0.0 || ifsvegfraclow > 1.0) { 
			dprintf("GUESS: Invalid VEG type or fraction after LPJG calls at timestep %i for cell %i \n",timestep, ifs_index);
			dprintf("GUESS: laiphen_high, VEGH type, and VEGH fraction: %f %i %f\n",laiphen_high, ifsvegtypehigh, ifsvegfrachigh);
			dprintf("GUESS: laiphen_low, VEGL type, and VEGL fraction: %f %i %f\n",laiphen_low, ifsvegtypelow, ifsvegfraclow);
		}

		// Scale the C fluxes before sending them to TM5 - done in guess_coupled
		dcflux[gg]=cfluxtoday;
		dnpp[gg]=npptoday;
		// NB! a positive cfluxtoday implies C RELEASE by this gridcell, a negative value implies C UPTAKE

		// ECEtest
		/*
		dprintf("GUESS: LAIH, FRACH, TYPEH at timestep %i for cell %i: %f, %f, %f\n",timestep,ifs_index,laiphen_high,
			lpjg_frach[gg],(double)ifsvegtypehigh);

		dprintf("GUESS: LAIL, FRACL, TYPEL at timestep %i for cell %i: %f, %f, %f\n",timestep,ifs_index,laiphen_low,
			lpjg_fracl[gg],(double)ifsvegtypelow);
		*/

	} 

	// Now get all slave processes to send their information to master process 0
	#ifdef HAVE_MPI

	if (localcomm != MPI_COMM_WORLD && islpjgspinup)
		dprintf("localcomm != MPI_COMM_WORLD \n");
	
	if (localcomm == MPI_COMM_WORLD && islpjgspinup)
		dprintf("localcomm == MPI_COMM_WORLD \n");

	// Synchronise first
	//dprintf("Before MPI_Barrier in runlpjguess_today, after guess_coupled on node %i \n",rank);
	MPI_Barrier(localcomm);
	//dprintf("After MPI_Barrier in runlpjguess_today on node %i \n",rank);

	if (!islpjgspinup) {

		int tag;

		if (num_procs > 1) {

			// Synchronise first
			// MPI_Barrier(localcomm);
        
			// Only do this aggregation if we're running with more than one processor
			double lail[MAXGRID];
			double laih[MAXGRID];
			double frh[MAXGRID];
			int tph[MAXGRID];
			double frl[MAXGRID];
			int tpl[MAXGRID];
			double cfl[MAXGRID];
			double cfl_npp[MAXGRID];

			if (rank == 0) {

				// Master process - receive LPJ-GUESS results from slave processes

				MPI_Status stat;
				for (int pr = 1; pr < num_procs; pr++) {
					// LAIL
					tag = pr*100+1;
					MPI_Recv(lail,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
					// LAIH
					tag = pr*100+2;
					MPI_Recv(laih,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
					// TYPH
					tag = pr*100+3;
					MPI_Recv(tph,MAXGRID,MPI_INT,pr,tag,localcomm,&stat);
					// FRACH
					tag = pr*100+4;
					MPI_Recv(frh,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
					// TYPL
					tag = pr*100+5;
					MPI_Recv(tpl,MAXGRID,MPI_INT,pr,tag,localcomm,&stat);
					// FRACL
					tag = pr*100+6;
					MPI_Recv(frl,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
					// DCFLUX
					tag = pr*100+7;
					MPI_Recv(cfl,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
					// NPP
					tag = pr*100+8;
					MPI_Recv(cfl_npp,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);

					// Report the cells this node will simulate
					dprintf("Master node %i has received the results from slave process %i.\n", rank, pr);

					// T255 - cellspernode goes from 2579 to 2457, but the final node gets more (lastcell is MAXGRID)
					// as it has less to do on account of Antarctica and always finishes most quickly otherwise.
					// The cell limits for this process:
					//CLN int cells = (int)(MAXGRID / (double)(num_procs + 0.5));
					//CLN int cell1 = pr * cells;
					//CLN int cell2 = cell1 + cells;
					//CLN if (pr == num_procs-1)
					//CLN 	cell2 = MAXGRID;
					//CLN 
					//CLN // Now populate the arrays on process 0 ahead of OASIS SEND calls
					//CLN for (int cl = cell1; cl < cell2; cl++) {

					//CLN Longitudinal distribution of HERE SAME translation as above!

					for (int cl = pr; cl<MAXGRID; cl+=num_procs) {
						lailow[cl] = lail[cl];
						laihigh[cl] = laih[cl];
						lpjg_typeh[cl] = tph[cl]; 
						lpjg_frach[cl] = frh[cl];
						lpjg_typel[cl] = tpl[cl]; 
						lpjg_fracl[cl] = frl[cl];
						dcflux[cl] = cfl[cl];
						dnpp[cl] = cfl_npp[cl];
					}	
				} 

			} else {

				// Slave processes - send LPJ-GUESS results to master process

				// LAIL
				tag = rank*100+1;
				MPI_Send(lailow,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
				// LAIH
				tag = rank*100+2;
				MPI_Send(laihigh,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
				// TYPH
				tag = rank*100+3;
				MPI_Send(lpjg_typeh,MAXGRID,MPI_INT,0,tag,localcomm);
				// FRACH
				tag = rank*100+4;
				MPI_Send(lpjg_frach,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
				// TYPL
				tag = rank*100+5;
				MPI_Send(lpjg_typel,MAXGRID,MPI_INT,0,tag,localcomm);
				// FRACH
				tag = rank*100+6;
				MPI_Send(lpjg_fracl,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
				// DCFLUX
				tag = rank*100+7;
				MPI_Send(dcflux,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
				// NPP
				tag = rank*100+8;
				MPI_Send(dnpp,MAXGRID,MPI_DOUBLE,0,tag,localcomm);

				// Report the cells this node will simulate
				if (timestep < 2) dprintf("Node %i has sent the results from its cells to the master.\n", rank);

			}	

			// Synchronise again before moving to the next timestep
			MPI_Barrier(localcomm);
	
		}

	}

	#endif

}

// ecev3
void runlpjguess(bool islpjgspinup, bool isparallel, bool error_flag ) {

	// ecev3 - this is the main time loop.
	// We establish a connection with OASIS, then loop through the days 

	int timestep;
	
	//// Received from IFS
	
	double temp[MAXGRID],prec[MAXGRID],swrad[MAXGRID],lwrad[MAXGRID];
	double vegl[MAXGRID],vegh[MAXGRID],snowc[MAXGRID],snowd[MAXGRID];


	/* HEAP approach 
	double* temp = new double[MAXGRID];
	double* prec = new double[MAXGRID];
	double* swrad = new double[MAXGRID];
	double* lwrad = new double[MAXGRID];
	double* vegl = new double[MAXGRID];
	double* vegh = new double[MAXGRID];
	double* snowc = new double[MAXGRID];
	double* snowd = new double[MAXGRID];
	*/

	double stl[MAXGRID][NHTESSELSOILLAYERS];
	double sml[MAXGRID][NHTESSELSOILLAYERS]; 

	int vegl_type[MAXGRID];

	// Sent to IFS
	double lailow[MAXGRID];
	double laihigh[MAXGRID];
	
	/* HEAP approach
	double* lailow = new double[MAXGRID];
	double* laihigh = new double[MAXGRID];
	*/

	int lpjg_typeh[MAXGRID];
	int lpjg_typel[MAXGRID];


	/*
	// LPJG vegetation mapped onto H-TESSEL vegetation types.
	// 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
	*/

	// Fraction of these high and low vegetation types that occupy the gridcell
	double lpjg_frach[MAXGRID];
	double lpjg_fracl[MAXGRID];

	// if printOutputToFile is true (see config.h)
	FILE *ofp;
	float* mlailow = NULL;
	float* mlaihigh = NULL;
	float* mlpjg_frach = NULL;
	float* mlpjg_fracl = NULL;
	unsigned int* lpjg_typeh_yesterday = NULL;
	unsigned int* lpjg_typel_yesterday = NULL;

	if (printOutputToFile) {

		mlailow = new float[12*MAXGRID];
		mlaihigh = new float[12*MAXGRID];
		mlpjg_frach = new float[12*MAXGRID];
		mlpjg_fracl = new float[12*MAXGRID];
		lpjg_typeh_yesterday = new unsigned int[MAXGRID];
		lpjg_typel_yesterday = new unsigned int[MAXGRID];

		ofp = fopen("LPJ-GUESS_monthlyoutput.txt", "a");
	}


	// For TM5 coupling
	double co2_tm5[MAXGRID];
	double co2curr;
	double dcflux[MAXGRID];
	double dnpp[MAXGRID]; // dcflux does not include dnpp, Total, daily flux is dcflux + dnpp

	// Local storage for the date
	struct {
		int year;
		int sim_year;
		int month;
		int day;
		int time; // time in seconds per 24 hour (0-based from 0h00) 
	} eceDate;

	int globalprocs = 1;
	int local_procs = 1;
	int myrank = 0;

	#ifdef HAVE_MPI
	if (islpjgspinup) {
		globalprocs = get_num_processes();
		local_procs = globalprocs; // and myrank is still = 0 from above
		myrank = get_rank();
	}
	#endif

	dprintf("LPJ-GUESS rank %i: %i local process(es), and %i global process(es) = %i \n",myrank, local_procs, globalprocs);

	// Initialise and define OASIS coupling. 
	// Only called when this is a parallel run AND when islpjgspinup is false 
	// (and on non-Windows platforms) 
	if(!islpjgspinup && isparallel){ // All processes must call oasis_init_comp - see OASIS-MCT documentation 

		// *** OASIS-MCT - initialisation *** 

		dprintf("LPJ-GUESS: OASIS initialising... \n");
		localcomm = -99;
		int ierror=OasisCoupler::init(localcomm);
		if (ierror == 0)
			dprintf("LPJ-GUESS: OASIS initialised, returned localcomm = %i \n",localcomm);
		else 
			dprintf("LPJ-GUESS: OASIS NOT initialised, returned ierror = %i \n",ierror);

		// Get rank for this process
		myrank = get_rank_specific(localcomm);

		if ( error_flag ) {
			dprintf("LPJ-GUESS aborts due to an error during setup. Please see error messages above.\n");
			char routine_name[] = "runlpjguess.cpp";
			char abort_message[] ="Error during init";
			int fin_OK = OasisCoupler::abort(localcomm, routine_name,abort_message,666);
			return;
		}
		// *** OASIS-MCT - exchange fields with root only *** 

		// For coupling to root only
		dprintf("LPJ-GUESS: process %i calling OasisCoupler::create_couplcomm \n", myrank);
		ierror=OasisCoupler::create_couplcomm(myrank,localcomm,couplcomm);
		if (ierror == 0) {
			dprintf("LPJ-GUESS: OASIS create_couplcomm OK for rank %i\n",myrank);
			dprintf("LPJ-GUESS: couplcomm = %i for rank %i\n",couplcomm,myrank);
		} else {
			dprintf("LPJ-GUESS: OASIS create_couplcomm FAILED for rank %i, so terminating OASIS (returned ierror = %i) \n",myrank, ierror);
			int fin_OK = OasisCoupler::abort(localcomm, "runlpjguess.cpp","Error during OASIS-init",-1);
			return;
		}

		// *** OASIS-MCT - partition and variable definitions *** 

		// Partition and variables (dummy partition array {0,0,0} if not root) 
		ierror=OasisCoupler::init_part_defvar(NX_ATMO,NY_ATMO,activateTM5coupling,myrank);
		if (ierror == 0)
			dprintf("LPJ-GUESS: OASIS init_part_defvar OK for rank %i\n",myrank);
		else {
			dprintf("LPJ-GUESS: OASIS init_part_defvar FAILED for rank %i , so terminating OASIS (returned ierror = %i) \n",myrank, ierror);
			int fin_OK = OasisCoupler::abort(localcomm, "runlpjguess.cpp","Error during OASIS-init-part-def-var",-1);
			return;
		}

        // Wait until all processes, including 0, reach this point, so OASIS data will be available
		//dprintf("Before MPI_Barrier 0 on node %i \n",myrank);
		#ifdef HAVE_MPI
		if (local_procs>1)
			MPI_Barrier(localcomm);
		//dprintf("After MPI_Barrier 0 on node %i \n",myrank);
        #endif

	} else {

		dprintf("LPJ-GUESS: No OASIS initialisation. \n");
		localcomm = 0;
		couplcomm = 0;
 		#ifdef HAVE_MPI
			localcomm = MPI_COMM_WORLD; 
		#endif

	}

	//dprintf("LPJ-GUESS: localcomm %i \n",localcomm);

	globalprocs = get_num_global_processes();
	if (!islpjgspinup)
		local_procs = get_num_local_processes(localcomm);
	else
		local_procs = get_num_processes();

	// dprintf("LPJ-GUESS: %i local processes, and %i global processes = %i, myrank %i \n",local_procs, globalprocs, myrank);

    if (isparallel && !islpjgspinup) {

		xtring path;
        path.printf("./run%d", myrank+1); // ecev3 - was rank+1
		dprintf("\nMoving directory on node %d\n",myrank);

		if (change_directory(path) != 0) {
			fprintf(stderr, "Failed to change to run directory\n");
			return;
		}
	}

	// Wait until all processes, including 0, reach this point, so OASIS data will be available
	//dprintf("Before MPI_Barrier 1 (after OASIS-MCT initialisation) on node %i \n",myrank); 
	#ifdef HAVE_MPI
	if (isparallel && !islpjgspinup && local_procs>1)
		MPI_Barrier(localcomm);
	#endif
	//dprintf("After MPI_Barrier 1 on node %i \n",myrank);


	dprintf("EC-Earth - LPJ-GUESS coupling ");
	dprintf("\nfor %d timesteps of %d seconds ",NTIMESTEP,TIMESTEP);
	dprintf("starting from %04d-%02d-%02d 0h00\n\n",STARTDATE.year,STARTDATE.month,STARTDATE.day);

	// Set date and time for first timestep
	eceDate.year=STARTDATE.year;
	eceDate.month=STARTDATE.month;
	eceDate.day=STARTDATE.day;
	eceDate.time=0;
	eceDate.sim_year=0;



	// MAIN LOOP THROUGH TIMESTEPS

	int isfinal = 0; // set flag for the final step (used for saving the LPJG state)

	timestep=0;

	// Ensure we only run for one timestep when spinning up
	int timestepstorun = 0;

	if (islpjgspinup) {
		timestepstorun = 1;
		isfinal = 1;
	} else
		timestepstorun = NTIMESTEP;


	// Main time loop
	while (timestep<timestepstorun) {

		// Final step?
		if (timestep == NTIMESTEP-1) {
			// Usually on the last day of the year
			isfinal = 1;
			dprintf("runlpjguess: entered final timestep\n");
		}

		if (timestep < 2)
			dprintf("runlpjguess: timestep %i of %i (%i s)\n",timestep,NTIMESTEP,timestep*TIMESTEP);


		// Initialise CO2 concentration from file
		// Dataset starts in 1 BC, i.e. year 0, so index is year

		if (ifs_FIXEDYEARCO2 <= 0 && !ifs_A4xCO2 && !ifs_1PCTCO2) {
			// Read CO2 values from the array iff this isn't a fixed-year DECK experiment,
			// a 4*CO2 experiment, or a 1%/year CO2 experiment

			if (eceDate.year<FIRSTYEAR_CO2) {
				co2curr=co2[0]; // before 1 BC
			} 
			else if ( eceDate.year > NYEAR_CO2 ) {
				fail("No CO2 data available beyond %d \n",NYEAR_CO2);
			}
			else {
				co2curr=co2[eceDate.year]; // from year 0+
			}
		}

		// reset CO2 for CMIP6 experiments
		// Fixed CO2 of a certain year
		if ( ifs_FIXEDYEARCO2 > 0 ) {
			co2curr=co2[ifs_FIXEDYEARCO2];
		} 
		// abrupt 4xCo2 in 1850
		else if ( ifs_A4xCO2 ) {
			if ( eceDate.year >= CMIP6STARTYEAR ) {
				co2curr = 4. * co2[CMIP6STARTYEAR];
			}
		}
		// 1% / year (exponential) increase from 1850 on
		else if ( ifs_1PCTCO2 ) {
			if ( eceDate.year >= CMIP6STARTYEAR ) {
				co2curr = pow(1.01,eceDate.year-CMIP6STARTYEAR + 1 ) * co2[CMIP6STARTYEAR];  
			}
		} 
		
		// mid-Holocene value - ecev3 - T159
		// co2curr = 264.4; 

		// Reset on Jan 1?
		if (printOutputToFile && eceDate.day-1 == 0 && eceDate.month - 1 == 0) {
			for (int ii = 0; ii < MAXGRID; ii++) {
				// reset mLAI
				for (int mm = 0; mm < 12; mm++) {
					mlailow[ii * 12 + mm] = 0.0;
					mlaihigh[ii * 12 + mm] = 0.0;
					mlpjg_fracl[ii * 12 + mm] = 0.0;
					mlpjg_frach[ii * 12 + mm] = 0.0;
				}
			}
		}


		// dprintf("Before OASIS GET on node %i \n",myrank);
		// Call OASIS GET (if not spinning up)
		if (!islpjgspinup && myrank==0) {			
			call_coupler_get(timestep,temp,prec,vegl,vegl_type,vegh,snowc,snowd,stl,sml,swrad,lwrad,co2_tm5);
			if (timestep<20) dprintf("runlpjguess: node %i of %i returned from call_coupler_get \n",myrank,local_procs);
		}
		
		//dprintf("Before MPI_Barrier 2 (after call_coupler_get) on node %i \n",myrank);
        #ifdef HAVE_MPI
		if (!islpjgspinup && local_procs>1)
			MPI_Barrier(localcomm);
        #endif
		//dprintf("After MPI_Barrier 2 on node %i \n",myrank);


		// ecev3 - populate the co2_tm5 file with CO2 values from file if we are not coupled to TM5
		if (!activateTM5coupling)
			for (int jj=0; jj<MAXGRID; jj++) co2_tm5[jj]=co2curr;

		// RUN LPJ-GUESS today, for ALL cells!
		if (timestep < 2)
			dprintf("before runlpjguess_today: exchanging fields: timestep %i of %i (%i s)\n", timestep,NTIMESTEP,timestep*TIMESTEP);


		// ecev3 - could also wrap this function with the loop through cells and call GUESS directly. At the moment it's done in call_guess
		// NB! We send in month-1 and day-1, as the LPJG Date class has a base 0.
		runlpjguess_today(islpjgspinup,timestep,isfinal,eceDate.year,eceDate.sim_year,
				  eceDate.month-1,eceDate.day-1,eceDate.time,temp,prec,swrad,lwrad,snowc,
				  snowd,stl,sml,lailow,laihigh,lpjg_typeh,lpjg_frach,lpjg_typel,
				  lpjg_fracl,co2_tm5,dcflux,dnpp);
		
		//dprintf("Before MPI_Barrier 3 (after runlpjguess_today) on node %i \n",myrank);
		#ifdef HAVE_MPI
		if (!islpjgspinup && local_procs>1)
			MPI_Barrier(localcomm);
		#endif
		//dprintf("After MPI_Barrier 3 on node %i \n",myrank);

		// Call OASIS PUT (if not spinning up)
		if (!islpjgspinup && myrank==0) {			
			call_coupler_put(timestep,lailow,laihigh,lpjg_typeh,lpjg_frach,lpjg_typel,lpjg_fracl,dcflux,dnpp);
			if (timestep<20) dprintf("runlpjguess: node %i of %i returned from call_coupler_put \n",myrank,local_procs);
		}


		///////////////////////////////////////////////////////////////////////////////

		// Advance clock for next timestep
		// ecev3 - simplify perhaps??? Or move to a new subroutine

		if (timestep < 2)
			dprintf("runlpjguess: before clock advance on timestep %i at time %i on day %i of month % i of year %i\n",
				timestep,eceDate.time,eceDate.day,eceDate.month,eceDate.year); 


		// Need this in case the FIRST year is a leap year
		if (IFLEAPYEARS) {
        
			int yy = eceDate.year;

			if (!(yy%400)) 
				NDAYMONTH[1]=29; // e.g. year 2000
			else if (!(yy%100)) 
				NDAYMONTH[1]=28; // e.g. year 1900 is not a leap year
			else if (!(yy%4)) 
				NDAYMONTH[1]=29; // 
			else 
				NDAYMONTH[1]=28;
		}

		// Print? Only by rank = 0 on Dec 31 in a parallel run
		if (printOutputToFile && isparallel && !islpjgspinup && myrank==0) {
			
			// Run through this loop each day
			for (int ii = 0; ii < MAXGRID; ii++) {
				
				// Save mLAI
				mlailow[ii * 12 + eceDate.month - 1] += lailow[ii] / NDAYMONTH[eceDate.month - 1];
				mlaihigh[ii * 12 + eceDate.month - 1] += laihigh[ii] / NDAYMONTH[eceDate.month - 1];

				mlpjg_frach[ii * 12 + eceDate.month - 1] += lpjg_frach[ii] / NDAYMONTH[eceDate.month - 1];
				mlpjg_fracl[ii * 12 + eceDate.month - 1] += lpjg_fracl[ii] / NDAYMONTH[eceDate.month - 1];
				
				// Save veg. type and fraction 
				if (eceDate.day == 1) { 
					// Avoid doing this on Dec 31, as the type and fractions are updated then
					lpjg_typeh_yesterday[ii] = lpjg_typeh[ii];
					lpjg_typel_yesterday[ii] = lpjg_typel[ii];
				}

				// Dec 31st? Print!
				if (eceDate.month == 12 && eceDate.day==31) {

					EceCoord& c = gridlist[ii];
					double alat = (float)c.lat;
					double alon = (float)c.lon;

					fprintf(ofp, "%8.6f\t %8.6f\t %d\t", alon, alat, eceDate.year);
					int mm;

					for (mm = 0; mm < 12; mm++) {
						fprintf(ofp, "%8.5f\t", mlailow[ii * 12 + mm]);
					}
					for (mm = 0; mm < 12; mm++) {
						fprintf(ofp, "%8.5f\t", mlaihigh[ii * 12 + mm]);
					}

					for (mm = 0; mm < 12; mm++) {
						fprintf(ofp, "%8.5f\t", mlpjg_fracl[ii * 12 + mm]);
					}
					for (mm = 0; mm < 12; mm++) {
						fprintf(ofp, "%8.5f\t", mlpjg_frach[ii * 12 + mm]);
					}

					fprintf(ofp, "%d\t %d\t\n", lpjg_typel_yesterday[ii], lpjg_typeh_yesterday[ii]);
				}

			} // for
			fflush(ofp);
		}

		eceDate.time+=TIMESTEP;
		if (eceDate.time>=24*3600) { // Next day
			eceDate.time-=24*3600;
			eceDate.day+=1;
			if (eceDate.day>NDAYMONTH[eceDate.month-1]) {
				eceDate.day=1;
				eceDate.month+=1;
				if (eceDate.month>12) {
					eceDate.month=1;
					eceDate.year++;
					eceDate.sim_year++; 
					if (IFLEAPYEARS) {

						int yy = eceDate.year;

						if (!(yy%400)) 
							NDAYMONTH[1]=29; // e.g. year 2000
						else if (!(yy%100)) 
							NDAYMONTH[1]=28;
						else if (!(yy%4)) 
							NDAYMONTH[1]=29;
						else 
							NDAYMONTH[1]=28;		
					
					}
				}	
			}
		}

		timestep++;	
	}

	/* HEAP approach
	delete[] temp;
	delete[] prec;
	delete[] swrad;
	delete[] lwrad;
	delete[] vegl;
	delete[] vegh;
	delete[] snowc;
	delete[] snowd;
	delete[] lailow;
	delete[] laihigh;
	*/

	if (printOutputToFile) {

		delete[] mlailow;
		delete[] mlaihigh;
		delete[] mlpjg_frach;
		delete[] mlpjg_fracl;
		delete[] lpjg_typeh_yesterday;
		delete[] lpjg_typel_yesterday;

		fclose(ofp);
	}

	if(!islpjgspinup){ // All processes must call OASIS-MCT terminate
		// termination OASIS
		dprintf("LPJ-GUESS: OASIS terminating...\n");
		OasisCoupler::finalize();
		dprintf("LPJ-GUESS: OASIS terminated!\n");
	}

	dprintf ("Terminating ...\n");
	
}
 
// ecev3 - called from main in trunk
// int ecemain(int argc,char* argv[]) {
int ecemain(const CommandLineArguments& args) {

	dprintf("Running LPJ-GUESS - in ecemain() in eceframework.cpp ...\n");
	bool error_flag = false;
	
	// Read grids.nc, masks.nc, soilcd, timesteps and CO2
	read_input_data(error_flag);    

	if (!error_flag) {
		if (args.get_islpjgspinup())
			dprintf("ecemain(): read_input_data OK. Now calling runlpjguess for a spinup...\n");
		else
			dprintf("ecemain(): read_input_data OK. Now calling runlpjguess - NO spinup\n");
	}
	else {
		dprintf("ecemain(): error in read_input_data. Now calling runlpjguess for a spinup...Aborting after OASIS start\n");
	}
	// Run LPJ-GUESS
	runlpjguess(args.get_islpjgspinup(), args.get_parallel(),error_flag);
        
	dprintf("LPJ-GUESS: exiting - ecemain() in eceframework.cpp...\n");

	return 0;
}