ecearth3/sources/lpjg/cru/guessio/cruinput.cpp

///////////////////////////////////////////////////////////////////////////////////////
/// \file cruinput.cpp
/// \brief LPJ-GUESS input module for CRU-NCEP data set
///
/// This input module reads in CRU-NCEP climate data in a customised binary format.
/// The binary files contain CRU-NCEP half-degree global historical climate data
/// for 1901-2015.
///
/// \author Ben Smith
/// $Date: 2013-11-22 10:56:18 +0100 (Fri, 22 Nov 2013) $
///
///////////////////////////////////////////////////////////////////////////////////////

#include "config.h"
#include "cruinput.h"

#include "driver.h"
#include "parameters.h"
#include <stdio.h>
#include <utility>
#include <vector>
#include <algorithm>
#include <netcdf.h>

REGISTER_INPUT_MODULE("cru_ncep", CRUInput)

// Anonymous namespace for variables and functions with file scope
namespace {

xtring file_cru;
xtring file_cru_misc;

/// Interpolates monthly data to quasi-daily values.
void interp_climate(double* mtemp, double* mprec, double* msun, double* mdtr,
					double* dtemp, double* dprec, double* dsun, double* ddtr) {
	interp_monthly_means_conserve(mtemp, dtemp, date.is_leap(date.year));
	interp_monthly_totals_conserve(mprec, dprec, date.is_leap(date.year), 0);
	interp_monthly_means_conserve(msun, dsun, date.is_leap(date.year), 0);
	interp_monthly_means_conserve(mdtr, ddtr, date.is_leap(date.year), 0);
}

} // namespace

CRUInput::CRUInput()
	: searchradius(0),
	  spinup_mtemp(NYEAR_SPINUP_DATA),
	  spinup_mprec(NYEAR_SPINUP_DATA),
	  spinup_msun(NYEAR_SPINUP_DATA),
	  spinup_mfrs(NYEAR_SPINUP_DATA),
	  spinup_mwet(NYEAR_SPINUP_DATA),
	  spinup_mdtr(NYEAR_SPINUP_DATA) {

	// Declare instruction file parameters

	declare_parameter("searchradius", &searchradius, 0, 100,
		"If specified, CRU data will be searched for in a circle");
}


void CRUInput::init() {

	// DESCRIPTION
	// Initialises input (e.g. opening files), and reads in the gridlist

	//
	// Reads list of grid cells and (optional) description text from grid list file
	// This file should consist of any number of one-line records in the format:
	//   <longitude> <latitude> [<description>]

	double dlon,dlat;
	bool eof=false;
	xtring descrip;

	// Read list of grid coordinates and store in global Coord object 'gridlist'

	// Retrieve name of grid list file as read from ins file
	xtring file_gridlist=param["file_gridlist"].str;

	FILE* in_grid=fopen(file_gridlist,"r");
	if (!in_grid) fail("initio: could not open %s for input",(char*)file_gridlist);

	file_cru=param["file_cru"].str;
	file_cru_misc=param["file_cru_misc"].str;

	gridlist.killall();
	first_call = true;	

	while (!eof) {

		// Read next record in file
		eof=!readfor(in_grid,"f,f,a#",&dlon,&dlat,&descrip);

		if (!eof && !(dlon==0.0 && dlat==0.0)) { // ignore blank lines at end (if any)
			Coord& c=gridlist.createobj(); // add new coordinate to grid list

			c.lon=dlon;
			c.lat=dlat;
			c.descrip=descrip;
		}
	}


	fclose(in_grid);

	// Read CO2 data from file
	co2.load_file(param["file_co2"].str);

	// Open landcover files
	landcover_input.init();
	// Open management files
	management_input.init();

	date.set_first_calendar_year(FIRSTHISTYEAR - nyear_spinup);
	// Set timers
	tprogress.init();
	tmute.init();

	tprogress.settimer();
	tmute.settimer(MUTESEC);
}


void CRUInput::get_monthly_ndep(int calendar_year,
                                double* mndrydep,
                                double* mnwetdep) {

	ndep.get_one_calendar_year(calendar_year,
	                           mndrydep, mnwetdep);
}


void CRUInput::adjust_raw_forcing_data(double lon,
                                       double lat,
                                       double hist_mtemp[NYEAR_HIST][12],
                                       double hist_mprec[NYEAR_HIST][12],
                                       double hist_msun[NYEAR_HIST][12]) {

	// The default (base class) implementation does nothing here.
}


bool CRUInput::getgridcell(Gridcell& gridcell) {

	// See base class for documentation about this function's responsibilities

	int soilcode;
	int elevation;

	// Make sure we use the first gridcell in the first call to this function,
	// and then step through the gridlist in subsequent calls.
	if (first_call) {
		gridlist.firstobj();

		// Note that first_call is static, so this assignment is remembered
		// across function calls.
		first_call = false;
	}
	else gridlist.nextobj();

	if (gridlist.isobj) {

		bool gridfound = false;
		bool LUerror = false;
		double lon;
		double lat;

		while(!gridfound) {

			if(gridlist.isobj) {

				lon = gridlist.getobj().lon;
				lat = gridlist.getobj().lat;
				gridfound = CRU_TS30::findnearestCRUdata(searchradius, file_cru, lon, lat, soilcode,
				                                         hist_mtemp, hist_mprec, hist_msun);
			  
				if (gridfound) // Get more historical CRU data for this grid cell
					gridfound = CRU_TS30::searchcru_misc(file_cru_misc, lon, lat, elevation,
					                                     hist_mfrs, hist_mwet, hist_mdtr);

        // back to gridlist lon lat, since the generated LU input are for the ece lon lats
				double lon2 = gridlist.getobj().lon;
				double lat2 = gridlist.getobj().lat;
				if (run_landcover && gridfound) {
					LUerror = landcover_input.loadlandcover(lon2, lat2);
					if(!LUerror)
						LUerror = management_input.loadmanagement(lon2, lat2);
				}

				if(!gridfound || LUerror) {
					if(!gridfound)
						dprintf("\nError: could not find stand at (%g,%g) in climate data files\n\n", gridlist.getobj().lon,gridlist.getobj().lat);
					else if(LUerror)
						dprintf("\nError: could not find stand at (%g,%g) in landcover/management data file(s)\n\n", gridlist.getobj().lon,gridlist.getobj().lat);
					gridfound = false;
					gridlist.nextobj();
				}
			}
			else return false;
		}

		// Give sub-classes a chance to modify the data
		adjust_raw_forcing_data(gridlist.getobj().lon,
		                        gridlist.getobj().lat,
		                        hist_mtemp, hist_mprec, hist_msun);

		// Build spinup data sets
		spinup_mtemp.get_data_from(hist_mtemp);
		spinup_mprec.get_data_from(hist_mprec);
		spinup_msun.get_data_from(hist_msun);

		// Detrend spinup temperature data
		spinup_mtemp.detrend_data();

		// guess2008 - new spinup data sets
		spinup_mfrs.get_data_from(hist_mfrs);
		spinup_mwet.get_data_from(hist_mwet);
		spinup_mdtr.get_data_from(hist_mdtr);

		// We wont detrend dtr for now. Partly because dtr is at the moment only
		// used for BVOC, so what happens during the spinup is not affecting
		// results in the period thereafter, and partly because the detrending
		// can give negative dtr values.
		//spinup_mdtr.detrend_data();


		dprintf("\nCommencing simulation for stand at (%g,%g)",gridlist.getobj().lon,
			gridlist.getobj().lat);
		if (gridlist.getobj().descrip!="") dprintf(" (%s)\n\n",
			(char*)gridlist.getobj().descrip);
		else dprintf("\n\n");
		
		// Tell framework the coordinates of this grid cell
		gridcell.set_coordinates(gridlist.getobj().lon, gridlist.getobj().lat);
		
		// Get nitrogen deposition data. 
		/* Since the historic data set does not reach decade 2010-2019, 
		 * we need to use the RCP data for the last decade. */
		ndep.getndep(param["file_ndep"].str, lon, lat, Lamarque::RCP60);

		// The insolation data will be sent (in function getclimate, below)
		// as incoming shortwave radiation, averages are over 24 hours
		
		gridcell.climate.instype = SWRAD_TS;

		// Tell framework the soil type of this grid cell
		soilparameters(gridcell.soiltype,soilcode);

		// For Windows shell - clear graphical output
		// (ignored on other platforms)
		
		clear_all_graphs();

		return true; // simulate this stand
	}

	return false; // no more stands
}


void CRUInput::getlandcover(Gridcell& gridcell) {

	landcover_input.getlandcover(gridcell);
	landcover_input.get_land_transitions(gridcell);
}

int getndepindex(double lon, double lat) {

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

	const char* filevar = "drynhx2.nc";

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

	// AMIP
	full_path += filevar;

	// Open NC file
	int status, ncid;
	status = nc_open(full_path.c_str(), NC_NOWRITE, &ncid);
	if (status != NC_NOERR)
		fail("Error reading dry NHx deposition data before LPJ-GUESS spinup. Quitting... \n");

	// Get the ID for lon
	int lon_id;
	status = nc_inq_varid(ncid, "lon", &lon_id);
	if (status != NC_NOERR)
		fail("Error reading longitude in dry NHx deposition data before LPJ-GUESS spinup. Quitting... \n");

	// Get the ID for lat
	int lat_id;
	status = nc_inq_varid(ncid, "lat", &lat_id);
	if (status != NC_NOERR)
		fail("Error reading latitude in dry NHx deposition data before LPJ-GUESS spinup. Quitting... \n");

	/*size_t latlength;
	status = nc_inq_dimlen(ncid, lat_id, &latlength);
	if (status != NC_NOERR)
	return false;*/

	int ndepindex;

	// Read data
	for (ndepindex = 0; ndepindex < 59191; ndepindex++) {
		size_t index[] = { ndepindex };

		double latval, lonval;

		status = nc_get_var1_double(ncid, lat_id, index, &latval);
		status = nc_get_var1_double(ncid, lon_id, index, &lonval);

		if (abs(latval-lat) < 0.000001 && abs(lonval-lon) < 0.000001)
			break;
	}
	if (ndepindex == 59191)
		fail("Could not find Lon %g and Lat %g in dry NHx deposition data. Quitting... \n", lon, lat);

	return ndepindex;
}


bool getmonthlyNdepforcing(const char* varname, int year, double data[12], int ndepindex) {

	// 1st (".nc") or 2nd order remapping file?
	const char* 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;
	full_path += 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 mm = 0; mm < 12; mm++) {
		size_t index[] = { year * 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[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

  /// ecev3 - Get nitrogen deposition data for this gridcell
  /**
  *  Determine N deposition forcing for this cell and year.
  *
  *  \param ndepfile			netcdf file containing the N dep. data. Read from guess.ins
  */
bool getndepCMIP6(const char* ndepfile, int ndep_index, int calendar_year, double mndrydep[12], double mnwetdep[12]) {

	double mdrynhx[12], mdrynoy[12], mwetnhx[12], mwetnoy[12];
	int maxyears = CMIP6ENDYEAR - CMIP6STARTYEAR;

	int ndep_year = min(maxyears, max(0, calendar_year - CMIP6STARTYEAR));

	// Monthly dry NH4 deposition (kg m-2 s-1)
	const char* vartoread = "drynhx";
	bool dataOK = getmonthlyNdepforcing(vartoread, ndep_year, mdrynhx, 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";
	dataOK = getmonthlyNdepforcing(vartoread, ndep_year, mwetnhx, 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";
	dataOK = getmonthlyNdepforcing(vartoread, ndep_year, mdrynoy, 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";
	dataOK = getmonthlyNdepforcing(vartoread, ndep_year, mwetnoy, ndep_index);
	if (!dataOK) {
		dprintf("Error reading wet NOy deposition data before LPJ-GUESS spinup. Quitting... \n");
		return false;
	}

	// Add up dry and wet deposition
	for (int m = 0; m < 12; m++) {
		mndrydep[m] = mdrynhx[m] + mdrynoy[m];
		mnwetdep[m] = mwetnhx[m] + mwetnoy[m];
	}

	return dataOK;
}


bool CRUInput::getclimate(Gridcell& gridcell) {

	// See base class for documentation about this function's responsibilities

	double progress;

	Climate& climate = gridcell.climate;

	if (date.day == 0) {

		// First day of year ...

		// Extract N deposition to use for this year,
		// monthly means to be distributed into daily values further down
		double mndrydep[12], mnwetdep[12];

		if (!ECEARTHWITHCRUNCEP) {
			ndep.get_one_calendar_year(date.year - nyear_spinup + FIRSTHISTYEAR,
				mndrydep, mnwetdep);
		}
		else {
			if (date.year == 0)
				gridcell.ndep_index = getndepindex(gridcell.get_lon(), gridcell.get_lat());

			bool Ndepdatafound = getndepCMIP6(param["ndep_cmip6_dir"].str, gridcell.ndep_index, date.year - nyear_spinup + FIRSTHISTYEAR, mndrydep, mnwetdep);
			if (!Ndepdatafound) {
				dprintf("Could not read N deposition files for cell with lon=%f, lat=%f\n", gridcell.get_lon(), gridcell.get_lat());
				fail();
			}
		}

		if (date.year < nyear_spinup) {

			// During spinup period

			if(date.year == state_year && restart) {

				int year_offset = state_year % NYEAR_SPINUP_DATA;

				for (int y=0;y<year_offset;y++) {
					spinup_mtemp.nextyear();
					spinup_mprec.nextyear();
					spinup_msun.nextyear();
					spinup_mfrs.nextyear();
					spinup_mwet.nextyear();
					spinup_mdtr.nextyear();
				}
			}

			int m;
			double mtemp[12],mprec[12],msun[12];
			double mfrs[12],mwet[12],mdtr[12];

			for (m=0;m<12;m++) {
				mtemp[m] = spinup_mtemp[m];
				mprec[m] = spinup_mprec[m];
				msun[m]	 = spinup_msun[m];

				mfrs[m] = spinup_mfrs[m];
				mwet[m] = spinup_mwet[m];
				mdtr[m] = spinup_mdtr[m];
			}

			// Interpolate monthly spinup data to quasi-daily values
			interp_climate(mtemp,mprec,msun,mdtr,dtemp,dprec,dsun,ddtr);

			// Only recalculate precipitation values using weather generator
			// if rainonwetdaysonly is true. Otherwise we assume that it rains a little every day.
			if (ifrainonwetdaysonly) {
				// (from Dieter Gerten 021121)
				prdaily(mprec, dprec, mwet, gridcell.seed);
			}

			spinup_mtemp.nextyear();
			spinup_mprec.nextyear();
			spinup_msun.nextyear();

			spinup_mfrs.nextyear();
			spinup_mwet.nextyear();
			spinup_mdtr.nextyear();

		}
		else if (date.year < nyear_spinup + NYEAR_HIST) {

			// Historical period

			// Interpolate this year's monthly data to quasi-daily values
			interp_climate(hist_mtemp[date.year-nyear_spinup],
				hist_mprec[date.year-nyear_spinup],hist_msun[date.year-nyear_spinup],
					   hist_mdtr[date.year-nyear_spinup],
				       dtemp,dprec,dsun,ddtr);

			// Only recalculate precipitation values using weather generator
			// if ifrainonwetdaysonly is true. Otherwise we assume that it rains a little every day.
			if (ifrainonwetdaysonly) {
				// (from Dieter Gerten 021121)
				prdaily(hist_mprec[date.year-nyear_spinup], dprec, hist_mwet[date.year-nyear_spinup], gridcell.seed);
			}
		}
		else {
			// Return false if last year was the last for the simulation
			return false;
		}

		for (int dd = 0; dd < Date::MAX_YEAR_LENGTH; dd++) dndep[dd] = 0.0;

		// Distribute N deposition
		if (!ECEARTHWITHCRUNCEP) {
			distribute_ndep(mndrydep, mnwetdep, dprec, date.ndaymonth, dndep);
		}
		else {
			// assuming it rains a bit every day
			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_ndep(mndrydep, mnwetdep, drizzle, date.ndaymonth, dndep);
		}
	}

	// Send environmental values for today to framework

	climate.co2 = co2[FIRSTHISTYEAR + date.year - nyear_spinup];

	if (date.day == 0) {
		for (int mth = 0; mth < 12; mth++) {
			gridcell.climate.mco2[mth] = 0.0;
		}
	}
	gridcell.climate.mco2[date.month] += climate.co2 / date.ndaymonth[date.month];

	climate.temp  = dtemp[date.day];
	climate.prec  = dprec[date.day];
	climate.insol = dsun[date.day];

	// Nitrogen deposition
	climate.dndep = dndep[date.day];

	// bvoc
	if(ifbvoc){
	  climate.dtr = ddtr[date.day];
	}

	// First day of year only ...

	if (date.day == 0) {

		// Progress report to user and update timer

		if (tmute.getprogress()>=1.0) {
			progress=(double)(gridlist.getobj().id*(nyear_spinup+NYEAR_HIST)
				+date.year)/(double)(gridlist.nobj*(nyear_spinup+NYEAR_HIST));
			tprogress.setprogress(progress);
			dprintf("%3d%% complete, %s elapsed, %s remaining\n",(int)(progress*100.0),
				tprogress.elapsed.str,tprogress.remaining.str);
			tmute.settimer(MUTESEC);
		}
	}

	return true;
}


CRUInput::~CRUInput() {

	// Performs memory deallocation, closing of files or other "cleanup" functions.

}


///////////////////////////////////////////////////////////////////////////////////////
// REFERENCES
// Lamarque, J.-F., Kyle, G. P., Meinshausen, M., Riahi, K., Smith, S. J., Van Vuuren,
//   D. P., Conley, A. J. & Vitt, F. 2011. Global and regional evolution of short-lived
//   radiatively-active gases and aerosols in the Representative Concentration Pathways.
//   Climatic Change, 109, 191-212.
// Nakai, T., Sumida, A., Kodama, Y., Hara, T., Ohta, T. (2010). A comparison between
//   various definitions of forest stand height and aerodynamic canopy height.
//   Agricultural and Forest Meteorology, 150(9), 1225-1233