123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760 |
- #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"
- static char file_co2[] = "mole_fraction_of_carbon_dioxide_in_air_input4MIPs_lpjg.nc";
- #ifdef GRID_T159
- static const int MAXGRID = 10407;
-
- static const int NX_ATMO=35718;
-
- static char parn_lon[]="L080.lon";
-
- static char parn_lat[]="L080.lat";
-
- static char parn_mask[]="L080.msk";
-
- static char parn_areas[] = "L080.srf";
-
- static char filen_soilcd[]="slt_T159.nc";
-
- static char parn_lon_cor[] = "L080.clo";
-
- static char parn_lat_cor[] = "L080.cla";
-
- #endif
- #ifdef GRID_T255
- static const int MAXGRID=25799;
-
- static const int NX_ATMO=88838;
-
- static char parn_lon[]="L128.lon";
-
- static char parn_lat[]="L128.lat";
-
- static char parn_mask[]="L128.msk";
-
- static char parn_areas[]="L128.srf";
-
- static char filen_soilcd[]="slt_T255.nc";
-
- static char parn_lon_cor[] = "L128.clo";
-
- static char parn_lat_cor[] = "L128.cla";
-
- #endif
- static char filen_grid[] = "grids.nc";
-
- static char filen_areas[]="areas.nc";
-
- static char filen_mask[] = "masks.nc";
-
-
- static char filen_cvl[]="cvl.nc";
-
- static char parn_cvl[]="var27";
-
- static char filen_cvh[]="cvh.nc";
-
- static char parn_cvh[]="var28";
-
- static char filen_tvl[]="tvl.nc";
-
- static char parn_tvl[]="var29";
-
- static char filen_tvh[]="tvh.nc";
-
- static char parn_tvh[]="var30";
-
- static char parn_soilcd[] = "var43";
-
- static char file_timesteps[]=".//lpjg_steps.rc";
-
- int resolution;
-
- bool activateTM5coupling;
-
- const int NYEAR_ECEARTH=271;
-
-
- const int FIRSTHISTYEAR=1840;
-
- const int NYEAR_CO2=2501;
- const int FIRSTYEAR_CO2=0;
- bool ifs_A4xCO2 = false;
-
- bool bgc_1pctCO2 = false;
-
- int ifs_FIXEDYEARCO2 = -1;
-
- int fixedNDepafter = -1;
-
- int fixedLUafter = -1;
-
- int TIMESTEP = 86400;
- struct{
- int year;
- int month;
- int day;
- } STARTDATE={1850,1,1};
- int NTIMESTEP=1;
- const bool IFLEAPYEARS=true;
- static int NDAYMONTH[12]={31,28,31,30,31,30,31,31,30,31,30,31};
-
- static int ngridcell;
-
- static double co2[NYEAR_CO2];
-
- static int soilcode[MAXGRID];
-
-
- static const int NY_ATMO=1;
- static float coup_lsmask[NX_ATMO][NY_ATMO];
- 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;
- 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];
- 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];
- static double coup_co2_field[NX_ATMO*NY_ATMO];
- static double coup_dcflux_nat[NX_ATMO*NY_ATMO];
- static double coup_dcflux_ant[NX_ATMO*NY_ATMO];
- static double coup_dnpp[NX_ATMO*NY_ATMO];
- static double coup_co2;
- struct EceCoord {
-
-
-
-
- int id;
- float lon;
- float lat;
- int soilcode;
- int IFSindex;
- float area;
- float cvl, cvh;
- int tvl, tvh;
- int ndep_index;
-
-
- 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));
-
- }
- }
- 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";
-
- status = nc_inq_varid (ncid, lonstr, &lon_id);
- if (status != NC_NOERR) handle_error(status);
-
- status = nc_inq_varid (ncid, latstr, &lat_id);
- if (status != NC_NOERR) handle_error(status);
-
-
- status = nc_get_var_double(ncid, lon_id, lon);
- if (status != NC_NOERR) handle_error(status);
-
- 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;
-
- status = nc_inq_varid (ncid, varname, &var_id);
- if (status != NC_NOERR) handle_error(status);
-
-
-
-
- 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;
-
- status = nc_inq_varid (ncid, varname, &var_id);
- if (status != NC_NOERR) handle_error(status);
- static size_t start[] = {0, 0, 0};
- static size_t count[] = {corners, lons, lats};
-
- 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;
-
- status = nc_inq_varid (ncid, varname, &var_id);
- if (status != NC_NOERR) handle_error(status);
-
-
- 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;
- }
- }
- static ListArray_id<EceCoord> gridlist;
-
-
- 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;
-
- eof=!readfor(in_ts,"i",&TIMESTEP);
-
- if (TIMESTEP>24*3600) {
- dprintf("Maximum allowable timestep is %d seconds (24 hours)\n",24*3600);
- fclose(in_ts);
- return false;
- }
- if (!eof) {
-
- readfor(in_ts,"i",&STARTDATE.year);
- readfor(in_ts,"i",&STARTDATE.month);
- readfor(in_ts,"i",&STARTDATE.day);
-
- 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)) {
- if (STARTDATE.day>NDAYMONTH[1]+1)
- valid=false;
- }
- else if (STARTDATE.day>NDAYMONTH[STARTDATE.month-1])
- valid=false;
- if (!valid) {
- dprintf("ERROR: Invalid start date %d-%02d-%02d\n", STARTDATE.year,STARTDATE.month,STARTDATE.day);
- fclose(in_ts);
- is_err = true;
- return false;
- }
- if (STARTDATE.year<0 ) {
- dprintf("ERROR: Specified start date %d-%02d-%02d is outside range of available data\n",
- STARTDATE.year,STARTDATE.month,STARTDATE.day);
- fclose(in_ts);
- is_err = true;
- return false;
- }
-
- eof=!readfor(in_ts,"i",&NTIMESTEP);
-
- readfor(in_ts,"i",&resolution);
- if (resolution != 159 && resolution != 255) {
- dprintf("ERROR: Invalid resolution %d\n", resolution);
- fclose(in_ts);
- is_err = true;
- return false;
- }
-
- int TM5couple = -1;
- readfor(in_ts,"i",&TM5couple);
- if (TM5couple==1)
- activateTM5coupling = true;
- else if (TM5couple==0)
- activateTM5coupling = false;
- else {
- dprintf("ERROR: Invalid TM5 coupling option %d\n", TM5couple);
- fclose(in_ts);
- is_err = true;
- return false;
- }
-
-
- xtring sdum[1];
- string is_true_l = "t";
- string is_true_u = "T";
- eof = !readfor(in_ts,"a1",&sdum);
- if ( eof ) {
- dprintf("ERROR: 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);
- }
-
-
- eof = !readfor(in_ts,"a1",&sdum);
- if ( eof ) {
- dprintf("ERROR: 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]) ) {
- bgc_1pctCO2 = 1;
- }
- else {
- bgc_1pctCO2 = 0;
- }
- dprintf("1PCTCO2 set to %i \n",bgc_1pctCO2);
- }
-
- eof = !readfor(in_ts,"i",&ifs_FIXEDYEARCO2);
- if ( eof ) {
- dprintf("ERROR: Entry for FIXEDYEARCO2 not found in %s \n", file_timesteps);
- is_err = true;
- }
- else {
- dprintf("CO2 is set fix at %d \n",ifs_FIXEDYEARCO2);
- }
-
- eof = !readfor(in_ts,"i",&fixedNDepafter);
- if ( eof ) {
- dprintf("ERROR: Entry for fixedNDepafter year not found in %s \n", file_timesteps);
- is_err = true;
- }
-
- eof = !readfor(in_ts,"i",&fixedLUafter);
- if ( eof ) {
- dprintf("ERROR: Entry for fixedLUafter year not found in %s \n", file_timesteps);
- is_err = true;
- }
-
- if ( ifs_A4xCO2 && bgc_1pctCO2 ) {
- dprintf("ERROR: Both ifs_A4xCO2 && bgc_1pctCO2 selected. Please adjust in config-run.xml\n");
- is_err = true;
- }
- int base_year = CMIP6STARTYEAR;
- if ( ifs_FIXEDYEARCO2 > 0) {
- base_year = ifs_FIXEDYEARCO2;
- }
-
- if ( ifs_A4xCO2 || bgc_1pctCO2 ) {
- fixedNDepafter = base_year;
- fixedLUafter = base_year;
- }
-
- }
- dprintf("fixedNDepafter set to %i \n",fixedNDepafter);
- dprintf("fixedLUafter set to %i \n",fixedLUafter);
- dprintf("End experiment Setup ======= \n");
- fclose(in_ts);
- if ( is_err ) {
- return false;
- }
- else {
- return true;
- }
- }
- 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++){
-
-
- ifsindex++;
- if(nmask[ix]){
- ngridcell++;
- coup_lsmask[0][ix] = nmask[ix];
- EceCoord& c=gridlist.createobj();
- c.id=0;
- c.lon=nlon[ix];
- if(c.lon>=180.)c.lon=c.lon-360.;
- c.lat=nlat[ix];
- if (readcorners) {
-
- }
- c.IFSindex=ifsindex;
- c.area = nareas[ix];
- 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;
-
- }
-
- delete[] nlon;
- delete[] nlat;
- delete[] nmask;
- delete[] nareas;
-
- if (readcorners) {
- delete[] nlon_cor;
- delete[] nlat_cor;
- }
-
- 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) {
-
-
- 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;
-
- 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) {
-
- 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);
-
- 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);
-
- 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);
-
- 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);
- }
-
-
- FILE *ofp;
- if (printgridlist)
- ofp = fopen("ece_gridlist_T255.txt", "w");
-
-
-
- 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;
- }
-
- 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;
-
- if (printgridlist)
- fclose(ofp);
-
-
-
- if (!soilfileOK || !soilOK)
- ok_flag = false;
-
- return ok_flag;
- }
- static void readco2_cmip6(bool& error_flag) {
-
-
-
-
- int status;
- int file_id = 0;
- status = nc_open(file_co2, NC_NOWRITE, &file_id);
-
- if (file_id <= 0) {
- dprintf("!!!ERROR!!! readco_cmip6: could not open CO2 file %s for input\n", file_co2);
- error_flag = true;
- }
-
- int time_id;
- status = nc_inq_unlimdim(file_id, &time_id);
- if (status != NC_NOERR) {
- handle_error(status);
- error_flag = true;
- }
-
- size_t time_size;
- status = nc_inq_dimlen(file_id, time_id, &time_size);
- if (status != NC_NOERR) {
- handle_error(status);
- error_flag = true;
- }
-
- 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;
- }
-
- static size_t start[] = { 0, 0 };
- static size_t count[] = { time_size, 1 };
- float co2_in[NYEAR_CO2];
-
- 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];
- }
-
- status = nc_close(file_id);
- if (status != NC_NOERR) {
- handle_error(status);
- error_flag = true;
- }
- }
- void read_input_data(bool& error_flag) {
-
-
- bool gridOK = read_grid_info(false);
- if (!gridOK) {
- dprintf("LPJ-GUESS: Problem with atmosphere model's mask and grid information. \n");
- error_flag = true;
-
- }
-
- bool rcfileOK = read_ifs_timesteps();
- if (!rcfileOK) {
- dprintf("Problem with EC-Earth .rc file. Quitting. \n");
- error_flag = true;
-
- }
-
-
-
-
-
- readco2_cmip6(error_flag);
-
-
- bool soilOK = read_soil_veg_info(false);
- if (!soilOK) {
- dprintf("Problem with EC-Earth Soil info. \n");
- error_flag = true;
- }
-
-
-
- }
- void call_oasis_get(int timestep){
-
-
-
- if (timestep < 20) {
- dprintf("call_oasis_get on timestep %i - LPJ-GUESS: OASIS communicating\n",timestep);
- }
-
-
- 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){
-
-
-
- int ix;
-
- 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_nat = 0.0;
- float mean_cflux_ant = 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_nat += (float)coup_dcflux_nat[ix]/MAXGRID;
- mean_cflux_ant += (float)coup_dcflux_ant[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_NAT, mean CFLUX_ANT, mean NPP are %12.6f, %12.6f, %12.6f, %12.6f, %12.6f, %12.6f, %12.6f\n",
- timestep,max_llai,max_hlai,mean_llai,mean_hlai,mean_cflux_nat,mean_cflux_ant, mean_npp);
-
- }
-
- 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_nat, coup_dcflux_ant, coup_dnpp);
-
-
- }
- 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;
-
- call_oasis_get(timestep);
-
-
-
-
-
-
-
-
- g=0;
-
- 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++){
- if(coup_lsmask[ix][iy]){
-
- temp_avg += coup_temper[ix];
- prec_avg += coup_precip[ix];
- swrad_avg += coup_sw_radiat[ix];
- lwrad_avg += coup_lw_radiat[ix];
- snowc_avg += coup_snowc[ix];
- snowd_avg += coup_snowd[ix];
- st1l_avg += coup_st1l[ix];
- st2l_avg += coup_st2l[ix];
- st3l_avg += coup_st3l[ix];
- st4l_avg += coup_st4l[ix];
- sm1l_avg += coup_sm1l[ix];
- sm2l_avg += coup_sm2l[ix];
- sm3l_avg += coup_sm3l[ix];
- sm4l_avg += coup_sm4l[ix];
- co2_avg += coup_co2_field[ix];
-
-
-
-
- temp[g]=coup_temper[ix]+TEMP_KTOC;
-
-
- prec[g]=coup_precip[ix]*1.e3*TIMESTEP/1000;
- if (prec[g] < 0.000001) prec[g] = 0.0;
- snowc[g] = coup_snowc[ix];
- snowd[g] = coup_snowd[ix];
- stl[g][0] = coup_st1l[ix]+TEMP_KTOC;
- 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];
- sml[g][1] = coup_sm2l[ix];
- sml[g][2] = coup_sm3l[ix];
- sml[g][3] = coup_sm4l[ix];
-
- swrad[g] = coup_sw_radiat[ix];
- lwrad[g] = coup_lw_radiat[ix];
- if (swrad[g] < 0.000001) swrad[g] = 0.0;
-
- co2_field[g] = coup_co2_field[ix];
- g++;
- }
- }
- }
- if (timestep < 20 && g > 0) {
-
- dprintf("ECEFRAMEWORK: checksum for timestep %d :\n",timestep);
- dprintf("ECEFRAMEWORK: (checksum TEMP) %12.6f\n",temp_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum PREC) %g\n",prec_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum SNOWC) %12.6f\n",snowc_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum SNOWD) %12.6f\n",snowd_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum ST1L) %12.6f\n",st1l_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum ST2L) %12.6f\n",st2l_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum ST3L) %12.6f\n",st3l_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum ST4L) %12.6f\n",st4l_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum SM1L) %12.6f\n",sm1l_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum SM2L) %12.6f\n",sm2l_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum SM3L) %12.6f\n",sm3l_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum SM4L) %12.6f\n",sm4l_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum SW RAD) %12.6f\n", swrad_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum LR RAD) %12.6f\n", lwrad_avg/(double)g);
- dprintf("ECEFRAMEWORK: (checksum CO2) %12.6f\n",co2_avg/(double)g);
- }
-
- return;
- }
- void call_coupler_put(int timestep,double lailow[MAXGRID], double laihigh[MAXGRID],
- int lpjg_typeh[MAXGRID], double lpjg_frach[MAXGRID],
- int lpjg_typel[MAXGRID], double lpjg_fracl[MAXGRID],
- double dcfluxnat[MAXGRID], double dcfluxant[MAXGRID], double dnpp[MAXGRID]) {
-
- int g;
- int ix;
- int iy;
-
-
- g=0;
- for(iy=0;iy<NY_ATMO;iy++) {
- for(ix=0;ix<NX_ATMO;ix++) {
-
- 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_nat[ix]=0.;
- coup_dcflux_ant[ix]=0.;
- coup_dnpp[ix]=0.;
-
- if(coup_lsmask[ix][iy]){
-
- coup_lowlai[ix]=lailow[g];
- coup_highlai[ix]=laihigh[g];
-
- coup_lpjg_typeh[ix]=(double)lpjg_typeh[g];
- coup_lpjg_frach[ix]=lpjg_frach[g];
-
- coup_lpjg_typel[ix]=(double)lpjg_typel[g];
- coup_lpjg_fracl[ix]=lpjg_fracl[g];
-
- coup_dcflux_nat[ix]=dcfluxnat[g];
- coup_dcflux_ant[ix]=dcfluxant[g];
- coup_dnpp[ix]=dnpp[g];
- g++;
- }
- }
- }
-
- call_oasis_put(timestep);
- return;
- }
- 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) {
-
-
- 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;
- }
- void runlpjguess_today(
-
-
- bool islpjgspinup,
-
- int timestep,int isfinal,int yearc,int sim_yr,int monthc,int dayc,int time,
-
- 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],
-
- double lailow[MAXGRID], double laihigh[MAXGRID],
- int lpjg_typeh[MAXGRID], double lpjg_frach[MAXGRID],
- int lpjg_typel[MAXGRID], double lpjg_fracl[MAXGRID],
-
- double co2_ppm[MAXGRID], double dcfluxnat[MAXGRID], double dcfluxant[MAXGRID], double dnpp[MAXGRID]) {
-
- 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;
-
- 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;
-
- float cfluxnattoday=0.0;
- float cfluxanttoday=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;
-
- #ifdef HAVE_MPI
- MPI_Status status;
- #endif
- int rank;
- int num_procs;
-
- hourc=time/3600;
-
-
- rank = get_rank_specific(localcomm);
- num_procs = get_num_local_processes(localcomm);
-
-
- 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);
-
-
-
-
- for (g = rank; g<MAXGRID; g+=num_procs) {
-
- dcfluxnat[g]=0.0;
- dcfluxant[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;
-
-
- }
-
-
-
- double stl1[MAXGRID],stl2[MAXGRID],stl3[MAXGRID],stl4[MAXGRID];
- double sml1[MAXGRID],sml2[MAXGRID],sml3[MAXGRID],sml4[MAXGRID];
- if (rank==0) {
-
- 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];
- }
- }
-
-
- 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
-
-
-
- if (!islpjgspinup) {
- if (timestep<2) dprintf("Broadcasting on node %i\n", rank);
-
- 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
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- 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 };
-
-
-
- for (int gg = rank; gg<MAXGRID; gg+=num_procs) {
-
-
-
- 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);
-
-
-
- cfluxnattoday=0.0;
- cfluxanttoday=0.0;
- npptoday=0.0;
- laiphen_high=0.0;
- laiphen_low=0.0;
-
- EceCoord& c=gridlist[gg];
- alat=(float)c.lat;
- alon=(float)c.lon;
- ifs_soilcd = c.soilcode;
- ifs_index = c.IFSindex;
- ndep_index = c.ndep_index;
- temp_2m=(float)temp[gg];
- temp_soil_1=(float)stl1[gg];
- temp_soil_2=(float)stl2[gg];
- temp_soil_3=(float)stl3[gg];
- temp_soil_4=(float)stl4[gg];
-
- temp_soil = ifs_to_lpjg_soiltemp(temp_soil_2,temp_soil_3);
-
- precip = (float)prec[gg];
-
- co2 = co2_ppm[gg];
-
- ifs_to_lpjg_soilwater(sml1[gg], sml2[gg], sml3[gg], sml4[gg], soilw_surf, soilw_deep);
-
- swrad_net_inst = (float)swrad[gg];
-
- lwrad_net_inst = (float)lwrad[gg];
-
- vegl_cell = c.cvl;
- vegl_type_cell = c.tvl;
- vegh_cell = c.cvh;
- vegh_type_cell = c.tvh;
-
-
-
-
-
-
-
-
-
-
- if (false && timestep < 2) {
-
- 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,
- vegl_cell, vegl_type_cell,
-
- cfluxnattoday, cfluxanttoday, npptoday,
- laiphen_high, laiphen_low,
- ifsvegfrachigh, ifsvegtypehigh, ifsvegfraclow, ifsvegtypelow,
- islpjgspinup, fixedNDepafter, fixedLUafter);
-
- lailow[gg] = laiphen_low;
- laihigh[gg] = laiphen_high;
- if (ifsvegtypehigh == 0)
- laihigh[gg] = 0.0;
- if (ifsvegtypelow == 0)
- lailow[gg] = 0.0;
-
- lpjg_typeh[gg] = (int)ifsvegtypehigh;
- lpjg_frach[gg] = (double)ifsvegfrachigh;
- lpjg_typel[gg] = (int)ifsvegtypelow;
- lpjg_fracl[gg] = (double)ifsvegfraclow;
-
- 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);
- }
-
- dcfluxnat[gg]=cfluxnattoday;
- dcfluxant[gg]=cfluxanttoday;
- dnpp[gg]=npptoday;
-
-
-
- }
-
- #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");
-
-
- MPI_Barrier(localcomm);
-
- if (!islpjgspinup) {
- int tag;
- if (num_procs > 1) {
-
-
-
-
- double lail[MAXGRID];
- double laih[MAXGRID];
- double frh[MAXGRID];
- int tph[MAXGRID];
- double frl[MAXGRID];
- int tpl[MAXGRID];
- double cfl_nat[MAXGRID];
- double cfl_ant[MAXGRID];
- double cfl_npp[MAXGRID];
- if (rank == 0) {
-
- MPI_Status stat;
- for (int pr = 1; pr < num_procs; pr++) {
-
- tag = pr*100+1;
- MPI_Recv(lail,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
-
- tag = pr*100+2;
- MPI_Recv(laih,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
-
- tag = pr*100+3;
- MPI_Recv(tph,MAXGRID,MPI_INT,pr,tag,localcomm,&stat);
-
- tag = pr*100+4;
- MPI_Recv(frh,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
-
- tag = pr*100+5;
- MPI_Recv(tpl,MAXGRID,MPI_INT,pr,tag,localcomm,&stat);
-
- tag = pr*100+6;
- MPI_Recv(frl,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
-
- tag = pr*100+7;
- MPI_Recv(cfl_nat,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
-
- tag = pr*100+8;
- MPI_Recv(cfl_ant,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
-
- tag = pr*100+9;
- MPI_Recv(cfl_npp,MAXGRID,MPI_DOUBLE,pr,tag,localcomm,&stat);
-
- dprintf("Master node %i has received the results from slave process %i.\n", rank, pr);
-
-
-
-
-
-
-
-
-
-
-
-
- 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];
- dcfluxnat[cl] = cfl_nat[cl];
- dcfluxant[cl] = cfl_ant[cl];
- dnpp[cl] = cfl_npp[cl];
- }
- }
- } else {
-
-
- tag = rank*100+1;
- MPI_Send(lailow,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
-
- tag = rank*100+2;
- MPI_Send(laihigh,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
-
- tag = rank*100+3;
- MPI_Send(lpjg_typeh,MAXGRID,MPI_INT,0,tag,localcomm);
-
- tag = rank*100+4;
- MPI_Send(lpjg_frach,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
-
- tag = rank*100+5;
- MPI_Send(lpjg_typel,MAXGRID,MPI_INT,0,tag,localcomm);
-
- tag = rank*100+6;
- MPI_Send(lpjg_fracl,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
-
- tag = rank*100+7;
- MPI_Send(dcfluxnat,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
-
- tag = rank*100+8;
- MPI_Send(dcfluxant,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
-
- tag = rank*100+9;
- MPI_Send(dnpp,MAXGRID,MPI_DOUBLE,0,tag,localcomm);
-
- if (timestep < 2) dprintf("Node %i has sent the results from its cells to the master.\n", rank);
- }
-
- MPI_Barrier(localcomm);
-
- }
- }
- #endif
- }
- void runlpjguess(bool islpjgspinup, bool isparallel, bool error_flag ) {
-
-
- int timestep;
-
-
-
- double temp[MAXGRID],prec[MAXGRID],swrad[MAXGRID],lwrad[MAXGRID];
- double vegl[MAXGRID],vegh[MAXGRID],snowc[MAXGRID],snowd[MAXGRID];
-
- double stl[MAXGRID][NHTESSELSOILLAYERS];
- double sml[MAXGRID][NHTESSELSOILLAYERS];
- int vegl_type[MAXGRID];
-
- double lailow[MAXGRID];
- double laihigh[MAXGRID];
-
-
- int lpjg_typeh[MAXGRID];
- int lpjg_typel[MAXGRID];
-
-
- double lpjg_frach[MAXGRID];
- double lpjg_fracl[MAXGRID];
-
- 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");
- }
-
- double co2_tm5[MAXGRID];
- double co2curr;
- double dcfluxnat[MAXGRID];
- double dcfluxant[MAXGRID];
- double dnpp[MAXGRID];
-
- struct {
- int year;
- int sim_year;
- int month;
- int day;
- int time;
- } eceDate;
- int globalprocs = 1;
- int local_procs = 1;
- int myrank = 0;
- #ifdef HAVE_MPI
- if (islpjgspinup) {
- globalprocs = get_num_processes();
- local_procs = globalprocs;
- myrank = GuessParallel::get_rank();
- }
- #endif
- dprintf("LPJ-GUESS rank %i: %i local process(es), and %i global process(es) = %i \n",myrank, local_procs, globalprocs);
-
-
-
- if(!islpjgspinup && isparallel){
-
- 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);
-
- 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;
- }
-
-
- 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;
- }
-
-
- 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;
- }
-
-
- #ifdef HAVE_MPI
- if (local_procs>1)
- MPI_Barrier(localcomm);
-
- #endif
- } else {
- dprintf("LPJ-GUESS: No OASIS initialisation. \n");
- localcomm = 0;
- couplcomm = 0;
- #ifdef HAVE_MPI
- localcomm = MPI_COMM_WORLD;
- #endif
- }
-
- globalprocs = get_num_global_processes();
- if (!islpjgspinup)
- local_procs = get_num_local_processes(localcomm);
- else
- local_procs = get_num_processes();
-
-
-
- if (isparallel) {
- xtring path;
- path.printf("./run%d", myrank+1);
- dprintf("\nMoving directory to %s on node %d\n",(const char*)path, myrank);
- if (change_directory(path) != 0) {
- fprintf(stderr, "Failed to change to run directory\n");
- return;
- }
- }
-
-
- #ifdef HAVE_MPI
- if (isparallel && !islpjgspinup && local_procs>1)
- MPI_Barrier(localcomm);
- #endif
-
- 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);
-
- eceDate.year=STARTDATE.year;
- eceDate.month=STARTDATE.month;
- eceDate.day=STARTDATE.day;
- eceDate.time=0;
- eceDate.sim_year=0;
-
- int isfinal = 0;
- timestep=0;
-
- int timestepstorun = 0;
- if (islpjgspinup) {
- timestepstorun = 1;
- isfinal = 1;
- } else
- timestepstorun = NTIMESTEP;
-
- while (timestep<timestepstorun) {
-
- if (timestep == NTIMESTEP-1) {
-
- isfinal = 1;
- dprintf("runlpjguess: entered final timestep\n");
- }
-
-
- if (ifs_FIXEDYEARCO2 <= 0 && !ifs_A4xCO2 && !bgc_1pctCO2) {
-
-
- if (eceDate.year<FIRSTYEAR_CO2) {
- co2curr=co2[0];
- }
- else if ( eceDate.year > NYEAR_CO2 ) {
- fail("No CO2 data available beyond %d \n",NYEAR_CO2);
- }
- else {
- co2curr=co2[eceDate.year];
- }
- }
- int base_year = CMIP6STARTYEAR;
-
-
- if ( ifs_FIXEDYEARCO2 > 0 ) {
- co2curr = co2[ifs_FIXEDYEARCO2];
- base_year = ifs_FIXEDYEARCO2;
- }
-
- if ( ifs_A4xCO2 ) {
- if ( eceDate.year >= base_year ) {
- co2curr = 4. * co2[base_year];
- }
- }
-
- if ( bgc_1pctCO2 ) {
- if ( eceDate.year > base_year ) {
- co2curr = pow(1.01,eceDate.year-base_year) * co2[base_year];
- }
- }
- if (timestep < 2) {
- dprintf("fixedYearCO2 set to %d \n",ifs_FIXEDYEARCO2);
- dprintf("ifs_A4xCO2 set to %d \n",ifs_A4xCO2);
- dprintf("bgc_1pctCO2 set to %d \n",bgc_1pctCO2);
- dprintf("fixedLUafter set to %d \n",fixedLUafter);
- dprintf("fixedNDepafter set to %d \n",fixedNDepafter);
- dprintf("runlpjguess: timestep %i of %i (%i s)\n",timestep,NTIMESTEP,timestep*TIMESTEP);
- }
- if (timestep < 5) {
- dprintf("CO2 at timestep %i: %d \n",timestep,co2curr);
- }
-
-
-
- if (printOutputToFile && eceDate.day-1 == 0 && eceDate.month - 1 == 0) {
- for (int ii = 0; ii < MAXGRID; ii++) {
-
- 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;
- }
- }
- }
-
-
- 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);
- }
-
-
- #ifdef HAVE_MPI
- if (!islpjgspinup && local_procs>1)
- MPI_Barrier(localcomm);
- #endif
-
-
- if (!activateTM5coupling)
- for (int jj=0; jj<MAXGRID; jj++) co2_tm5[jj]=co2curr;
-
- if (timestep < 2)
- dprintf("before runlpjguess_today: exchanging fields: timestep %i of %i (%i s)\n", timestep,NTIMESTEP,timestep*TIMESTEP);
-
-
- 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,dcfluxnat,dcfluxant,dnpp);
-
-
- #ifdef HAVE_MPI
- if (!islpjgspinup && local_procs>1)
- MPI_Barrier(localcomm);
- #endif
-
-
- if (!islpjgspinup && myrank==0) {
- call_coupler_put(timestep,lailow,laihigh,lpjg_typeh,lpjg_frach,lpjg_typel,lpjg_fracl,dcfluxnat,dcfluxant,dnpp);
- if (timestep<20) dprintf("runlpjguess: node %i of %i returned from call_coupler_put \n",myrank,local_procs);
- }
-
-
-
- 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);
-
- if (IFLEAPYEARS) {
-
- int yy = eceDate.year;
- if (!(yy%400))
- NDAYMONTH[1]=29;
- else if (!(yy%100))
- NDAYMONTH[1]=28;
- else if (!(yy%4))
- NDAYMONTH[1]=29;
- else
- NDAYMONTH[1]=28;
- }
-
- if (printOutputToFile && isparallel && !islpjgspinup && myrank==0) {
-
-
- for (int ii = 0; ii < MAXGRID; ii++) {
-
-
- 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];
-
-
- if (eceDate.day == 1) {
-
- lpjg_typeh_yesterday[ii] = lpjg_typeh[ii];
- lpjg_typel_yesterday[ii] = lpjg_typel[ii];
- }
-
- 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]);
- }
- }
- fflush(ofp);
- }
- eceDate.time+=TIMESTEP;
- if (eceDate.time>=24*3600) {
- 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;
- else if (!(yy%100))
- NDAYMONTH[1]=28;
- else if (!(yy%4))
- NDAYMONTH[1]=29;
- else
- NDAYMONTH[1]=28;
-
- }
- }
- }
- }
- timestep++;
- }
-
- if (printOutputToFile) {
- delete[] mlailow;
- delete[] mlaihigh;
- delete[] mlpjg_frach;
- delete[] mlpjg_fracl;
- delete[] lpjg_typeh_yesterday;
- delete[] lpjg_typel_yesterday;
- fclose(ofp);
- }
- guess_coupled_finish();
- if(!islpjgspinup){
-
- dprintf("LPJ-GUESS: OASIS terminating...\n");
- OasisCoupler::finalize();
- dprintf("LPJ-GUESS: OASIS terminated!\n");
- }
- dprintf ("Terminating ...\n");
-
- }
-
- int ecemain(const CommandLineArguments& args) {
- dprintf("Running LPJ-GUESS - in ecemain() in eceframework.cpp ...\n");
- bool error_flag = false;
-
-
- 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");
- }
-
- runlpjguess(args.get_islpjgspinup(), args.get_parallel(),error_flag);
- dprintf("LPJ-GUESS: exiting - ecemain() in eceframework.cpp...\n");
- return 0;
- }
|