gdrs.c

Go to the documentation of this file.

/***************************************************
Apache License Version 2.0

Copyright 2007-2011 EcoEquity and Stockholm Environment Institute

Licensed under the Apache License, Version 2.0 (the "License"); you may not use this
file except in compliance with the License. You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software distributed under
the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
KIND, either express or implied. See the License for the specific language governing
permissions and limitations under the License.
****************************************************/

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include "gdrs.h"
#ifdef DMALLOC
#include <dmalloc.h>
#endif


void handle_sql3error(sqlite3 *db, sqlite3_stmt *stmt) {
    fprintf(stderr, "Error: %s\n", sqlite3_errmsg(db));
    sqlite3_finalize(stmt);
    sqlite3_close(db);
}

int loadgdrs(struct gdrs_struct *gdrs, const char *fname) {
    int retcode, sql3ret, i, emergpath_idx;
    struct datablock tempdb, tempdb2;
    time_t now;
    struct tm *timep;
    double val;
    sqlite3 *db;
    sqlite3_stmt *stmt;
    char querybuff[1024], subquery[255], param_name[255];

    // The return value is 0 if the access is permitted, and -1 otherwise
    if (!access(fname, F_OK)) {
        if ((sql3ret = sqlite3_open(fname, &db))) {
            if (db == NULL) {
                return GDRS_NOMEM;
            } else {
                sqlite3_close(db);
                return GDRS_SQLITE_OFFSET + sql3ret;
            }
        }
    } else {
        return GDRS_FREADERROR;
    }

    now = time(NULL);
    timep = localtime(&now);
    gdrs->thisyear = timep->tm_year + 1900;
    
    // *** Get total number of records
    sprintf(querybuff, "SELECT COUNT(iso3) FROM country;");;
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    sqlite3_step(stmt);
    gdrs->ad.numrecs = sqlite3_column_int(stmt, 0);
    sqlite3_finalize(stmt);

    // *** Get min & max years for entire data set
    // These will be updated later as data sets are read in
    sprintf(querybuff, "SELECT MIN(year) AS min_year, MAX(year) AS max_year FROM core;");
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    sqlite3_step(stmt);
    gdrs->ad.minyear = sqlite3_column_int(stmt, 0);
    gdrs->ad.maxyear = sqlite3_column_int(stmt, 1);
    sqlite3_finalize(stmt);

    // *** Get iso3 codes
    sprintf(querybuff, "SELECT iso3 FROM country ORDER BY iso3;");
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    gdrs->iso3 = malloc(gdrs->ad.numrecs * sizeof(char *));
    for (i = 0; sqlite3_step(stmt) == SQLITE_ROW; i++) {
        gdrs->iso3[i] = strdup((const char*)sqlite3_column_text(stmt, 0));
    }
    sqlite3_finalize(stmt);

    // *** Get values for parameters from data table
    sprintf(querybuff, "SELECT param_id, int_val FROM params WHERE int_val IS NOT NULL ORDER BY param_id;");
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    while (sqlite3_step(stmt) == SQLITE_ROW) {
        strcpy(param_name, sqlite3_column_text(stmt, 0));
        if (!strcmp(param_name, "a1_ref_year")) {
            gdrs->a1refyear = sqlite3_column_int(stmt, 1); // The year compared to which A1 reduces emissions (e.g., 1990 or 2005)
        } else if (!strcmp(param_name, "assign_mit_gap_to")) {
             gdrs->assign_mitgap_to = sqlite3_column_int(stmt, 1); // If sequencing: either 1 (Annex 1) or 2 (Annex 2) for who bears responsibility for mitigation gap
        } else if (!strcmp(param_name, "cumsince")) {
            gdrs->cumsince = sqlite3_column_int(stmt, 1); // Calculate responsibility from this year
        } else if (!strcmp(param_name, "do_luxcap")) {
            gdrs->do_luxcap = sqlite3_column_int(stmt, 1); // Calculate responsibility from this year
        } else if (!strcmp(param_name, "emerg_path_id")) {
            gdrs->emerg_path_id = sqlite3_column_int(stmt, 1); // Start year of "emergency progam"
        } else if (!strcmp(param_name, "emergstart")) {
            gdrs->emergstart = sqlite3_column_int(stmt, 1); // Start year of "emergency progam"
        } else if (!strcmp(param_name, "interp_between_thresholds")) {
            gdrs->interp_between_thresholds = sqlite3_column_int(stmt, 1); // Interpolate between thresholds
        } else if (!strcmp(param_name, "kab_only_ratified")) {
            gdrs->kab_only_ratified = sqlite3_column_int(stmt, 1); // KABs--only include countries that ratified Kyoto
        } else if (!strcmp(param_name, "use_lag")) {
            gdrs->use_lag = sqlite3_column_int(stmt, 1); // Apply a (varying) time lag between mitigation investment and realized reductions
        } else if (!strcmp(param_name, "sequenceyear")) {
            gdrs->sequenceyear = sqlite3_column_int(stmt, 1); // End of the first committment period with sequencing: only A1 up to here
        } else if (!strcmp(param_name, "use_kab")) {
            gdrs->use_kab = sqlite3_column_int(stmt, 1); // Use Kyoto-adjusted baselines (KABs)
        } else if (!strcmp(param_name, "use_lulucf")) {
            gdrs->use_lulucf = sqlite3_column_int(stmt, 1); // Flag whether to include land-use change emissions
        } else if (!strcmp(param_name, "use_netexports")) {
            gdrs->use_netexports = sqlite3_column_int(stmt, 1); // Flag whether to include emissions embodied in traded goods
        } else if (!strcmp(param_name, "use_nonco2")) {
            gdrs->use_nonco2 = sqlite3_column_int(stmt, 1); // Flag whether to include non-CO2 gases
        } else if (!strcmp(param_name, "usesequence")) {
            gdrs->usesequence = sqlite3_column_int(stmt, 1); // Boolean: true if Annex 1 countries act first (sequencing)
        }
    }
    sqlite3_finalize(stmt);
    // Set the "cumsince" year as the minyear for the whole analysis -- but not if > 1990, since the 1990 values are used elsewhere
    gdrs->ad.minyear = gdrs->cumsince < 1990 ? gdrs->cumsince : 1990;

    sprintf(querybuff, "SELECT param_id, real_val FROM params WHERE real_val IS NOT NULL ORDER BY param_id;");
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    
    gdrs->billpercgwp = 0; // this will be the sum of the adaptation and mitigation parts
    while (sqlite3_step(stmt) == SQLITE_ROW) {
        strcpy(param_name, sqlite3_column_text(stmt, 0));
        if (!strcmp(param_name, "a1_perc_rdxn")) {
            gdrs->a1_perc_rdxn = sqlite3_column_double(stmt, 1); // If sequencing, % below ref year that A1 must reduce it's emissions by sequence year
        } else if (!strcmp(param_name, "a1_shape_param")) {
            gdrs->a1_shape_param = sqlite3_column_double(stmt, 1); // If sequencing, exponent on tfrac^a1_shape_param that affects how rapidly A1 countries approach target
        } else if (!strcmp(param_name, "billpercgwp_adapt")) {
            gdrs->billpercgwp += sqlite3_column_double(stmt, 1); // Bill as percent of GWP - adaptation
        } else if (!strcmp(param_name, "billpercgwp_mit")) {
            gdrs->billpercgwp += sqlite3_column_double(stmt, 1); // Bill as percent of GWP - mitigation
        } else if (!strcmp(param_name, "dev_thresh")) {
            gdrs->dt_low = sqlite3_column_double(stmt, 1); // Development threshold
        } else if (!strcmp(param_name, "emisselast")) {
            gdrs->emisselast = sqlite3_column_double(stmt, 1); // Emissions elasticity
        } else if (!strcmp(param_name, "lux_thresh")) {
            gdrs->dt_high = sqlite3_column_double(stmt, 1); // Luxury threshold
        } else if (!strcmp(param_name, "lux_thresh_mult")) {
            gdrs->lux_thresh_mult = sqlite3_column_double(stmt, 1); // Responsibility weight
        } else if (!strcmp(param_name, "respweight")) {
            gdrs->respweight = sqlite3_column_double(stmt, 1); // Responsibility weight
        }
    }

    sqlite3_finalize(stmt);
    
    // ************************************************************
    // Emergency pathway
    // ************************************************************
    // First, min & max years
    sprintf(querybuff, "SELECT MIN(year) AS min_year, MAX(year) AS max_year FROM pathways WHERE emergpath_GtC IS NOT NULL;");
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    sqlite3_step(stmt);
    gdrs->emergpath.y_start = sqlite3_column_int(stmt, 0);
    gdrs->emergpath.y_end = sqlite3_column_int(stmt, 1);
    sqlite3_finalize(stmt);
    allocts(&gdrs->emergpath, gdrs->ad);
    emergpath_idx = gdrs->use_lulucf + 2 * gdrs->use_nonco2;
    switch (emergpath_idx) {
    case 0: {
        sprintf(subquery, "source = 'fossil'");
        break;
    }
    case 1: {
        sprintf(subquery, "source = 'fossil' OR source ='lulucf'");
        break;
    }
    case 2: {
        sprintf(subquery, "source = 'fossil' OR source ='nonco2'");
        break;
    }
    default: {
        sprintf(subquery, "source = 'fossil' OR source ='lulucf' OR source ='nonco2'");
    }
    }
    sprintf(querybuff,
            "SELECT year, 1000 * SUM(emergpath_GtC) as emergpath_MtC FROM pathways, pathway_names WHERE pathway_names.pathway_id = %d AND pathway_names.name_short = pathways.pathway AND (%s) GROUP BY year ORDER BY year;",
            gdrs->emerg_path_id, subquery);
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    for (i = 0; sqlite3_step(stmt) == SQLITE_ROW; i++) {
        gdrs->emergpath.data[i] = sqlite3_column_double(stmt, 1);
    }
    sqlite3_finalize(stmt);

    // ************************************************************
    // Fixed factors
    // ************************************************************

    if ((retcode = allocff(&gdrs->annex1, gdrs->ad))) return retcode;
    if ((retcode = allocff(&gdrs->annex2, gdrs->ad))) return retcode;
    if ((retcode = allocff(&gdrs->frozen, gdrs->ad))) return retcode;
    if ((retcode = allocff(&gdrs->kyoto_startep_frac_1990, gdrs->ad))) return retcode;
    if ((retcode = allocff(&gdrs->kyoto_ratified, gdrs->ad))) return retcode;
    if ((retcode = allocff(&gdrs->rci_lagged, gdrs->ad))) return retcode;

    sprintf(querybuff, "SELECT iso3, value FROM flags WHERE flag=\"annex_1\" ORDER BY iso3;");
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    for (i = 0; i < gdrs->ad.numrecs; i++) {
        sqlite3_step(stmt);
        val = sqlite3_column_int(stmt, 1);
        gdrs->annex1.data[i] = (double) val;
    }
    sqlite3_finalize(stmt);

    sprintf(querybuff, "SELECT iso3, value FROM flags WHERE flag=\"annex_2\" ORDER BY iso3;");
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    for (i = 0; i < gdrs->ad.numrecs; i++) {
        sqlite3_step(stmt);
        val = sqlite3_column_int(stmt, 1);
        gdrs->annex2.data[i] = (double) val;
    }
    sqlite3_finalize(stmt);
    
    // Note that SQLite3 returns 0.0 for NULL floats and 0 for NULL ints
    sprintf(querybuff, "SELECT country.iso3, commitment_percent FROM country LEFT OUTER JOIN kyoto_info ON country.iso3 = kyoto_info.iso3 ORDER BY country.iso3;");
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    for (i = 0; i < gdrs->ad.numrecs; i++) {
        sqlite3_step(stmt);
        val = sqlite3_column_double(stmt, 1);
        gdrs->kyoto_startep_frac_1990.data[i] = 0.01 * val;
    }
    sqlite3_finalize(stmt);

    sprintf(querybuff, "SELECT country.iso3, ratified FROM country LEFT OUTER JOIN kyoto_info ON country.iso3 = kyoto_info.iso3 ORDER BY country.iso3;");
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    for (i = 0; i < gdrs->ad.numrecs; i++) {
        sqlite3_step(stmt);
        val = sqlite3_column_int(stmt, 1);
        gdrs->kyoto_ratified.data[i] = (double) val;
    }
    sqlite3_finalize(stmt);

    // ************************************************************
    // Data files
    // ************************************************************

    // *** Population
    if ((retcode = readtable(db, "pop_person", &gdrs->pop, gdrs->ad))) {
        sqlite3_close(db);
        return retcode;
    }
    if (gdrs->pop.y_start > gdrs->ad.minyear) {
        gdrs->ad.minyear = gdrs->pop.y_start;
    }
    if (gdrs->pop.y_end < gdrs->ad.maxyear) {
        gdrs->ad.maxyear = gdrs->pop.y_end;
    }
    // writedb2delim(gdrs->pop, "testing/pop_temp.csv", ',', gdrs->ad);

    // *** GDP
    if ((retcode = readtable(db, "gdp_blnUSDMER", &gdrs->gdp, gdrs->ad))) {
        sqlite3_close(db);
        return retcode;
    }
    if (gdrs->gdp.y_start > gdrs->ad.minyear) {
        gdrs->ad.minyear = gdrs->gdp.y_start;
    }
    if (gdrs->gdp.y_end < gdrs->ad.maxyear) {
        gdrs->ad.maxyear = gdrs->gdp.y_end;
    }
    // writedb2delim(gdrs->gdp, "testing/gdp_temp.csv", ',', gdrs->ad);

    // *** PPP to MER ratio
    if ((retcode = readtable(db, "ppp2mer", &gdrs->ppp2mer, gdrs->ad))) {
        sqlite3_close(db);
        return retcode;
    }
    if (gdrs->ppp2mer.y_start > gdrs->ad.minyear) {
        gdrs->ad.minyear = gdrs->ppp2mer.y_start;
    }
    if (gdrs->ppp2mer.y_end < gdrs->ad.maxyear) {
        gdrs->ad.maxyear = gdrs->ppp2mer.y_end;
    }
    // writedb2delim(gdrs->ppp2mer, "testing/ppp2mer_temp.csv", ',', gdrs->ad);

    // *** Gini coefficient
    if ((retcode = readtable(db, "gini", &gdrs->gini, gdrs->ad))) {
        sqlite3_close(db);
        return retcode;
    }
    if (gdrs->gini.y_start > gdrs->ad.minyear) {
        gdrs->ad.minyear = gdrs->gini.y_start;
    }
    if (gdrs->gini.y_end < gdrs->ad.maxyear) {
        gdrs->ad.maxyear = gdrs->gini.y_end;
    }
    // writedb2delim(gdrs->gini, "testing/gini.csv", ',', gdrs->ad);

    // *** Voluntary reductions
    if ((retcode = readtable(db, "vol_rdxn_MtC", &gdrs->volrdxn, gdrs->ad))) {
        sqlite3_close(db);
        return retcode;
    }
    if (gdrs->usesequence) {
        if (gdrs->volrdxn.y_start > gdrs->ad.minyear) {
            gdrs->ad.minyear = gdrs->volrdxn.y_start;
        }
        if (gdrs->volrdxn.y_end < gdrs->ad.maxyear) {
            gdrs->ad.maxyear = gdrs->volrdxn.y_end;
        }
        // writedb2delim(gdrs->volrdxn, "testing/volrdxn.csv", ',', gdrs->ad);
    }

    // *** Baseline fossil CO2 emissions trajectory
    if ((retcode = readtable(db, "fossil_CO2_MtC", &gdrs->bl_fossil, gdrs->ad))) {
        sqlite3_close(db);
        return retcode;
    }
    if (gdrs->bl_fossil.y_start > gdrs->ad.minyear) {
        gdrs->ad.minyear = gdrs->bl_fossil.y_start;
    }
    if (gdrs->bl_fossil.y_end < gdrs->ad.maxyear) {
        gdrs->ad.maxyear = gdrs->bl_fossil.y_end;
    }
    // Make a copy for "baseline"
    gdrs->baseline = scalmult(1, gdrs->bl_fossil, gdrs->ad);
    // writedb2delim(gdrs->bl_fossil, "testing/bl_fossil.csv", ',', gdrs->ad);

    // *** Baseline land use, land-use change & forestry trajectory
    if ((retcode = readtable(db, "LULCF_MtC", &gdrs->bl_lulucf, gdrs->ad))) {
        sqlite3_close(db);
        return retcode;
    }
    if (gdrs->use_lulucf) {
        if (gdrs->bl_lulucf.y_start > gdrs->ad.minyear) {
            gdrs->ad.minyear = gdrs->bl_lulucf.y_start;
        }
        if (gdrs->bl_lulucf.y_end < gdrs->ad.maxyear) {
            gdrs->ad.maxyear = gdrs->bl_lulucf.y_end;
        }
        tempdb = binop(gdrs->baseline, gdrs->bl_lulucf, '+', gdrs->ad);
        cleardata(gdrs->baseline);
        gdrs->baseline = scalmult(1, tempdb, gdrs->ad);
        cleardata(tempdb);
    }
    // writedb2delim(gdrs->use_lulucf, "testing/use_lulucf.csv", ',', gdrs->ad);

    // *** Baseline non-CO2 emissions
    if ((retcode = readtable(db, "NonCO2_MtCe", &gdrs->bl_nonco2, gdrs->ad))) {
        sqlite3_close(db);
        return retcode;
    }
    if (gdrs->use_nonco2) {
        if (gdrs->bl_nonco2.y_start > gdrs->ad.minyear) {
            gdrs->ad.minyear = gdrs->bl_nonco2.y_start;
        }
        if (gdrs->bl_nonco2.y_end < gdrs->ad.maxyear) {
            gdrs->ad.maxyear = gdrs->bl_nonco2.y_end;
        }
        tempdb = binop(gdrs->baseline, gdrs->bl_nonco2, '+', gdrs->ad);
        cleardata(gdrs->baseline);
        gdrs->baseline = scalmult(1, tempdb, gdrs->ad);
        cleardata(tempdb);
    }
    // writedb2delim(gdrs->bl_nonco2, "testing/bl_nonco2.csv", ',', gdrs->ad);
    // writedb2delim(gdrs->baseline, "testing/baseline.csv", ',', gdrs->ad);

    // *** Baseline embodied emissions
    // TODO: Change code on this and similar so that text for "use_netexports" is done
    //    before this, then put the check on years into readtable
    if ((retcode = readtable(db, "net_export_ktCO2", &gdrs->bl_netexports, gdrs->ad))) {
        sqlite3_close(db);
        return retcode;
    }
    if (gdrs->use_netexports) {
        if (gdrs->bl_netexports.y_start > gdrs->ad.minyear) {
            gdrs->ad.minyear = gdrs->bl_netexports.y_start;
        }
        if (gdrs->bl_netexports.y_end < gdrs->ad.maxyear) {
            gdrs->ad.maxyear = gdrs->bl_netexports.y_end;
        }
        // Convert from ktCO2 to MtC
        tempdb2 = scalmult(1.0e-3 * 3.0/11.0, gdrs->bl_netexports, gdrs->ad);
        // Subtract netexports to get corrected value
        tempdb = binop(gdrs->baseline, tempdb2, '-', gdrs->ad);
        cleardata(gdrs->baseline);
        gdrs->baseline = scalmult(1, tempdb, gdrs->ad);
        cleardata(tempdb);
        cleardata(tempdb2);
    }
    
    // *** Tax levels
    sprintf(querybuff, "SELECT seq_no, ppp, value, quintile, global FROM tax_levels;");
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    // Have to loop through twice; first get the number of tax table entries
    for (gdrs->num_tax_levels = 0; sqlite3_step(stmt) == SQLITE_ROW; gdrs->num_tax_levels++);
    if (gdrs->num_tax_levels > TAX_MAX_ENTRIES) {
        return GDRS_TOO_MANY_TAX_ENTRIES;
    }
    sqlite3_reset(stmt);
    for (i = 0; sqlite3_step(stmt) == SQLITE_ROW; i++) {
        gdrs->tax_levels[i].seq_no = sqlite3_column_int(stmt, 0);
        gdrs->tax_levels[i].is_ppp = sqlite3_column_int(stmt, 1);
        if (sqlite3_column_type(stmt, 2) != SQLITE_NULL) {
            gdrs->tax_levels[i].value = sqlite3_column_double(stmt, 2);
            gdrs->tax_levels[i].type = TAX_TYPE_VALUE;
        } else if (sqlite3_column_type(stmt, 3) != SQLITE_NULL) {
            gdrs->tax_levels[i].quintile = sqlite3_column_double(stmt, 3);
            if (sqlite3_column_int(stmt, 4)) {
                gdrs->tax_levels[i].type = TAX_TYPE_GLOBAL_QUINTILE;
            } else {
                gdrs->tax_levels[i].type = TAX_TYPE_NATL_QUINTILE;
            }
        }
    }
    sqlite3_finalize(stmt);
    
    //
    // Done with the database--can close
    //
    sqlite3_close(db);

    // *** Create scaled population
    gdrs->scalpop = scalmult(1.0e-9, gdrs->pop, gdrs->ad);
    // writedb2delim(gdrs->scalpop, "testing/scalpop.csv", ',', gdrs->ad);

    // *** Calculate sd of log (sigma) corresponding to Gini
    gdrs->sdlog = transform(gdrs->gini, *gini2sdlog, gdrs->ad);
    // writedb2delim(gdrs->sdlog, "testing/sdlog.csv", ',', gdrs->ad);
    
    if ((retcode = allocff(&gdrs->kyoto_startep_frac, gdrs->ad))) return retcode;
    for (i = 0; i < gdrs->ad.numrecs; i++) {
        val = getval(gdrs->baseline, i, 1990, gdrs->ad)/getval(gdrs->baseline, i, gdrs->emergstart, gdrs->ad);
        val *= gdrs->kyoto_startep_frac_1990.data[i];
        gdrs->kyoto_startep_frac.data[i] = val > 1.0 ? 1.0 : val;
    }
    // writeff2delim(gdrs->kyoto_startep_frac, "testing/kab_value.csv", ',', gdrs->ad);

    // *** Allocate memory for calculated results
    gdrs->r.y_start = gdrs->ad.minyear;
    gdrs->c.y_start = gdrs->ad.minyear;
    gdrs->rci.y_start = gdrs->ad.minyear;
    gdrs->alloc_gdrs.y_start = gdrs->ad.minyear;
    gdrs->pop_above_dl.y_start = gdrs->ad.minyear;
    gdrs->a1domrdxn.y_start = gdrs->ad.minyear;
    gdrs->lux_emiss.y_start = gdrs->ad.minyear;
    gdrs->zu_adj.y_start = gdrs->ad.minyear;
    gdrs->yu_adj.y_start = gdrs->ad.minyear;
    gdrs->pop_above_lux.y_start = gdrs->ad.minyear;
    gdrs->lux_emiss_applied.y_start = gdrs->ad.minyear;
    gdrs->kab_gap.y_start = gdrs->ad.minyear;

    gdrs->r.y_end = gdrs->ad.maxyear;
    gdrs->c.y_end = gdrs->ad.maxyear;
    gdrs->rci.y_end = gdrs->ad.maxyear;
    gdrs->alloc_gdrs.y_end = gdrs->ad.maxyear;
    gdrs->pop_above_dl.y_end = gdrs->ad.maxyear;
    gdrs->a1domrdxn.y_end = gdrs->ad.maxyear;
    gdrs->lux_emiss.y_end = gdrs->ad.maxyear;
    gdrs->zu_adj.y_end = gdrs->ad.maxyear;
    gdrs->yu_adj.y_end = gdrs->ad.maxyear;
    gdrs->pop_above_lux.y_end = gdrs->ad.maxyear;
    gdrs->lux_emiss_applied.y_end = gdrs->ad.maxyear;
    gdrs->kab_gap.y_end = gdrs->ad.maxyear;

    if ((retcode = allocdata(&gdrs->r, gdrs->ad))) return retcode;
    if ((retcode = allocdata(&gdrs->c, gdrs->ad))) return retcode;
    if ((retcode = allocdata(&gdrs->rci, gdrs->ad))) return retcode;
    if ((retcode = allocdata(&gdrs->alloc_gdrs, gdrs->ad))) return retcode;
    if ((retcode = allocdata(&gdrs->pop_above_dl, gdrs->ad))) return retcode;
    if ((retcode = allocdata(&gdrs->a1domrdxn, gdrs->ad))) return retcode;
    if ((retcode = allocdata(&gdrs->lux_emiss, gdrs->ad))) return retcode;
    if ((retcode = allocdata(&gdrs->zu_adj, gdrs->ad))) return retcode;
    if ((retcode = allocts(&gdrs->yu_adj, gdrs->ad))) return retcode; // yu_adj is a timeseries
    if ((retcode = allocdata(&gdrs->pop_above_lux, gdrs->ad))) return retcode;
    if ((retcode = allocdata(&gdrs->lux_emiss_applied, gdrs->ad))) return retcode;
    if ((retcode = allocdata(&gdrs->kab_gap, gdrs->ad))) return retcode;
    
    for (i = 0; i < gdrs->num_tax_levels; i++) {
        gdrs->tax_revenue_mer_dens[i].y_start = gdrs->ad.minyear;
        gdrs->tax_revenue_mer_dens[i].y_end = gdrs->ad.maxyear;
        if ((retcode = allocdata(&(gdrs->tax_revenue_mer_dens[i]), gdrs->ad))) return retcode;
        gdrs->tax_revenue_ppp_dens[i].y_start = gdrs->ad.minyear;
        gdrs->tax_revenue_ppp_dens[i].y_end = gdrs->ad.maxyear;
        if ((retcode = allocdata(&(gdrs->tax_revenue_ppp_dens[i]), gdrs->ad))) return retcode;
        gdrs->tax_income_mer_dens[i].y_start = gdrs->ad.minyear;
        gdrs->tax_income_mer_dens[i].y_end = gdrs->ad.maxyear;
        if ((retcode = allocdata(&(gdrs->tax_income_mer_dens[i]), gdrs->ad))) return retcode;
        gdrs->tax_income_ppp_dens[i].y_start = gdrs->ad.minyear;
        gdrs->tax_income_ppp_dens[i].y_end = gdrs->ad.maxyear;
        if ((retcode = allocdata(&(gdrs->tax_income_ppp_dens[i]), gdrs->ad))) return retcode;
        gdrs->tax_pop_dens[i].y_start = gdrs->ad.minyear;
        gdrs->tax_pop_dens[i].y_end = gdrs->ad.maxyear;
        if ((retcode = allocdata(&(gdrs->tax_pop_dens[i]), gdrs->ad))) return retcode;
        gdrs->tax_pop_below[i].y_start = gdrs->ad.minyear;
        gdrs->tax_pop_below[i].y_end = gdrs->ad.maxyear;
        if ((retcode = allocdata(&(gdrs->tax_pop_below[i]), gdrs->ad))) return retcode;
    }

    calcybar(gdrs);

    // Create the "z-values" for thresholds
    gdrs->zl = z(gdrs->dt_low, gdrs, XR_PPP);
    gdrs->zu = z(gdrs->dt_high, gdrs, XR_MER);

    calcallocs(gdrs);
    
    return GDRS_OK;
}

int readtable(sqlite3 *db, const char *column, struct datablock *table, struct appdata ad) {
    sqlite3_stmt *stmt;
    int i, year, sql3ret, err;
    double val;
    char querybuff[1024];

    // First, min & max years
    sprintf(querybuff, "SELECT MIN(year) AS min_year, MAX(year) AS max_year FROM core WHERE %s IS NOT NULL;", column);
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    sqlite3_step(stmt);
    table->y_start = sqlite3_column_int(stmt, 0);
    table->y_end = sqlite3_column_int(stmt, 1);
    sqlite3_finalize(stmt);

    // Initialize the variable
    if ((err = allocdata(table, ad))) {
        return err;
    }

    // Next, whole table
    sprintf(querybuff, "SELECT iso3, year, %s FROM core WHERE %s IS NOT NULL ORDER BY iso3, year;", column, column);
    if ((sql3ret = sqlite3_prepare_v2(db, querybuff, strlen(querybuff), &stmt, NULL))) {
        handle_sql3error(db, stmt);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    for (i = 0; i < ad.numrecs; i++) {
        for (year = table->y_start; year <= table->y_end; year++) {
            sqlite3_step(stmt);
            val = sqlite3_column_double(stmt, 2);
            setval(*table, val, i, year, ad);
        }
    }
    sqlite3_finalize(stmt);

    return GDRS_OK;
}

void cleargdrs(struct gdrs_struct gdrs) {
    int i;

    cleardata(gdrs.bl_netexports);
    cleardata(gdrs.bl_fossil);
    cleardata(gdrs.bl_lulucf);
    cleardata(gdrs.bl_nonco2);
    cleardata(gdrs.baseline);
    cleardata(gdrs.gdp);
    cleardata(gdrs.pop);
    cleardata(gdrs.scalpop);
    cleardata(gdrs.ppp2mer);
    cleardata(gdrs.ybar_ppp);
    cleardata(gdrs.gini);
    cleardata(gdrs.sdlog);
    cleardata(gdrs.r);
    cleardata(gdrs.c);
    cleardata(gdrs.rci);
    cleardata(gdrs.alloc_gdrs);
    cleardata(gdrs.pop_above_dl);
    cleardata(gdrs.volrdxn);
    cleardata(gdrs.a1domrdxn);
    cleardata(gdrs.lux_emiss);
    cleardata(gdrs.zu_adj);
    cleardata(gdrs.pop_above_lux);
    cleardata(gdrs.lux_emiss_applied);
    cleardata(gdrs.kab_gap);
    cleardata(gdrs.zu);
    cleardata(gdrs.zl);

    clearts(gdrs.yu_adj);
    clearts(gdrs.emergpath);

    clearff(gdrs.annex1);
    clearff(gdrs.annex2);
    clearff(gdrs.frozen);
    clearff(gdrs.rci_lagged);

    for (i = 0; i < gdrs.ad.numrecs; i++) {
        free(gdrs.iso3[i]);
    }
    free(gdrs.iso3);
    
    for (i = 0; i < gdrs.num_tax_levels; i++) {
        cleardata(gdrs.tax_revenue_mer_dens[i]);
        cleardata(gdrs.tax_revenue_ppp_dens[i]);
        cleardata(gdrs.tax_income_mer_dens[i]);
        cleardata(gdrs.tax_income_ppp_dens[i]);
        cleardata(gdrs.tax_pop_dens[i]);
        cleardata(gdrs.tax_pop_below[i]);
    }
}

double globaldist(double y, int year, struct gdrs_struct *gdrs, int ppp_or_mer, double *dist_deriv) {
    int i;
    double ybar, P, sigma, z;
    double cum, deriv, cumP;
    
    if (y < 0) {
        return NAN;
    } else if (y == 0) {
        return 0.0;
    } else if (isinf(y) == 1) {
        return 1.0;
    }
    cum = 0; // Cumulative distribution up to specified income level
    deriv = 0; // Derivative of the cumulative distribution
    cumP = 0;
    for (i = 0; i < gdrs->ad.numrecs; i++) {
        // Store values from tables into local variables for convenience
        ybar = getval(gdrs->ybar_ppp, i, year, gdrs->ad);
        sigma = getval(gdrs->sdlog, i, year, gdrs->ad);
        P = getval(gdrs->scalpop, i, year, gdrs->ad);
        if (ppp_or_mer == XR_MER) {
            ybar /= getval(gdrs->ppp2mer, i, year, gdrs->ad);
        }
        z = ztransf(y, ybar, sigma);
        
        // Note that dz/dy = 1/(sigma * y)
        if (dist_deriv) {
            deriv += P * dnorm(z)/(sigma * y);
        }
        cum += P * pnorm(z);
        cumP += P;
    }
    cum /= cumP;
    deriv /= cumP;
    
    if (dist_deriv) *dist_deriv = deriv;
    
    return cum;
}

double globalquintile(double q, int year, struct gdrs_struct *gdrs, int ppp_or_mer) {
    double y, qcalc, qderiv;
    int niter;
    
    y = INIT_Y;
    
    for (niter = 0; niter < MAXITER; niter++) {
        qcalc = globaldist(y, year, gdrs, ppp_or_mer, &qderiv);
        if (fabs(qcalc - q) < EPS_NEWTON) break;
        
        y += (q - qcalc)/qderiv;
    }
    
    if (niter == MAXITER) {
        y = NAN;
    }
    
    return y;
}

double nationalquintile(double q, double ybar, double sigma) {
    return ybar * exp(sigma * (qnorm(q) - 0.5 * sigma));
}

void combine_gap(struct fixedfactor* gap, int i, double newgap, int combine) {
    switch (combine) {
        case GAP_COMBINE_NONE:
            gap->data[i] = newgap;
            break;
        case GAP_COMBINE_ADD:
            gap->data[i] += newgap;
            break;
        case GAP_COMBINE_MAX:
            if (newgap > gap->data[i]) gap->data[i] = newgap;
            break;
        case GAP_COMBINE_MIN:
            if (newgap < gap->data[i]) gap->data[i] = newgap;
            break;
        default:
            break;
    }
}

void calc_kab(struct fixedfactor* gap, double* totgap, int year, int period, int combine, struct gdrs_struct *gdrs) {
    int i, emergstart;
    double bau_at_start_ep, bau_now, tot_bau, ep;
    double kyoto_startep_frac, filter;
    double totgap_temp, newgap, ep_gap;
    struct fixedfactor newgaps;
    
    
    if (period == PER_BEFORE_EP) {
        *totgap = 0.0;
        for (i = 0; i < gdrs->ad.numrecs; i++) {
            setval(gdrs->kab_gap, 0.0, i, year, gdrs->ad);
            if (combine == GAP_COMBINE_NONE) {
                // Make sure initialize to zero, if this is the first one
                combine_gap(gap, i, 0.0, combine);
            }
        }
    } else {
        allocff(&newgaps, gdrs->ad);

        ep = gettsval(gdrs->emergpath, year, gdrs->ad);
        emergstart = gdrs->emergstart;

        totgap_temp = 0;
        tot_bau = 0;
        for (i = 0; i < gdrs->ad.numrecs; i++) {
            bau_at_start_ep = getval(gdrs->baseline, i, emergstart, gdrs->ad);
            kyoto_startep_frac = gdrs->kyoto_startep_frac.data[i];
            if (gdrs->kab_only_ratified) {
                filter = gdrs->kyoto_ratified.data[i];
            } else {
                if (gdrs->kyoto_ratified.data[i] > 0.5) {
                    filter = gdrs->kyoto_ratified.data[i];
                } else {
                    // Any countries without a Kyoto commitment get assigned a value of 0.0
                    filter = fabs(kyoto_startep_frac) < EPSILON ? 0.0 : 1.0;
                }
            }
            tot_bau += bau_now = getval(gdrs->baseline, i, year, gdrs->ad);
            newgap = filter * bau_at_start_ep * (1.0 - kyoto_startep_frac); // If "bau_now", the slab expands; if "bau_at_start_ep" it's a fixed slab
            if (newgap < 0.0) newgap = 0.0;
            newgaps.data[i] = newgap;
            totgap_temp += newgap;
        }

        // Check if the correction exceeds the current gap between BAU & EP, and adjust if necessary
        *totgap = 0;
        ep_gap = tot_bau > ep ? tot_bau - ep : 0.0;
        for (i = 0; i < gdrs->ad.numrecs; i++) {
            if (totgap_temp > ep_gap) {
                newgap = fabs(totgap_temp) > EPSILON ? newgaps.data[i] * ep_gap/totgap_temp : 0.0;
            } else {
                newgap = newgaps.data[i];
            }
            setval(gdrs->kab_gap, newgap, i, year, gdrs->ad);
            if (gdrs->use_kab) {
                combine_gap(gap, i, newgap, combine);
                *totgap += newgap;
            }
        }

        clearff(newgaps);
    }
    
}


void calc_luxemiss(struct fixedfactor* gap, double* totgap, int year, int period, int combine, struct gdrs_struct *gdrs) {
    int i;
    double lux_emiss, ybar, yu, yu_adj, yu_ppp, f1;
    double blval, gamma, sigma, zu, filter, newgap;
    
    gamma = gdrs->emisselast;
    yu = gdrs->dt_high; // The user-defined, fixed luxury threshold
    if (gdrs->do_luxcap) {
        yu_adj = gettsval(gdrs->yu_adj, year, gdrs->ad);
    } else {
        yu_adj = yu;
    }

    *totgap = 0;
    for (i = 0; i < gdrs->ad.numrecs; i++) {
        // Always report luxury emissions (not net of frozen emissions), whether capping baselines or not, using user-defined line
        blval = getval(gdrs->baseline, i, year, gdrs->ad);
        sigma = getval(gdrs->sdlog, i, year, gdrs->ad);
        // Fixed line...
        zu = getval(gdrs->zu, i, year, gdrs->ad);
        ybar = getval(gdrs->ybar_ppp, i, year, gdrs->ad);
        // Convert fixed lux threshold to PPP for consistent comparison to ybar
        yu_ppp = yu * getval(gdrs->ppp2mer, i, year, gdrs->ad);
        f1 = pow(yu_ppp/ybar, gamma) * exp(-0.5 * gamma * (gamma - 1) * sigma * sigma);
        lux_emiss = blval * (1 - pnorm(zu - gamma * sigma) - f1 * (1 - pnorm(zu)));
        setval(gdrs->lux_emiss, lux_emiss, i, year, gdrs->ad);
        if (gdrs->do_luxcap && (period != PER_BEFORE_EP)) {
            // Correct baseline value for "frozen" emissions
            blval -= gdrs->frozen.data[i];
            zu = getval(gdrs->zu_adj, i, year, gdrs->ad);
            // This captures the PPP conversion factor:
            //   f1' = f1 * (yu_adj * ppp2mer/yu * ppp2mer)^gamma = f1 * (yu_adj/yu)^gamma
            f1 *= pow(yu_adj/yu, gamma);
            // When sequencing, only apply luxury emissions to Annex 1
            if (gdrs->usesequence && period == PER_SEQUENCE) {
                filter = gdrs->annex1.data[i];
            } else {
                filter = 1;
            }
            newgap = filter * blval * (1 - pnorm(zu - gamma * sigma) - f1 * (1 - pnorm(zu)));
        } else {
            newgap = 0;
        }
        combine_gap(gap, i, newgap, combine);
        *totgap += gap->data[i];
        setval(gdrs->lux_emiss_applied, gap->data[i], i, year, gdrs->ad);
    }
}

void calcallocs(struct gdrs_struct *gdrs) {
    int i, year;
    // Basic RCI/GDRs
    double bltot, globalrdxn, alloc, rci;
    int lag_start, lag_ndx;
    // For adjusted baselines
    double totgap;
    struct fixedfactor gap;
    // For sequencing
    int period;
    double a1, a2, a1rci, a2rci, a1natlrdxn;
    double tfrac, a1rdxn;
    double blref, blval, f;

    //--------------------------------------
    //
    // Initialize variables
    //
    //--------------------------------------
    
    // Sequencing
    f = 0.01 * gdrs->a1_perc_rdxn;
    for (i = 0; i < gdrs->ad.numrecs; i++) {
        gdrs->frozen.data[i] = 0.0;
    }

    // To store results from "gaps", such as luxury-capped baselines or Kyoto-adjusted baselines
    allocff(&gap, gdrs->ad);
    
    // Find start year for mitigation lag/smoothing
    lag_start = gdrs->thisyear > gdrs->cumsince ? gdrs->thisyear : gdrs->cumsince;

    //--------------------------------------
    //
    // Loop over years
    //
    //--------------------------------------
    for (year = gdrs->ad.minyear; year <= gdrs->ad.maxyear; year++) {
        // First, get RCI based on allocations from previous (or initial) year; this will also calculate yu_adj if doing luxury-capping
        nextrci(year, gdrs);
        
        if (gdrs->use_lag && year > lag_start) {
            lag_ndx = year - lag_start;
        } else {
            lag_ndx = 0;
        }
        
        for (i = 0; i < gdrs->ad.numrecs; i++) {
            rci = getval(gdrs->rci, i, year, gdrs->ad);
            if (gdrs->use_lag && lag_ndx > 0) {
                rci = (rci + lag_ndx * getffval(gdrs->rci_lagged, i, gdrs->ad))/(lag_ndx + 1);
                setval(gdrs->rci, rci, i, year, gdrs->ad);
            }
            setffval(gdrs->rci_lagged, rci, i, gdrs->ad);
        }
        
        // Figure out the sequencing "period" -- period can equal PER_SEQUENCE only if using sequencing
        if (year < gdrs->emergstart) {
            period = PER_BEFORE_EP;
        } else {
            if (gdrs->usesequence && year <= gdrs->sequenceyear) {
                period = PER_SEQUENCE;
            } else {
                period = PER_GDRS;
            }
        }

        //------ Prep: if sequencing
        a1rci = a2rci = a1rdxn = 0.0;
        tfrac = -999.0; // Unused if not doing sequencing, so initialize to something wrong but recognizable
        if (period == PER_SEQUENCE) {
            tfrac = (double)(year - gdrs->emergstart) / (double)(gdrs->sequenceyear - gdrs->emergstart);
            for (i = 0; i < gdrs->ad.numrecs; i++) {
                rci = getffval(gdrs->rci_lagged, i, gdrs->ad);
                a1 = gdrs->annex1.data[i];
                a2 = gdrs->annex2.data[i];
                a1rci += a1 * rci;
                a2rci += a2 * rci;
                // a1 is either 0.0 or 1.0. This will definitely select for 1.0 (Annex 1)
                if (a1 > 0.5) {
                    // TODO: This could give an error if a1refyear is less than start year for baseline
                    blref = getval(gdrs->baseline, i, gdrs->a1refyear, gdrs->ad);
                    // Annex 1 domestic reduction commitment -- interpolate between emerg start and end of 1st commitment period
                    a1rdxn += (getval(gdrs->baseline, i, year, gdrs->ad) - (1 - f) * blref) * pow(tfrac, gdrs->a1_shape_param);
                }
            }
        }

        //------ Calculate luxury emissions and (if doing luxury-capped baselines), update gap & totgap
        calc_luxemiss(&gap, &totgap, year, period, GAP_COMBINE_NONE, gdrs);
        calc_kab(&gap, &totgap, year, period, GAP_COMBINE_MAX, gdrs);

        //----------------------------------
        // Calculate allocations
        //----------------------------------
        if (period == PER_BEFORE_EP) {
            globalrdxn = 0.0;
        } else {
            bltot = 0.0;
            for (i = 0; i < gdrs->ad.numrecs; i++) {
                // Note: frozen emissions are zero if not sequencing
                bltot += getval(gdrs->baseline, i, year, gdrs->ad) - gdrs->frozen.data[i];
            }
            // Note: totgap is zero if not using adjusted baselines
            globalrdxn = bltot - totgap - gettsval(gdrs->emergpath, year, gdrs->ad);
        }
        for (i = 0; i < gdrs->ad.numrecs; i++) {
            rci = getffval(gdrs->rci_lagged, i, gdrs->ad);
            blval = getval(gdrs->baseline, i, year, gdrs->ad) - gdrs->frozen.data[i];

            //------ Sequencing
            a1natlrdxn = 0.0;
            if (period == PER_SEQUENCE) {
                a1 = gdrs->annex1.data[i];
                a2 = gdrs->annex2.data[i];
                // Non-Annex 1
                if (a1 < 0.5) {
                    alloc = blval;
                } else {
                    // First, subtract off share of A1 domestic reductions
                    a1natlrdxn = (rci/a1rci) * a1rdxn;
                    alloc = blval - gap.data[i] - a1natlrdxn;
                    // Next, assign remaining:
                    if (gdrs->assign_mitgap_to == ANNEX1) {
                        alloc -= (rci/a1rci) * (globalrdxn - a1rdxn);
                    } else {
                        alloc -= a2 * (rci/a2rci) * (globalrdxn - a1rdxn);
                    }
                }
                // Store "frozen obligations" if at end of sequencing period
                if (year == gdrs->sequenceyear && gdrs->usesequence) {
                    if (a1 > 0.5) {
                        gdrs->frozen.data[i] = blval - alloc;
                    }
                }
            } else {
                // Note: gap = 0 if not using adjusted baselines
                // Note: frozen emissions have already been subtracted off of baseline emissions
                alloc = blval - gap.data[i] - globalrdxn * rci;
            }
            // Update data series
            setval(gdrs->a1domrdxn, a1natlrdxn, i, year, gdrs->ad);
            setval(gdrs->alloc_gdrs, alloc, i, year, gdrs->ad);
        }
    }
    clearff(gap);
}

void nextrci(int year, struct gdrs_struct *gdrs) {
    int i, j, yprev, do_interp, ppp_or_mer;
    double sigma, gamma, popdens, z_tax, zl, zu, yl, yu, yu_ppp, ybar, E, P, A;
    double c, r, rci, rci_indiv, ppp2mer, pi, tax_mer;
    double csum, rsum, rprev;
    double gwp, tot_tax, ytax, ytax_ppp;
    struct fixedfactor r_ann;

    // TODO: Check return value and abort if nonzero
    allocff(&r_ann, gdrs->ad);
    
    do_interp = gdrs->interp_between_thresholds;
    gamma = gdrs->emisselast;
    yl = gdrs->dt_low;

    if (gdrs->do_luxcap) {
        yu = luxcap(year, gdrs);
    } else {
        yu = gdrs->dt_high;
    }

    yprev = year - 1;
    if (yprev < gdrs->baseline.y_start) yprev = gdrs->baseline.y_start;

    csum = rsum = gwp = 0.0;
    for (i = 0; i < gdrs->ad.numrecs; i++) {
        // Store values from tables into local variables for convenience
        ybar = getval(gdrs->ybar_ppp, i, year, gdrs->ad);
        sigma = getval(gdrs->sdlog, i, year, gdrs->ad);
        P = getval(gdrs->scalpop, i, year, gdrs->ad);
        ppp2mer = getval(gdrs->ppp2mer, i, year, gdrs->ad);
        pi = 1/ppp2mer;
        // Calculate gross world product in MER terms
        gwp += P * pi * ybar;
        // The luxury threshold is in MER terms, but some calculations here are in PPP,
        // so offer scaled version for each country
        yu_ppp = yu/pi;
        zl = getval(gdrs->zl, i, year, gdrs->ad);
        if (gdrs->do_luxcap) {
            zu = getval(gdrs->zu_adj, i, year, gdrs->ad);
        } else {
            zu = getval(gdrs->zu, i, year, gdrs->ad);
        }
        if (year <= gdrs->cumsince) {
            E = getval(gdrs->baseline, i, yprev, gdrs->ad);
        } else {
            E = getval(gdrs->alloc_gdrs, i, yprev, gdrs->ad);
        }
        A = rhat(E, gamma, sigma, ybar);
        // Calculate capacity and annual responsibility
        if (!do_interp) {
            c = P * capbetween(-INFINITY, INFINITY, zl, zl, yl, yl, ybar, sigma, pi, gdrs->lux_thresh_mult);
            r_ann.data[i] = respbetween(gamma, -INFINITY, INFINITY, zl, zl, yl, yl, ybar, sigma, A);
        } else {
            c = P * capbetween(-INFINITY, INFINITY, zl, zu, yl, yu_ppp, ybar, sigma, pi, gdrs->lux_thresh_mult);
            r_ann.data[i] = respbetween(gamma, -INFINITY, INFINITY, zl, zu, yl, yu_ppp, ybar, sigma, A);
        }
        // Responsibility is cumulative
        if (year <= gdrs->cumsince) {
            rprev = 0.0;
        } else {
            rprev = getval(gdrs->r, i, yprev, gdrs->ad);
        }
        r = r_ann.data[i] + rprev;
        // Don't allow negative values
        if (r < 0) r = 0;
        // Update global totals
        csum += c;
        rsum += r;
        // Store the non-normalized c & r values (these will be overwritten later) and pop above dl
        setval(gdrs->c, c, i, year, gdrs->ad);
        setval(gdrs->r, r, i, year, gdrs->ad);
        setval(gdrs->pop_above_dl, P * (1 - pnorm(zl)), i, year, gdrs->ad);
        // For population above luxury threshold, use the fixed, not adjusted, value
        zu = getval(gdrs->zu, i, year, gdrs->ad);
        setval(gdrs->pop_above_lux, P * (1 - pnorm(zu)), i, year, gdrs->ad);
    }
    // For tax, calculate any global quintiles for the current year
    for (i = 0; i < gdrs->num_tax_levels; i++) {
        if (gdrs->tax_levels[i].type == TAX_TYPE_GLOBAL_QUINTILE) {
            if (gdrs->tax_levels[i].is_ppp) {
                ppp_or_mer = XR_PPP;
            } else {
                ppp_or_mer = XR_MER;
            }
            gdrs->tax_levels[i].value = globalquintile(gdrs->tax_levels[i].quintile, year, gdrs, ppp_or_mer);
        }
    }
    // Calculate RCI and tax levels
    tot_tax = 0.01 * gdrs->billpercgwp * gwp;
    for (i = 0; i < gdrs->ad.numrecs; i++) {
        c = getval(gdrs->c, i, year, gdrs->ad);
        r = getval(gdrs->r, i, year, gdrs->ad);
        rci = (1 - gdrs->respweight) * c/csum + gdrs->respweight * r/rsum;
        setval(gdrs->rci, rci, i, year, gdrs->ad);
        
        // Tax calculation
        P = getval(gdrs->scalpop, i, year, gdrs->ad);
        ybar = getval(gdrs->ybar_ppp, i, year, gdrs->ad);
        sigma = getval(gdrs->sdlog, i, year, gdrs->ad);
        ppp2mer = getval(gdrs->ppp2mer, i, year, gdrs->ad);
        pi = 1/ppp2mer;
        if (year <= gdrs->cumsince) {
            E = getval(gdrs->baseline, i, yprev, gdrs->ad);
        } else {
            E = getval(gdrs->alloc_gdrs, i, yprev, gdrs->ad);
        }
        A = rhat(E, gamma, sigma, ybar);
        for (j = 0; j < gdrs->num_tax_levels; j++) {
            if (gdrs->tax_levels[j].type == TAX_TYPE_NATL_QUINTILE) {
                gdrs->tax_levels[j].value = nationalquintile(gdrs->tax_levels[i].quintile, ybar, sigma);
                if (!gdrs->tax_levels[j].is_ppp) {
                    gdrs->tax_levels[j].value *= pi;
                }
            }
            if (gdrs->tax_levels[j].is_ppp) {
                ytax_ppp = gdrs->tax_levels[j].value;
                ytax = pi * ytax_ppp;
            } else {
                ytax = gdrs->tax_levels[j].value;
                ytax_ppp = ytax/pi;
            }
            rci_indiv = (1 - gdrs->respweight) * P * c_indiv(ytax, pi * yl, yu, do_interp)/csum;
            rci_indiv += gdrs->respweight * (r/rsum) * r_indiv(gamma, ytax_ppp, yl, yu/pi, A, do_interp)/r_ann.data[i];
            // TODO: Make this work properly: shouldn't get NAN's
            if (isnan(rci_indiv)) {
                rci_indiv = -9999;
            }
            // These densities can be summed across countries and then divided by each other. So:
            //   tax rate = tax_revenue_mer_dens/tax_income_mer_dens
            //   tax per capita = tax_revenue_mer_dens/tax_pop_dens
            //   income at tax level (MER) = tax_income_mer_dens/tax_pop_dens
            //   income at tax level (PPP) = tax_income_ppp_dens/tax_pop_dens
            // etc.
            // Note that ybar is in PPP terms
            z_tax = ztransf(ytax_ppp, ybar, sigma);
            popdens = P * dnorm(z_tax)/(sigma * ytax_ppp);
            tax_mer = rci_indiv * tot_tax * popdens;
            setval(gdrs->tax_revenue_mer_dens[j], tax_mer, i, year, gdrs->ad);
            setval(gdrs->tax_revenue_ppp_dens[j], tax_mer/pi, i, year, gdrs->ad);
            setval(gdrs->tax_income_mer_dens[j], ytax * popdens, i, year, gdrs->ad);
            setval(gdrs->tax_income_ppp_dens[j], ytax_ppp * popdens, i, year, gdrs->ad);
            setval(gdrs->tax_pop_dens[j], popdens, i, year, gdrs->ad);
            setval(gdrs->tax_pop_below[j], P * pnorm(z_tax), i, year, gdrs->ad);
        }
    }
    
    clearff(r_ann);
}

double M(double a, double zc, double ybar, double sigma) {
    double fact1, fact2, fact3;

    fact1 = pow(ybar, a);
    fact2 = exp(0.5 * sigma * sigma * a * (a - 1));
    if (!isinf(zc)) {
        fact3 = 1 - pnorm(zc - a * sigma);
    } else if (zc > 0) {
        fact3 = 0;
    } else {
        fact3 = 1;
    }

    return fact1 * fact2 * fact3;
}

double a_gamma(double gamma, double y, double a, double b) {
    double r, t1, t2;
    
    if (fabs(a - b) < EPSILON) {
        return 0;
    }
    
    r = gamma/(1 + gamma);
    t1 = r * (pow(y,gamma + 1) - pow(a, gamma+1));
    t2 = a * (pow(y, gamma) - pow(a, gamma));
    
    return (t1 - t2)/(b - a);
}

double b_gamma(double gamma, double y, double b) {
    return pow(y, gamma) - pow(b, gamma);
}

double A_gamma(double gamma, double za, double zb, double yl, double yu, double ybar, double sigma) {
    double r, t1, t2, t3, z;
    
    if (fabs(yu - yl) < EPSILON) {
        return 0;
    }
    
    r = 1/(1 + gamma);
    t1 = -gamma * r * (M(gamma + 1, zb, ybar, sigma) - M(gamma + 1, za, ybar, sigma));
    t2 = yl * (M(gamma, zb, ybar, sigma) - M(gamma, za, ybar, sigma));
    t3 = -pow(yl, gamma + 1) * r * (M(0, zb, ybar, sigma) - M(0, za, ybar, sigma));
    
    return (t1 + t2 + t3)/(yu - yl);
}

double B_gamma(double gamma, double za, double zb, double yu, double ybar, double sigma) {
    double t1, t2;
    
    t1 = -(M(gamma, zb, ybar, sigma) - M(gamma, za, ybar, sigma));
    t2 = -pow(yu, gamma) * (M(0, zb, ybar, sigma) - M(0, za, ybar, sigma));
    
    return t1 - t2;
}

double capbetween(double za, double zb, double zl, double zu, double yl, double yu_ppp, double ybar, double sigma, double pi, double luxmult) {
    double retval;
    double gamma = 1;
    
    if (zb <= za) {
        return NAN;
    } else if (zb < zl) {
        return 0;
    } else if (za < zl) {
        if (zb < zu) {
            retval = A_gamma(gamma, zl, zb, yl, yu_ppp, ybar, sigma);
        } else {
            retval = A_gamma(gamma, zl, zu, yl, yu_ppp, ybar, sigma) + 
                    a_gamma(gamma, yu_ppp, yl, yu_ppp) * (pnorm(zb) - pnorm(zu)) +
                    luxmult * B_gamma(gamma, zu, zb, yu_ppp, ybar, sigma);
        }
    } else if (za < zu) {
        if (zb < zu) {
            retval = A_gamma(gamma, za, zb, yl, yu_ppp, ybar, sigma);
        } else {
            retval = A_gamma(gamma, za, zu, yl, yu_ppp, ybar, sigma) +
                    a_gamma(gamma, yu_ppp, yl, yu_ppp) * (pnorm(zb) - pnorm(zu)) +
                    luxmult * B_gamma(gamma, zu, zb, yu_ppp, ybar, sigma);
        }
    } else {
        retval = luxmult * B_gamma(gamma, za, zb, yu_ppp, ybar, sigma) +
                a_gamma(gamma, yu_ppp, yl, yu_ppp) * (pnorm(zb) - pnorm(za));
    }
    
    return pi * retval;
}

double respbetween(double gamma, double za, double zb, double zl, double zu, double yl, double yu_ppp, double ybar, double sigma, double A) {
    double retval;
    
    if (zb <= za) {
        return NAN;
    } else if (zb < zl) {
        return 0;
    } else if (za < zl) {
        if (zb < zu) {
            retval = A_gamma(gamma, zl, zb, yl, yu_ppp, ybar, sigma);
        } else {
            retval = A_gamma(gamma, zl, zu, yl, yu_ppp, ybar, sigma) + 
                    a_gamma(gamma, yu_ppp, yl, yu_ppp) * (pnorm(zb) - pnorm(zu)) +
                    B_gamma(gamma, zu, zb, yu_ppp, ybar, sigma);
        }
    } else if (za < zu) {
        if (zb < zu) {
            retval = A_gamma(gamma, za, zb, yl, yu_ppp, ybar, sigma);
        } else {
            retval = A_gamma(gamma, za, zu, yl, yu_ppp, ybar, sigma) + 
                    a_gamma(gamma, yu_ppp, yl, yu_ppp) * (pnorm(zb) - pnorm(zu)) +
                    B_gamma(gamma, zu, zb, yu_ppp, ybar, sigma);
        }
    } else {
        retval = B_gamma(gamma, za, zb, yu_ppp, ybar, sigma) +
                a_gamma(gamma, yu_ppp, yl, yu_ppp) * (pnorm(zb) - pnorm(za));
    }
    
    return A * retval;
}

double rhat(double E, double gamma, double sigma, double ybar) {
    return E * exp(-0.5 * gamma * (gamma-1) * sigma * sigma)/pow(ybar, gamma);
}

double c_indiv(double y, double yl, double yu, int do_interp) {
    if (y < yl) {
        return 0;
    } else if (do_interp) {
        if (y < yu)  {
            return a_gamma(1, y, yl, yu);
        } else {
            return a_gamma(1, yu, yl, yu) + b_gamma(1, y, yu);
        }
    } else {
        return y - yl;
    }
}

double r_indiv(double gamma, double y, double yl, double yu, double A, int do_interp) {
    double retval;
    
    if (y < yl) {
        retval = 0;
    } else if (do_interp) {
        if (y < yu)  {
            retval = a_gamma(gamma, y, yl, yu);
        } else {
            retval = a_gamma(gamma, yu, yl, yu) + b_gamma(gamma, y, yu);
        }
    } else {
        retval = pow(y,gamma) - pow(yl,gamma);
    }
    
    return A * retval;
}

double luxcap(int year, struct gdrs_struct *gdrs) {
    double sigma, gamma, zu, zushift, yu, ybar, E;
    double Ebau, Egap, Eep, Eabove, Fderiv;
    double f1, f2, f3, f4, f5, frozen_tot;
    int niter, i, recalc;
    // Sequencing
    int period;
    struct fixedfactor filter;

    // Assume line is not recalculated unless told otherwise
    recalc = 0;

    // The defined line & corresponding zu
    yu = gdrs->dt_high;

    gamma = gdrs->emisselast;
    // Emergency pathway
    Eep = gettsval(gdrs->emergpath, year, gdrs->ad);

    // Create a filter -- either A1, if sequencing and in the right period, or all 1s
    allocff(&filter, gdrs->ad);

    if (year < gdrs->emergstart) {
        period = PER_BEFORE_EP;
    } else {
        if (gdrs->usesequence && year <= gdrs->sequenceyear) {
            period = PER_SEQUENCE;
        } else {
            period = PER_GDRS;
        }
    }

    // Before emergency start date, the gap is (supposed to be) equal to zero, and
    // no work is done
    if (period != PER_BEFORE_EP) {
        Ebau = 0;
        Eabove = 0;
        frozen_tot = 0;
        for (i = 0; i < gdrs->ad.numrecs; i++) {
            // Set filter
            if (period == PER_SEQUENCE) {
                filter.data[i] = gdrs->annex1.data[i];;
            } else {
                filter.data[i] = 1;
            }
            // Calculate emissions above user-defined luxury threshold & global BAU emissions
            E = getval(gdrs->baseline, i, year, gdrs->ad);
            // Subtract off any "frozen" emissions from sequencing
            Ebau += E - gdrs->frozen.data[i];
            zu = getval(gdrs->zu, i, year, gdrs->ad);
             // Initialize the adjusted zu to the default value
            setval(gdrs->zu_adj, zu, i, year, gdrs->ad);
            // Add to Eabove only if the country is participating in GDRs
            if (filter.data[i]) {
                sigma = getval(gdrs->sdlog, i, year, gdrs->ad);
                ybar = getval(gdrs->ybar_ppp, i, year, gdrs->ad);
                // The luxury threshold is in MER
                ybar /= getval(gdrs->ppp2mer, i, year, gdrs->ad);
                // Only one factor, here, but have f1-f5 defined, so use f1
                f1 = pow(yu/ybar, gamma) * exp(-0.5 * gamma * (gamma - 1) * sigma * sigma);
                Eabove += E * (1 - pnorm(zu - gamma * sigma) - f1 * (1 - pnorm(zu)));
            }
            frozen_tot += gdrs->frozen.data[i];
        }

        // Calculate the gap that must be closed between BAU & EP
        Egap = Ebau + frozen_tot < Eep ? 0 : Ebau + frozen_tot - Eep;
        // If the luxury capping doesn't close the gap, then more work to do: return default

        if (Eabove > Egap) {
            recalc = 1;
            // Iterate to convergence, adjusting yu
            for (niter = 0; niter < MAXITER; niter++) {
                Eabove = 0;
                Fderiv = 0;
                for (i = 0; i < gdrs->ad.numrecs; i++) {
                    // Only summing over those countries that are participating in GDRs
                    if (!filter.data[i]) continue;
                    // Store values from tables into local variables for convenience
                    ybar = getval(gdrs->ybar_ppp, i, year, gdrs->ad);
                    // The luxury threshold is in MER
                    ybar /= getval(gdrs->ppp2mer, i, year, gdrs->ad);
                    sigma = getval(gdrs->sdlog, i, year, gdrs->ad);
                    // Have to recalculate zu for the altered yu
                    zu = (1/sigma) * log(yu/ybar) + 0.5 * sigma;
                    E = getval(gdrs->baseline, i, year, gdrs->ad) - gdrs->frozen.data[i];

                    f1 = pow(yu/ybar, gamma) * exp(-0.5 * gamma * (gamma - 1) * sigma * sigma);
                    zushift = zu - gamma * sigma;
                    Eabove += E * (1 - pnorm(zushift) - f1 * (1 - pnorm(zu)));
                    f2 = E/(sigma * yu);
                    f3 = dnorm(zushift);
                    f4 = gamma * sigma * (1 - pnorm(zu));
                    f5 = dnorm(zu);
                    Fderiv += f2 * (f3 + f1 * (f4 - f5));
                }
                if (fabs(Egap - Eabove) < EPS_NEWTON) {
                    break;
                }
                // Note that Egap - Ebove = (Ebau - Eep) - (Ebau - Elcbaseline) = Elcbaseline - Eep
                //    This is the same as the function f(yu) in the technical documentation
                yu = yu - (Egap - Eabove)/Fderiv;
            }

            if (niter == MAXITER) {
                yu = NAN;
            }
        }
    }

    // Assign adjusted zu's for this year, recalculating if yu has changed
    for (i = 0; i < gdrs->ad.numrecs; i++) {
        if (yu == NAN) {
            zu = NAN;
        } else {
            if (recalc) {
                ybar = getval(gdrs->ybar_ppp, i, year, gdrs->ad);
                // The luxury threshold is in MER
                ybar /= getval(gdrs->ppp2mer, i, year, gdrs->ad);
                sigma = getval(gdrs->sdlog, i, year, gdrs->ad);
                // Have to recalculate zu for the altered yu
                zu = (1/sigma) * log(yu/ybar) + 0.5 * sigma;
            } else {
                zu = getval(gdrs->zu, i, year, gdrs->ad);
            }
        }
        setval(gdrs->zu_adj, zu, i, year, gdrs->ad);
    }

    settsval(gdrs->yu_adj, yu, year, gdrs->ad);

    clearff(filter);

    return yu;
}

void calcybar(struct gdrs_struct *gdrs) {
    int i, year;
    double a[3], y;

    a[0] = gdrs->ppp2mer.y_start;
    a[1] = gdrs->gdp.y_start;
    a[2] = gdrs->scalpop.y_start;
    gdrs->ybar_ppp.y_start = a[0];
    for (i = 1; i < 3; i++) {
        if (a[i] > gdrs->ybar_ppp.y_start) gdrs->ybar_ppp.y_start = a[i];
    }
    a[0] = gdrs->ppp2mer.y_end;
    a[1] = gdrs->gdp.y_end;
    a[2] = gdrs->scalpop.y_end;
    gdrs->ybar_ppp.y_end = a[0];
    for (i = 1; i < 3; i++) {
        if (a[i] < gdrs->ybar_ppp.y_end) gdrs->ybar_ppp.y_end = a[i];
    }

    allocdata(&gdrs->ybar_ppp, gdrs->ad);

    for (i = 0; i < gdrs->ad.numrecs; i++) {
        for (year = gdrs->ybar_ppp.y_start; year <= gdrs->ybar_ppp.y_end; year++) {
            y = getval(gdrs->ppp2mer, i, year, gdrs->ad) * getval(gdrs->gdp, i, year, gdrs->ad);
            y /= getval(gdrs->scalpop, i, year, gdrs->ad);
            setval(gdrs->ybar_ppp, y, i, year, gdrs->ad);
        }
    }
}

struct datablock z(double y, struct gdrs_struct *gdrs, int ppp_or_mer) {
    int i, year;
    double ybar, sigma, z;
    struct datablock result;
    int do_calc = 1;
    
    if (y == 0) {
        z = -INFINITY;
        do_calc = 0;
    }
    if (isinf(y) == 1) {
        z = INFINITY;
        do_calc = 0;
    }

    // Use most constraining bounds
    result.y_start = gdrs->ybar_ppp.y_start;
    if (gdrs->sdlog.y_start > result.y_start) result.y_start = gdrs->sdlog.y_start;
    result.y_end = gdrs->ybar_ppp.y_end;
    if (gdrs->sdlog.y_end < result.y_end) result.y_end = gdrs->sdlog.y_end;
    allocdata(&result, gdrs->ad);

    for (i = 0; i < gdrs->ad.numrecs; i++) {
        for (year = result.y_start; year <= result.y_end; year++) {
            if (do_calc) {
                ybar = getval(gdrs->ybar_ppp, i, year, gdrs->ad);
                if (ppp_or_mer == XR_MER) {
                    ybar /= getval(gdrs->ppp2mer, i, year, gdrs->ad);
                }
                sigma = getval(gdrs->sdlog, i, year, gdrs->ad);
                z = (1/sigma) * log(y/ybar) + 0.5 * sigma;
            }
            setval(result, z, i, year, gdrs->ad);
        }
    }
    return result;
}

int dump2sql3(const char *fname, struct gdrs_struct gdrs) {
    int sql3ret, i, j, k, year, numvals;
    char valbuff[GDRS_SQLITE_HALFBUFSIZE], sqlbuff[GDRS_SQLITE_BUFSIZE], *errmsg;
    sqlite3 *db;
    sqlite3_stmt *stmt;
    double alloc, tmp;

    if ((sql3ret = sqlite3_open(fname, &db))) {
        if (db == NULL) {
            return GDRS_NOMEM;
        } else {
            sqlite3_close(db);
            return GDRS_SQLITE_OFFSET + sql3ret;
        }
    }
    
    sqlite3_exec(db, "PRAGMA journal_mode = MEMORY", NULL, NULL, &errmsg);

    if ((sql3ret = sqlite3_exec(db, "DROP TABLE IF EXISTS gdrs;", NULL, NULL, &errmsg))) {
        fprintf(stderr, "Error: %s\n", errmsg);
        sqlite3_free(errmsg);
        sqlite3_close(db);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }
    
    strcpy(sqlbuff, "CREATE TABLE gdrs (");
    strcat(sqlbuff, "iso3 TEXT NOT NULL");
    strcat(sqlbuff, ", year INT NOT NULL");
    strcat(sqlbuff, ", allocation_MtC REAL");
    strcat(sqlbuff, ", a1_dom_rdxn_MtC REAL");
    strcat(sqlbuff, ", r_MtC REAL");
    strcat(sqlbuff, ", c_blnUSDMER REAL");
    strcat(sqlbuff, ", rci REAL");
    strcat(sqlbuff, ", pop_above_dl REAL");
    strcat(sqlbuff, ", pop_above_lux REAL");
    strcat(sqlbuff, ", lux_emiss_MtC REAL");
    strcat(sqlbuff, ", lux_emiss_applied_MtC REAL");
    strcat(sqlbuff, ", kyoto_gap_MtC REAL");
    for (i = 0; i < gdrs.num_tax_levels; i++) {
        // Can't count on a POSIX environment, so don't use positional notation (eg., %1$d) in format
        sprintf(valbuff, ", tax_pop_below_%d REAL, tax_income_mer_dens_%d REAL, tax_income_ppp_dens_%d REAL, tax_revenue_mer_dens_%d REAL, tax_revenue_ppp_dens_%d REAL, tax_pop_dens_%d REAL",
                gdrs.tax_levels[i].seq_no, gdrs.tax_levels[i].seq_no, gdrs.tax_levels[i].seq_no, gdrs.tax_levels[i].seq_no, gdrs.tax_levels[i].seq_no, gdrs.tax_levels[i].seq_no);
        strcat(sqlbuff, valbuff);
    }
    strcat(sqlbuff, ");");
    
    if ((sql3ret = sqlite3_exec(db, sqlbuff, NULL, NULL, &errmsg))) {
        fprintf(stderr, "Error: %s\n", errmsg);
        sqlite3_free(errmsg);
        sqlite3_close(db);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }

    
    numvals = GDRS_NUM_PARAMS + GDRS_NUM_TAX_PARAMS * gdrs.num_tax_levels;
    sprintf(sqlbuff, "INSERT INTO gdrs VALUES(");
    for (i = 0; i < numvals; i++) {
        if (!i) {
            strcat(sqlbuff, "?");
        } else {
            strcat(sqlbuff, ",?");
        }
    }
    strcat(sqlbuff, ")");
    sqlite3_prepare_v2(db, sqlbuff, GDRS_SQLITE_BUFSIZE, &stmt, 0);
    
    sqlite3_exec(db, "BEGIN;", NULL, NULL, &errmsg);
    sqlite3_free(errmsg);
    for (i = 0; i < gdrs.ad.numrecs; i++) {
        for (year = gdrs.rci.y_start; year <= gdrs.rci.y_end; year++) {
            // If correcting for embodied emissions, get that correction to shift baseline back
            alloc = getval(gdrs.alloc_gdrs, i, year, gdrs.ad);
            if (gdrs.use_netexports) {
                alloc += 1.0e-3 * (3.0/11.0) * getval(gdrs.bl_netexports, i, year, gdrs.ad);
            }
            sqlite3_bind_text(stmt, 1, gdrs.iso3[i], -1, SQLITE_TRANSIENT);
            sqlite3_bind_int(stmt, 2, year);
            sqlite3_bind_double(stmt, 3, alloc);
            sqlite3_bind_double(stmt, 4, getval(gdrs.a1domrdxn, i, year, gdrs.ad));
            sqlite3_bind_double(stmt, 5, getval(gdrs.r, i, year, gdrs.ad));
            sqlite3_bind_double(stmt, 6, getval(gdrs.c, i, year, gdrs.ad));
            sqlite3_bind_double(stmt, 7, getval(gdrs.rci, i, year, gdrs.ad));
            sqlite3_bind_double(stmt, 8, 1e9 * getval(gdrs.pop_above_dl, i, year, gdrs.ad));
            sqlite3_bind_double(stmt, 9, 1e9 * getval(gdrs.pop_above_lux, i, year, gdrs.ad));
            sqlite3_bind_double(stmt, 10, getval(gdrs.lux_emiss, i, year, gdrs.ad));
            sqlite3_bind_double(stmt, 11, getval(gdrs.lux_emiss_applied, i, year, gdrs.ad));
            sqlite3_bind_double(stmt, 12, getval(gdrs.kab_gap, i, year, gdrs.ad));
            for (j = 0; j < gdrs.num_tax_levels; j++) {
                k = GDRS_NUM_PARAMS + GDRS_NUM_TAX_PARAMS * j;
                sqlite3_bind_double(stmt, 1 + k, 1e9 * getval(gdrs.tax_pop_below[j], i, year, gdrs.ad));
                sqlite3_bind_double(stmt, 2 + k, getval(gdrs.tax_income_mer_dens[j], i, year, gdrs.ad));
                sqlite3_bind_double(stmt, 3 + k, getval(gdrs.tax_income_ppp_dens[j], i, year, gdrs.ad));
                sqlite3_bind_double(stmt, 4 + k, getval(gdrs.tax_revenue_mer_dens[j], i, year, gdrs.ad));
                sqlite3_bind_double(stmt, 5 + k, getval(gdrs.tax_revenue_ppp_dens[j], i, year, gdrs.ad));
                sqlite3_bind_double(stmt, 6 + k, getval(gdrs.tax_pop_dens[j], i, year, gdrs.ad));
            }
            if (sqlite3_step(stmt) != SQLITE_DONE) {
                fprintf(stderr, "Error: %s\n\tSQL: %s\n", errmsg, sqlbuff);
                sqlite3_free(errmsg);
                sqlite3_close(db);
                return GDRS_SQLITE_OFFSET + sql3ret;
            }
            sqlite3_reset(stmt);
        }
    }
    sqlite3_exec(db, "END;", NULL, NULL, &errmsg);
    sqlite3_free(errmsg);

    if ((sql3ret = sqlite3_exec(db, "CREATE UNIQUE INDEX idx_gdrs ON gdrs(iso3 ASC, year ASC);", NULL, NULL, &errmsg))) {
        fprintf(stderr, "Error: %s\n", errmsg);
        sqlite3_free(errmsg);
        sqlite3_close(db);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }

    // Add the scheme
    if ((sql3ret = sqlite3_exec(db,
                                "INSERT OR IGNORE INTO schemes VALUES (\"Greenhouse Development Rights\", \"gdrs\");",
                                NULL, NULL, &errmsg))) {
        fprintf(stderr, "Error: %s\n", errmsg);
        sqlite3_free(errmsg);
        sqlite3_close(db);
        return GDRS_SQLITE_OFFSET + sql3ret;
    }

    sqlite3_close(db);

    return GDRS_OK;
}