eceinput.cpp 8.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320
  1. ///////////////////////////////////////////////////////////////////////////////////////
  2. /// \file eceinput.h
  3. /// \brief Input module copied from demoinput.h for EC-Earth data (ecev3). Paul Miller
  4. ///
  5. ///
  6. /// \author Ben Smith
  7. /// $Date: 2015-11-20 13:48:07 +0100 (fr, 20 nov 2015) $
  8. ///
  9. ///////////////////////////////////////////////////////////////////////////////////////
  10. #include "config.h"
  11. #include "eceinput.h"
  12. #include "driver.h"
  13. #include "outputchannel.h"
  14. #include <stdio.h>
  15. #include "eceframework.h" // ecev3 - for EceCoord & gridlist
  16. REGISTER_INPUT_MODULE("ece", EceInput)
  17. /*
  18. // Anonymous namespace for variables and functions with file scope
  19. namespace {
  20. // File names for temperature, precipitation, sunshine and soil code driver files
  21. xtring file_temp,file_prec,file_sun,file_soil;
  22. // LPJ soil code
  23. int soilcode;
  24. /// Interpolates monthly data to quasi-daily values.
  25. void interp_climate(double* mtemp, double* mprec, double* msun, double* mdtr,
  26. double* dtemp, double* dprec, double* dsun, double* ddtr) {
  27. interp_monthly_means_conserve(mtemp, dtemp);
  28. interp_monthly_totals_conserve(mprec, dprec, 0);
  29. interp_monthly_means_conserve(msun, dsun, 0, 100);
  30. interp_monthly_means_conserve(mdtr, ddtr, 0);
  31. }
  32. } // namespace
  33. */
  34. EceInput::EceInput()
  35. : nyear(1) {
  36. // Declare instruction file parameters
  37. // declare_parameter("nyear", &nyear, 1, 10000, "Number of simulation years to run after spinup");
  38. }
  39. bool EceInput::read_from_file(Coord coord, xtring fname, const char* format,
  40. double monthly[12], bool soil /* = false */) {
  41. /*
  42. ecev3 - unreferenced
  43. double dlon, dlat;
  44. int elev;
  45. */
  46. FILE* in = fopen(fname, "r");
  47. if (!in) {
  48. fail("readenv: could not open %s for input", (char*)fname);
  49. }
  50. bool foundgrid = false;
  51. /*
  52. while (!feof(in) && !foundgrid) {
  53. if (!soil) {
  54. readfor(in, format, &dlon, &dlat, &elev, monthly);
  55. } else {
  56. readfor(in, format, &dlon, &dlat, &soilcode);
  57. }
  58. foundgrid = equal(coord.lon, dlon) && equal(coord.lat, dlat);
  59. }
  60. fclose(in);
  61. if (!foundgrid) {
  62. dprintf("readenv: could not find record for (%g,%g) in %s",
  63. coord.lon, coord.lat, (char*)fname);
  64. }
  65. */
  66. return foundgrid;
  67. }
  68. bool EceInput::readenv(Coord coord, long& seed) {
  69. // Searches for environmental data in driver temperature, precipitation,
  70. // sunshine and soil code files for the grid cell whose coordinates are given by
  71. // 'coord'. Data are written to arrays mtemp, mprec, msun and the variable
  72. // soilcode, which are defined as global variables in this file
  73. // The temperature, precipitation and sunshine files (Cramer & Leemans,
  74. // unpublished) should be in ASCII text format and contain one-line records for the
  75. // climate (mean monthly temperature, mean monthly percentage sunshine or total
  76. // monthly precipitation) of a particular 0.5 x 0.5 degree grid cell. Elevation is
  77. // also included in each grid cell record.
  78. // The following sample record from the temperature file:
  79. // " -4400 8300 293-298-311-316-239-105 -7 26 -3 -91-184-239-277"
  80. // corresponds to the following data:
  81. // longitude 44 deg (-=W)
  82. // latitude 83 deg (+=N)
  83. // elevation 293 m
  84. // mean monthly temperatures (deg C) -29.8 (Jan), -31.1 (Feb), ..., -27.7 (Dec)
  85. // The following sample record from the precipitation file:
  86. // " 12750 -200 223 190 165 168 239 415 465 486 339 218 162 149 180"
  87. // corresponds to the following data:
  88. // longitude 127.5 deg (+=E)
  89. // latitude 20 deg (-=S)
  90. // elevation 223 m
  91. // monthly precipitation sum (mm) 190 (Jan), 165 (Feb), ..., 180 (Dec)
  92. // The following sample record from the sunshine file:
  93. // " 2600 7000 293 0 20 38 37 31 28 28 25 21 17 9 7"
  94. // corresponds to the following data:
  95. // longitude 26 deg (+=E)
  96. // latitude 70 deg (+=N)
  97. // elevation 293 m
  98. // monthly mean %age of full sunshine 0 (Jan), 20 (Feb), ..., 7 (Dec)
  99. // The LPJ soil code file is in ASCII text format and contains one-line records for
  100. // each grid cell in the form:
  101. // <lon> <lat> <soilcode>
  102. // where <lon> = longitude as a floating point number (-=W, +=E)
  103. // <lat> = latitude as a floating point number (-=S, +=N)
  104. // <soilcode> = integer in the range 0 (no soil) to 9 (see function
  105. // soilparameters in driver module)
  106. // The fields in each record are separated by spaces
  107. /*
  108. ecev3 - unreferenced
  109. double mtemp[12]; // monthly mean temperature (deg C)
  110. double mprec[12]; // monthly precipitation sum (mm)
  111. double msun[12]; // monthly mean percentage sunshine values
  112. */
  113. double mwet[12]={31,28,31,30,31,30,31,31,30,31,30,31}; // number of rain days per month
  114. /*
  115. double mdtr[12]; // monthly mean diurnal temperature range (oC)
  116. for(int m=0; m<12; m++) {
  117. mdtr[m] = 0.;
  118. if (ifbvoc) {
  119. dprintf("WARNING: No data available for dtr in sample data set!\nNo daytime temperature correction for BVOC calculations applied.");
  120. }
  121. }
  122. bool gridfound = read_from_file(coord, file_temp, "f6.2,f5.2,i4,12f4.1", mtemp);
  123. if(gridfound)
  124. gridfound = read_from_file(coord, file_prec, "f6.2,f5.2,i4,12f4", mprec);
  125. if(gridfound)
  126. gridfound = read_from_file(coord, file_sun, "f6.2,f5.2,i4,12f3", msun);
  127. if(gridfound)
  128. gridfound = read_from_file(coord, file_soil, "f,f,i", msun, true); // msun is not used here: just dummy
  129. if(gridfound) {
  130. // Interpolate monthly values for environmental drivers to daily values
  131. // (relevant daily values will be sent to the framework each simulation
  132. // day in function getclimate, below)
  133. interp_climate(mtemp, mprec, msun, mdtr, dtemp, dprec, dsun, ddtr);
  134. // Recalculate precipitation values using weather generator
  135. // (from Dieter Gerten 021121)
  136. prdaily(mprec, dprec, mwet, seed);
  137. }
  138. return gridfound;
  139. */
  140. return true; // ecev3
  141. }
  142. void EceInput::init() {
  143. // DESCRIPTION
  144. // Initialises input (e.g. opening files), and reads in the gridlist
  145. //
  146. // Reads list of grid cells and (optional) description text from grid list file
  147. // This file should consist of any number of one-line records in the format:
  148. // <longitude> <latitude> [<description>]
  149. double dlon,dlat;
  150. bool eof=false;
  151. xtring descrip;
  152. // Read list of grid coordinates and store in global Coord object 'gridlist'
  153. // Retrieve name of grid list file as read from ins file
  154. // This will be *created* in eceframework.cpp
  155. xtring file_ecegridlist=param["file_gridlist"].str;
  156. FILE* in_grid=fopen(file_ecegridlist,"r");
  157. if (!in_grid) fail("EceInput::init: could not open %s for input",(char*)file_ecegridlist);
  158. gridlist.killall();
  159. first_call = true;
  160. while (!eof) {
  161. // Read next record in file
  162. eof=!readfor(in_grid,"f,f,a#",&dlon,&dlat,&descrip);
  163. if (!eof && !(dlon==0.0 && dlat==0.0)) { // ignore blank lines at end (if any)
  164. Coord& c=gridlist.createobj(); // add new coordinate to grid list
  165. c.lon=(float)dlon;
  166. c.lat=(float)dlat;
  167. // c.descrip=descrip;
  168. }
  169. }
  170. fclose(in_grid);
  171. // Open landcover files
  172. landcover_input.init();
  173. // Open management files
  174. management_input.init();
  175. // ecev3 - set the first calendar year to 0
  176. date.set_first_calendar_year(0);
  177. // Set timers
  178. tprogress.init();
  179. tmute.init();
  180. tprogress.settimer();
  181. tmute.settimer(MUTESEC);
  182. }
  183. bool EceInput::getgridcell(Gridcell& gridcell) {
  184. // See base class for documentation about this function's responsibilities
  185. // Select coordinates for next grid cell in linked list
  186. bool gridfound = false;
  187. bool LUerror = false;
  188. if(run_landcover) {
  189. LUerror = landcover_input.loadlandcover(gridcell.get_lon(), gridcell.get_lat());
  190. if(!LUerror)
  191. LUerror = management_input.loadmanagement(gridcell.get_lon(), gridcell.get_lat());
  192. }
  193. if(!LUerror)
  194. return true;
  195. else
  196. return false;
  197. }
  198. void EceInput::getlandcover(Gridcell& gridcell) {
  199. landcover_input.getlandcover(gridcell);
  200. if (! ( gridcell.fixedLUafter >= 0 && date.year >= gridcell.fixedLUafter )){
  201. landcover_input.get_land_transitions(gridcell);
  202. }
  203. }
  204. bool EceInput::getclimate(Gridcell& gridcell) {
  205. // See base class for documentation about this function's responsibilities
  206. /*
  207. Climate& climate = gridcell.climate;
  208. // Send environmental values for today to framework
  209. climate.dndep = ndep / (365.0 * 10000.0);
  210. climate.co2 = co2;
  211. climate.temp = dtemp[date.day];
  212. climate.prec = dprec[date.day];
  213. climate.insol = dsun[date.day];
  214. // bvoc
  215. climate.dtr=ddtr[date.day];
  216. // First day of year only ...
  217. if (date.day == 0) {
  218. // Return false if last year was the last for the simulation
  219. if (date.year==nyear_spinup+nyear) return false;
  220. // Progress report to user and update timer
  221. if (tmute.getprogress()>=1.0) {
  222. progress=(double)(gridlist.getobj().id*(nyear_spinup+nyear)
  223. +date.year)/(double)(gridlist.nobj*(nyear_spinup+nyear));
  224. tprogress.setprogress(progress);
  225. dprintf("%3d%% complete, %s elapsed, %s remaining\n",(int)(progress*100.0),
  226. tprogress.elapsed.str,tprogress.remaining.str);
  227. tmute.settimer(MUTESEC);
  228. }
  229. }
  230. */
  231. return true;
  232. }
  233. EceInput::~EceInput() {
  234. // Performs memory deallocation, closing of files or other "cleanup" functions.
  235. }