import numpy as np
import pandas as pd
import os, sys
import re

MOL_WEIGHTS = {
    "CO2": 44.01,
    "C": 12.01,
}

rcmip_output_dir = "../RCMIP3_input_datafiles/"

def check_co2_component_sums(df, rtol=1e-6, atol=1e-8):
    """Check that Emissions|CO2 equals AFOLU + Energy components per Scenario.

    For each Scenario in ``df`` this helper looks for rows named
    ``Emissions|CO2``, ``Emissions|CO2|AFOLU`` and
    ``Emissions|CO2|Energy and Industrial Processes``. When a CO2 total row
    exists and at least one of the two component rows is present, the function
    sums the component rows (summing multiple rows if present) across all
    four-digit year columns and compares the result to the CO2 total.

    Parameters
    ----------
    df : pandas.DataFrame
        DataFrame containing at minimum ``Scenario`` and ``Variable`` columns
        and year columns named with four-digit strings (e.g. '1750').
    rtol, atol : float, optional
        Relative and absolute tolerance passed to :func:`numpy.allclose` when
        comparing the summed components to the CO2 total.

    Returns
    -------
    dict
        Mapping from scenario -> dict with keys ``ok`` (bool) and ``mismatches``
        (list of (year, co2, parts_sum, diff) tuples) for any years that
        differ beyond the tolerances. Scenarios where no check was applicable
        (no CO2 row or no component rows) are omitted from the result.
    """
    if df is None:
        raise ValueError("df is None")

    # detect column names for Scenario and Variable (case-insensitive)
    cols_lower = {c.lower(): c for c in df.columns}
    scen_col = cols_lower.get('scenario')
    var_col = cols_lower.get('variable')
    if scen_col is None or var_col is None:
        raise ValueError("Input DataFrame must contain 'Scenario' and 'Variable' columns")

    # year columns are 4-digit names
    year_cols = [c for c in df.columns if isinstance(c, str) and re.fullmatch(r"\d{4}", c)]
    if not year_cols:
        raise ValueError("No year columns found (4-digit names expected)")

    results = {}
    for scenario, group in df.groupby(scen_col):
        # locate rows
        co2_rows = group[group[var_col] == 'Emissions|CO2']
        afolu_rows = group[group[var_col] == 'Emissions|CO2|AFOLU']
        energy_rows = group[group[var_col] == 'Emissions|CO2|Energy and Industrial Processes']

        if co2_rows.empty:
            # nothing to check for this scenario
            continue
        if afolu_rows.empty and energy_rows.empty:
            # no component rows present; nothing to check
            continue

        # sum across possible multiple matching rows
        co2_sum = co2_rows[year_cols].astype(float).sum(axis=0).values
        parts_sum = np.zeros_like(co2_sum, dtype=float)
        if not afolu_rows.empty:
            parts_sum = parts_sum + afolu_rows[year_cols].astype(float).sum(axis=0).values
        if not energy_rows.empty:
            parts_sum = parts_sum + energy_rows[year_cols].astype(float).sum(axis=0).values

        # compare
        equal_mask = np.isclose(co2_sum, parts_sum, rtol=rtol, atol=atol, equal_nan=True)
        mismatches = []
        for col, ok, cval, pval in zip(year_cols, equal_mask, co2_sum, parts_sum):
            if not ok:
                mismatches.append((col, float(cval), float(pval), float(cval - pval)))

        results[scenario] = {
            'ok': len(mismatches) == 0,
            'mismatches': mismatches,
        }

    # Print a short summary
    for scen, info in results.items():
        if info['ok']:
            continue
            #print(f"Scenario {scen}: OK (CO2 equals AFOLU+Energy)")
        else:
            print(f"Scenario {scen}: MISMATCH in {len(info['mismatches'])} year(s)")
            print("  Year | CO2 total | AFOLU+Energy sum | Difference")
            print("  ---- | --------- | ----------------- | ----------")
            for col, cval, pval, diff in info['mismatches']:
                print(f"  {col} | {cval:.6g} | {pval:.6g} | {diff:.6g}")

    return results


# TODO cut data that go from 1700

def check_duration_matches_years(protocol_df, df):
    """Verify per-scenario year coverage matches the protocol duration.

    For each scenario present in ``df`` this helper locates the corresponding
    scenario row in ``protocol_df`` and reads a 'duration' column (case-
    insensitive match against any column name containing the substring
    ``'duration'``). It then counts the number of four-digit year columns in
    the grouped subframe that contain at least one non-NaN value and compares
    that count to the declared duration.

    Parameters
    ----------
    protocol_df : pandas.DataFrame
        The protocol DataFrame containing per-scenario metadata including a
        duration column and a ``Scenario`` column.
    df : pandas.DataFrame
        DataFrame to validate; must contain a ``Scenario`` column (case-
        insensitive) and year columns whose names are four-digit strings.

    Returns
    -------
    dict
        Mapping from scenario name -> dict with keys:
        - ``ok`` (bool) whether the counted years equal the protocol duration
        - ``duration`` (int or None) the value read from the protocol table
        - ``counted_years`` (int) number of 4-digit year columns with any
          non-NaN data across the scenario subframe
        - ``years_present`` (list of str) the year column names that had data
        - ``note`` (str) optional explanatory message when protocol lookup
          failed or other issues occurred.

    Notes
    -----
    - The function performs a case-insensitive search for the ``Scenario``
      column in both inputs. It will try exact scenario name matches in the
      protocol table; if not found it will also attempt to match by stripping
      a leading ``'esm-'`` prefix from the scenario name.
    - If no duration-like column is found in ``protocol_df`` the function
      raises ``ValueError``.
    """
    if protocol_df is None or df is None:
        raise ValueError("protocol_df and df must be provided")

    # detect scenario column names (case-insensitive)
    prot_cols_lower = {c.lower(): c for c in protocol_df.columns}
    df_cols_lower = {c.lower(): c for c in df.columns}
    prot_scen_col = prot_cols_lower.get('scenario')
    df_scen_col = df_cols_lower.get('scenario')
    if prot_scen_col is None or df_scen_col is None:
        raise ValueError("Both protocol_df and df must contain a 'Scenario' column")

    # find a duration-like column in protocol_df
    duration_col = None
    for c in protocol_df.columns:
        if 'duration' in str(c).lower():
            duration_col = c
            break
    if duration_col is None:
        raise ValueError('No duration-like column found in protocol_df (searching for "duration")')

    # locate year columns in df (4-digit names)
    year_cols = [c for c in df.columns if isinstance(c, str) and re.fullmatch(r"\d{4}", c)]
    if not year_cols:
        raise ValueError('No 4-digit year columns found in df')

    results = {}
    for scenario, group in df.groupby(df_scen_col):
        # attempt to find the scenario row in protocol_df
        prot_row = protocol_df[protocol_df[prot_scen_col] == scenario]
        note = ''
        if prot_row.empty:
            # try matching without leading 'esm-' prefix
            if scenario.startswith('scen7') and not scenario.endswith('C'):
                short = f"{scenario}C"
            else:
                short = scenario.split('esm-')[-1] if isinstance(scenario, str) else scenario
            prot_row = protocol_df[protocol_df[prot_scen_col] == short]
            if prot_row.empty:
                print(short)
                note = 'protocol lookup failed'
                duration_val = None
            else:
                try:
                    duration_val = int(prot_row.iloc[0][duration_col])
                except Exception:
                    duration_val = None
        else:
            try:
                duration_val = int(prot_row.iloc[0][duration_col])
            except Exception:
                duration_val = None

        # count years with any non-NaN value in the group
        years_present = [y for y in year_cols if group[y].notna().any()]
        counted = len(years_present)

        ok = (duration_val is not None and counted == duration_val)
        results[scenario] = {
            'ok': ok,
            'duration': duration_val,
            'counted_years': counted,
            'years_present': years_present,
            'note': note,
        }

    # print a short summary
    for scen, info in results.items():
        if info['ok']:
            continue
            # print(f"Scenario {scen}: OK (duration {info['duration']} matches {info['counted_years']} years present)")
        else:
            print(f"Scenario {scen}: MISMATCH (declared {info['duration']}, found {info['counted_years']}) {info.get('note','')}")

    return results


def cut_old_rcmip_data_and_add_meta(df, scenario_type, scenario_priority):
    """Trim old RCMIP rows and add metadata columns.

    This helper trims rows that are irrelevant for RCMIP Phase 3 generation and
    inserts two metadata columns, ``Type`` and ``Priority``, in the output
    DataFrame. The trimming rules are:

    - If a ``Region`` column exists, keep only rows with ``Region == 'World'``.
    - Identify year columns with four-digit names (strings like '1750'). If all
      year columns from 1750 onward are NA (or 0) for a row, drop the row
      except when the row's ``Scenario`` is 'esm-piControl' or the
      ``Variable`` is 'Emissions|CO2'.

    Parameters
    ----------
    df : pandas.DataFrame or None
        Input table (emissions/concentrations/forcing). If ``None`` the
        function returns ``None``.
    scenario_type : str
        Type label to insert into the resulting DataFrame's ``Type`` column.
    scenario_priority : int
        Priority value to insert into the resulting DataFrame's ``Priority`` column.

    Returns
    -------
    pandas.DataFrame or None
        A copy of ``df`` with rows trimmed and ``Type``/``Priority`` inserted.
    """
    if df is None:
        return None
    df = df.copy()
    # If a Region column exists, keep only rows for 'World'
    if 'Region' in df.columns:
        df = df[df['Region'] == 'World']

    # Identify year columns (e.g. '1750', '1751', ...) and drop rows that
    # are all-NaN from 1750 onwards, except for special cases.
    year_cols = [c for c in df.columns if re.fullmatch(r"\d{4}", str(c)) and int(c) >= 1750]
    #print(df.columns)
    if year_cols:
        # Rows where all year columns from 1750+ are NaN
        mask_all_nan = (df[year_cols].isna() | (df[year_cols] == 0)).all(axis=1)
        #mask_all_nan = df[year_cols].isna().all(axis=1)
        # Find actual column names for experiment and variable (case-insensitive)
        exp_col = next((c for c in df.columns if c.lower() == 'scenario'), None)
        var_col = next((c for c in df.columns if c.lower() == 'variable'), None)
        #print(exp_col, var_col)
        #sys.exit(4)
        # Only apply the drop if both columns exist; keep rows otherwise
        if exp_col and var_col:
            # Drop rows where years are all NaN AND experiment != 'esm-piControl' AND Variable != 'Emissions|CO2'
            drop_mask = mask_all_nan & (df[exp_col] != 'esm-piControl') & (df[var_col] != 'Emissions|CO2')
            if drop_mask.any():
                df = df.loc[~drop_mask]
    if 'Activity_Id' in df.columns:
        insert_idx = df.columns.get_loc('Activity_Id') + 1
        if 'Type' in df.columns:
            df['Type'] = scenario_type
        else:
            df.insert(insert_idx, 'Type', scenario_type)
        if 'Priority' in df.columns:
            df['Priority'] = scenario_priority
        else:
            df.insert(insert_idx + 1, 'Priority', scenario_priority)
    else:
        df['Type'] = scenario_type
        df['Priority'] = scenario_priority
    return df

def fill_in_nans(df, scenario_df):
    """Fill leading and trailing NaNs per-row across year columns.

    For each row the function locates year columns (4-digit names >= 1750) and
    applies three rules:

    1. Leading NaNs are replaced by the first non-NaN value in the row.
    2. Trailing NaNs are replaced by the last non-NaN value in the row.
    3. If NaNs occur in the middle (i.e., between two non-NaN values), the
       function will attempt to fill them by linear interpolation or constant
       filling when possible; if interpolation is ambiguous the function prints
       the problematic row and exits with ``sys.exit(4)`` (preserves original
       requested behaviour).

    Parameters
    ----------
    df : pandas.DataFrame or None
        DataFrame containing year columns to operate on. If ``None``, ``None``
        is returned.

    Returns
    -------
    pandas.DataFrame or None
        A copy of ``df`` with leading/trailing NaNs filled. May call
        ``sys.exit(4)`` if a non-trivial middle-NaN case is detected that
        requires manual inspection.
    """
    if df is None:
        return df
    scenario_df.set_index('Scenario', inplace=True)
    df = df.copy()
    # Identify year columns (4-digit strings >= 1750)
    year_cols = [c for c in df.columns if re.fullmatch(r"\d{4}", str(c)) and int(str(c)) >= 1750]
    if not year_cols:
        return df

    # Work on string column names for consistent indexing
    year_cols_str = [str(c) for c in year_cols]

    for idx in df.index:
        row = df.loc[idx, year_cols_str]
        # treat zeros as valid values (do not consider them NaN)
        is_nan = row.isna()
        if not is_nan.any():
            continue

        # Find first and last non-nan positions
        non_nan_idx = np.where(~is_nan.values)[0]
        if non_nan_idx.size == 0:
            # Entire row is NaN: leave as-is (or could fill with zeros)
            continue

        first_non_nan = non_nan_idx[0]
        last_non_nan = non_nan_idx[-1]

        # Leading NaNs: positions before first_non_nan
        if first_non_nan > 0 and is_nan.values[:first_non_nan].any():
            
            first_val = row.iloc[first_non_nan]
            if first_val != 0.0:
                if "esm-bell" in df.loc[idx, "Scenario"] or "esm-flat" in df.loc[idx, "Scenario"]:
                    first_val = 0.0
                elif "Concentrations" in df.loc[idx, "Variable"]:
                    first_val = first_val
                else:
                    print(first_val)
                    print(df.loc[idx, year_cols_str[first_non_nan-2: first_non_nan+2]])
                    print(df.loc[idx, :])
                    sys.exit(4)

            # broadcast first_val to leading positions
            df.loc[idx, year_cols_str[:first_non_nan]] = first_val

        # Trailing NaNs: positions after last_non_nan
        if last_non_nan < (len(year_cols_str) - 1) and is_nan.values[last_non_nan + 1 :].any():
            last_val = row.iloc[last_non_nan]
            scenario = df.loc[idx, "Scenario"]
            # For now to avoid weird behaviour for CMIP7 concentrations, lets leave trailing NaNs
            if "scen7" in scenario:
                break
            if last_non_nan >= scenario_df.loc[scenario, "Duration of scenario"]-1:
                pass
            elif "hist" in df.loc[idx, "Scenario"] and last_val == 0.0:
                df.loc[idx, year_cols_str[last_non_nan+1:scenario_df.loc[scenario, "Duration of scenario"]]] = last_val
            elif last_val != 0.0:
                if "hist" in df.loc[idx, "Scenario"]:
                    pass
                sys.exit(4)
            else:    
                df.loc[idx, year_cols_str[last_non_nan+1:]] = last_val

        # Middle NaNs: any NaNs between first_non_nan and last_non_nan
        if (is_nan.values[first_non_nan:last_non_nan + 1]).any():
            # Print the problematic row and exit as requested
            # Print row slice for debugging
            unique_non_nan_values = np.unique(row.iloc[first_non_nan:last_non_nan + 1].dropna().values)
            if len(unique_non_nan_values) == 1:
                # All non-nan values are the same, fill middle nans with that value
                fill_value = unique_non_nan_values[0]
                df.loc[idx, year_cols_str[first_non_nan+1:last_non_nan]] = fill_value
            else:
                
                prev_non_nan = non_nan_idx[0]
                for non_nan in non_nan_idx[1:]:
                    if non_nan == prev_non_nan + 1:
                        prev_non_nan = non_nan
                        continue
                    # Really this bit could be handled in one go
                    if row.iloc[prev_non_nan] == row.iloc[non_nan]:
                        # fill in between with this value
                        fill_value = row.iloc[prev_non_nan]
                        df.loc[idx, year_cols_str[prev_non_nan+1:non_nan]] = float(fill_value)
                        prev_non_nan = non_nan
                    else:
                        xp = np.array([prev_non_nan, non_nan])
                        fp = np.array([row.iloc[prev_non_nan], row.iloc[non_nan]])
                        x = np.arange(prev_non_nan + 1, non_nan)
                        df.loc[idx, year_cols_str[prev_non_nan+1:non_nan]] = np.interp(x, xp, fp)
                        prev_non_nan = non_nan

    return df

def fix_fgases_methanemip(methanemip_df: pd.DataFrame, correct_variable_names: list) -> pd.DataFrame:
    """Apply special handling to F-gases in Methanemip data.

    This function identifies rows in the provided ``methanemip_df`` that
    correspond to F-gas emissions and applies specific renaming rules to align
    them with RCMIP variable names. The function also checks for the presence of
    expected variables and prints warnings if they are missing.

    Parameters
    ----------
    methanemip_df : pandas.DataFrame
        DataFrame containing Methanemip emissions data, expected to have a
        ``Variable`` column and year columns.

    Returns
    -------
    pandas.DataFrame
        The updated DataFrame with F-gas variables renamed according to RCMIP
        conventions where applicable.
    """
    for idx, row in methanemip_df.iterrows():
        variable = row['Variable']
        if variable in correct_variable_names:
            continue
        variable_main = variable.split("|")[-1]
        candidate = "Emissions|F-Gases|" + variable_main.replace("-", "")
        if candidate in correct_variable_names:
            methanemip_df.at[idx, 'Variable'] = candidate
            continue
        candidate = "Emissions|F-Gases|HFC|" + variable_main.replace("-", "")
        if candidate in correct_variable_names:
            methanemip_df.at[idx, 'Variable'] = candidate
            continue
        candidate = "Emissions|F-Gases|PFC|" + variable_main.replace("-", "")
        if candidate in correct_variable_names:
            methanemip_df.at[idx, 'Variable'] = candidate
            continue
        candidate = "Emissions|Montreal Gases|" + variable_main.replace("-", "")
        if candidate in correct_variable_names:
            methanemip_df.at[idx, 'Variable'] = candidate
            continue
        candidate = "Emissions|Montreal Gases|CFC|" + variable_main.replace("-", "")
        if candidate in correct_variable_names:
            methanemip_df.at[idx, 'Variable'] = candidate
            continue
        sys.exit(4)
    return methanemip_df

def fix_units_methanemip(methanemip_subframe: pd.DataFrame, unique_couples:pd.DataFrame):
    """
    Fix methanemip units to match the units in the unique_couples DataFrame.

    This function checks the units of the Methanemip data against 
    the expected units in the unique_couples DataFrame. If a mismatch 
    is found for a variable, it adjusts the units in the Methanemip 
    subframe to match the expected units and applies a conversion 
    factor to the values if necessary (e.g., converting from M to k 
    or vice versa).

    Parameters
    ----------
    methanemip_subframe : pandas.DataFrame
        DataFrame containing the Methanemip data for a specific scenario, 
        expected to have a 'Variable' column and a 'Unit' column.
    unique_couples : pandas.DataFrame

    Returns
    -------
    pandas.DataFrame
        The updated Methanemip subframe with units fixed to match the 
        unique_couples DataFrame where applicable.
    """
    for idx, row in unique_couples.iterrows():
        variable = row['Variable']
        unit = row['Unit']
        mask_methane = methanemip_subframe['Variable'] == variable
        if len(mask_methane) > 0:
            if methanemip_subframe.loc[mask_methane, 'Unit'].values[0] != unit:
                old_unit = methanemip_subframe.loc[mask_methane, 'Unit'].unique()[0]
                print("Adjusting methanemip unit for variable", variable, "from", old_unit, "to", unit)
                methanemip_subframe.loc[mask_methane, 'Unit'] = unit
                if old_unit.replace("M", "k") == unit:
                    print("Also adjusting values for M/k difference")
                    for col in methanemip_subframe.columns:
                        if re.fullmatch(r"\d{4}", str(col)):
                            methanemip_subframe.loc[mask_methane, col] = methanemip_subframe.loc[mask_methane, col] * 1e3
                elif old_unit.replace("k", "M") == unit:
                    print("Also adjusting values for k/M difference")
                    for col in methanemip_subframe.columns:
                        if re.fullmatch(r"\d{4}", str(col)):
                            methanemip_subframe.loc[mask_methane, col] = methanemip_subframe.loc[mask_methane, col] / 1e3
    return methanemip_subframe

    

def map_in_cmip7_style_emissions_to_rcmip_format(input_em_df, cmip7_em_df, year_str_col):
    """Map CMIP7-style emissions into the RCMIP emissions template.

    This function iterates rows of an RCMIP-style ``emissions_df`` and
    attempts to populate the provided ``year_str_col`` columns from a
    CMIP7-format source frame ``cmip7_em_df`` (which contains a ``variable``
    column). Special-case handling includes summing AFOLU + Energy for CO2 and
    renaming of some F-gas variables.

    Parameters
    ----------
    emissions_df : pandas.DataFrame
        Target RCMIP-style DataFrame to fill. Must contain a ``Variable``
        column describing the emission variable names.
    cmip7_em_df : pandas.DataFrame
        Source DataFrame containing a ``variable`` column and columns for years
        (matching the entries in ``year_str_col``).
    year_str_col : list of str
        List of year column names (strings) to copy from ``cmip7_em_df`` into
        ``emissions_df``.

    Returns
    -------
    pandas.DataFrame
        The updated ``emissions_df`` with the requested year columns populated
        where matches were found.
    """
    substrings_to_check = ["F-Gases", "Montreal Gases", "CFC|"] 
    rename_dict_em = {
        "Emissions|F-Gases|HFC|HFC4310mee": "Emissions|HFC|HFC43-10"
    }
    if "Variable" in cmip7_em_df.columns:
        cmip7_em_df = cmip7_em_df.rename(columns={"Variable": "variable"})
    emissions_df = input_em_df.copy()
    for idx, row in input_em_df.iterrows():
        variable = row['Variable']
        data_source = cmip7_em_df.loc[cmip7_em_df['variable'] == variable]
        if not data_source.empty:
            emissions_df.loc[idx, year_str_col] = data_source[ year_str_col].values[0]
             # Initialize with NaN  # Initialize with NaN
            #sys.exit(4)
            # for year in range(1750, 2024):
            #     year_str = str(year)
            #     if year_str in data_source.columns:
            #         emissions_df.at[idx, year_str] = data_source[year_str].values[0]
        elif variable == "Emissions|CO2":
            # For Emissions|CO2, sum AFOLU and Energy and Industrial Processes
            data_source_afolu = cmip7_em_df.loc[cmip7_em_df['variable'] == "Emissions|CO2|AFOLU"]
            data_source_energy = cmip7_em_df.loc[cmip7_em_df['variable'] == "Emissions|CO2|Energy and Industrial Processes"]
            if data_source_afolu.empty or data_source_energy.empty:
                print(f"Warning: missing AFOLU or Energy data for variable {variable} in emissions source data.")
                sys.exit(4)
            emissions_df.loc[idx, year_str_col] = (
                data_source_afolu[year_str_col].values[0] +
                data_source_energy[year_str_col].values[0]
            )
        elif variable in rename_dict_em or any(sub in variable for sub in substrings_to_check):
            # Handle renamed variables
            if variable in rename_dict_em:
                new_variable = rename_dict_em[variable]
            elif "F-Gases|PFC" in variable:
                new_variable = variable.replace("|F-Gases|PFC|", "|")
            elif "Montreal Gases|CFC" in variable:
                new_variable = variable.replace("|Montreal Gases|CFC|", "|")
            elif "F-Gases" in variable:
                new_variable = variable.replace("|F-Gases|", "|")
            elif "Montreal Gases" in variable:
                new_variable = variable.replace("|Montreal Gases|", "|")
            else: 
                print("Why am I here?")
                print(variable)
                sys.exit(4)
            data_source_renamed = cmip7_em_df.loc[cmip7_em_df['variable'] == new_variable]
            if data_source_renamed.empty:
                print(cmip7_em_df["variable"].unique())
                print(f"Warning: missing data for renamed variable {new_variable} in emissions source data.")
                sys.exit(4)
            emissions_df.loc[idx, year_str_col] = data_source_renamed[year_str_col].values[0]
        else:
            print(cmip7_em_df["variable"].unique())
            print(f"Warning: variable {variable} not found in emissions source data.")
            sys.exit(4)
    emissions_df["Model"] = "emissions_harmonisation_pipeline"
    emissions_df["Mip_Era"] = "CMIP7"
    return emissions_df

def fix_columns_to_int_string(df):
    """Convert column names ending with '.0' to integer string format.
    
    This function renames DataFrame columns that end with '.0' (e.g., '1750.0')
    to their integer string equivalents (e.g., '1750').
    
    Parameters
    ----------
    df : pandas.DataFrame
        DataFrame with potentially float-formatted year column names.
    
    Returns
    -------
    pandas.DataFrame
        The same DataFrame with renamed columns.
    """
    rename_map = {}
    for col in df.columns:
        col_str = str(col)
        if col_str.endswith('.0'):
            rename_map[col] = col_str[:-2]
    
    if rename_map:
        df = df.rename(columns=rename_map)
    
    return df

def merge_cmip7_concentrations_data(conc_df, src_conc):
    """Merge CMIP7 concentration data from source into the provided DataFrame.

    This function takes an RCMIP-style concentrations DataFrame and merges in
    data from a CMIP7-format source DataFrame (which contains a 'variable'
    column and year columns). It handles variable name mapping/renaming as
    needed to align the two datasets.

    Parameters
    ----------
    conc_df : pandas.DataFrame
        Target RCMIP-style concentrations DataFrame to fill. Must contain a
        ``Variable`` column describing the concentration variable names.
    src_conc : pandas.DataFrame
        Source DataFrame containing a ``variable`` column and columns for years
        (matching the year columns in ``conc_df``).
    Returns
    ------- 
    pandas.DataFrame
        The updated ``conc_df`` with the requested year columns populated where
        matches were found.
    """
    rename_dict_conc = {
        "HFC4310mee": "HFC-43-10mee",
        "C4F10": "n-C4F10",
        "C5F12": "n-C5F12",
        "cC4F8": "c-C4F8"
    }
    for idx, row in conc_df.iterrows():
        variable = row['Variable'].split("|")[-1]  # Get the part after the '|'
        if variable.upper() in src_conc.columns:
            variable = variable.upper()
        elif variable not in src_conc.columns:
            if variable in rename_dict_conc:
                variable = rename_dict_conc[variable]
            elif variable.startswith("HCFC"):
                variable = variable.replace("HCFC", "HCFC-")
            elif variable.startswith("CFC"):
                variable = variable.replace("CFC", "CFC-")
            elif variable.startswith("HFC") :
                variable = variable.replace("HFC", "HFC-")
            elif variable.startswith("Halon"):
                variable = variable.replace("Halon", "Halon-")
        if variable in src_conc.columns:
            conc_df.loc[idx, src_conc["Year"].values.astype(str)] = src_conc[variable].values
        else:
            print(f"Warning: variable {variable} not found in concentrations source data.")
    return conc_df

def get_cmip7_historical_data(emissions_df, concentrations_df):
    """Populate CMIP7 historical emissions and concentrations from source CSVs.

    The function reads two external CSVs (historical emissions harmonisation and
    IGCC concentrations), maps/renames variables where required and fills the
    provided ``emissions_df`` and ``concentrations_df`` year columns for the
    CMIP7 'historical' period (1750 onwards).

    Parameters
    ----------
    emissions_df : pandas.DataFrame
        RCMIP-style emissions DataFrame (filtered to the 'historical' scenario).
    concentrations_df : pandas.DataFrame
        RCMIP-style concentrations DataFrame (filtered to the 'historical' scenario).

    Returns
    -------
    tuple
        Tuple ``(emissions_df, concentrations_df)`` with year columns filled
        from the external source files where available.

    Notes
    -----
    - If source CSVs are missing or unreadable the function returns the input
      frames unchanged and prints a warning (no exception is raised).
    - The function may call ``sys.exit(4)`` for unrecoverable mapping errors
      (preserves original behaviour).
    """
    # Paths to source data
    em_path = os.path.join('..', 'data', 'historical_emissions_from_emissions_harmonisation_Feb2026.csv')
    conc_path = os.path.join('..', 'data', 'CMIP-1750-2022_output_concentrations.csv')

    try:
        src_em = pd.read_csv(em_path)
    except Exception as e:
        # print(f"Warning: could not read emissions source file {em_path}: {e}")
        return emissions_df, concentrations_df
    # print(src_em.head())
    # print(src_em["variable"].unique())
    #sys.exit(4)
    try:
        src_conc = pd.read_csv(conc_path)
    except Exception as e:
        # print(f"Warning: could not read concentrations source file {conc_path}: {e}")
        return emissions_df, concentrations_df
    
    # Adding together concentrations for n-C6F14 and i-C6F14 into total C6F14
    # This is for concentrations from igcc, not needed for esgf-version
    if "n-C6F14" in src_conc.columns and "i-C6F14" in src_conc.columns:
        src_conc["C6F14"] = src_conc["n-C6F14"] + src_conc["i-C6F14"]
    year_str_col = [str(y) for y in range(1750, 2024)]
    emissions_df = map_in_cmip7_style_emissions_to_rcmip_format(emissions_df, src_em, year_str_col)
    # TODO: Add in gases that are not in the old rcmip dataset
    concentrations_df = merge_cmip7_concentrations_data(concentrations_df, src_conc) 

    return emissions_df, concentrations_df

def check_experiment_contents_of_generated_data(protocol_df, version="v2.0.0"):
    """Compare protocol scenarios to already-generated RCMIP outputs.

    Reads the three generated CSVs in ``../data/generated_inputs/`` and checks
    which protocol scenarios are already present in the generated tables. The
    function gathers lists of existing and missing scenarios for emissions and
    concentrations and is primarily used for QA/debugging.

    Parameters
    ----------
    protocol_df : pandas.DataFrame
        Protocol sheet dataframe containing at least a ``Scenario`` column.

    Returns
    -------
    None
        The function prints lists of existing/missing scenarios (side-effect).
    """
    emissions_generated = pd.read_csv(f"{rcmip_output_dir}rcmip_phase3_emissions_{version}.csv")
    concentrations_generated = pd.read_csv(f"{rcmip_output_dir}rcmip_phase3_concentrations_{version}.csv")
    # print(protocol_df.head())
    em_existing = []
    em_missing = []
    conc_existing = []
    conc_missing = []
    for scenario in protocol_df['Scenario'].unique():
        # print(scenario)
        if scenario in emissions_generated['Scenario'].values:
            em_existing.append(scenario)
        elif scenario.split("esm-")[-1] in emissions_generated['Scenario'].values:
            em_existing.append(scenario)
        elif scenario.split("esm-allGHG")[-1] in emissions_generated['Scenario'].values:
            em_existing.append(scenario)
        elif "scen7" not in scenario and not any([scenario.startswith(nostart) for nostart in ["1pct", "abrupt", "methanemip"]]) and "brch" not in scenario:
            em_missing.append(scenario) 
        if scenario in concentrations_generated['Scenario'].values:
            conc_existing.append(scenario)
        elif "esm-" not in scenario and "scen7" not in scenario and not scenario.startswith("methanemip"):
            conc_missing.append(scenario)

    print("Checking CO2 component sums...")
    check_co2_component_sums(emissions_generated)
    print("Checking durations for emissions...")
    check_duration_matches_years(protocol_df, emissions_generated)
    print("Checking durations for concentrations...")
    check_duration_matches_years(protocol_df, concentrations_generated)


def simplify_fgas_and_montreal_gases(df):
    """Simplify variable names for F-gases and Montreal protocol gases.

    Many source files use nested groupings such as ``F-Gases|HFC|...`` or
    ``Montreal Gases|CFC|...``. This function flattens those names by removing
    the intermediate grouping so downstream code can match a consistent
    variable naming convention.

    Parameters
    ----------
    df : pandas.DataFrame
        DataFrame with a ``Variable`` column to be simplified in-place.

    Returns
    -------
    pandas.DataFrame
        The same DataFrame instance with the ``Variable`` column modified.
    """
    substrings_to_check = ["F-Gases|", "Montreal Gases|", "CFC|"] 
    for idx, row in df.iterrows():
        variable = row['Variable']
        if any(sub in variable for sub in substrings_to_check):
            if "F-Gases|PFC" in variable:
                new_variable = variable.replace("|F-Gases|PFC|", "|")
            if "F-Gases|HFC" in variable:
                new_variable = variable.replace("|F-Gases|HFC|", "|")
            elif "Montreal Gases|CFC" in variable:
                new_variable = variable.replace("|Montreal Gases|CFC|", "|")
            elif "F-Gases" in variable:
                new_variable = variable.replace("|F-Gases|", "|")
            elif "Montreal Gases" in variable:
                new_variable = variable.replace("|Montreal Gases|", "|")
            else:
                print(df.head())
                print(variable) 
                print("Why am I here?")
                sys.exit(4)
            df.at[idx, 'Variable'] = new_variable
    return df

def _prefix_emissions(val):
    """Normalise a short variable label into an RCMIP-style variable name.

    The helper converts a compact variable string into the long pipe-separated
    RCMIP style used throughout the repository. Specific special cases are
    handled for common short names used in input files.

    Parameters
    ----------
    val : object
        A value from a dataframe Variable column. May be a string or NaN.

    Returns
    -------
    str or object
        If ``val`` is NaN the original value is returned. Otherwise a string
        is returned in the form ``Emissions|<...>`` or an already-normalised
        string is returned unchanged.

    Examples
    --------
    >>> _prefix_emissions('CO2 AFOLU')
    'Emissions|CO2|AFOLU'
    >>> _prefix_emissions('Emissions|CH4')
    'Emissions|CH4'
    """
    if pd.isna(val):
        return val
    s = str(val).strip()
    # If already contains a pipe-grouped name, leave unchanged; otherwise prefix
    if '|' in s or s.startswith('Emissions|'):
        return s
    if s == "CO2 AFOLU":
        return "Emissions|CO2|AFOLU"
    if s == "CO2 FFI":
        return "Emissions|CO2|Energy and Industrial Processes"
    return f"Emissions|{s}"

def read_and_restructure_methanemip_data(methanemip_file_path= "../data/methane-emissions-scenarios.csv"):
    """Read and restructure MethaneMIP data from the provided CSV file.

    This function reads a MethaneMIP CSV file and restructures it into an
    RCMIP-style DataFrame with appropriate columns and metadata.

    Parameters
    ----------
    methanemip_file_path : str
        Path to the MethaneMIP CSV file.

    Returns
    -------
    pandas.DataFrame
        Restructured DataFrame suitable for RCMIP input generation.
    """
    methanemip_df = pd.read_csv(methanemip_file_path)
    #print(methanemip_df.head(10))
    #print(methanemip_df.shape)
    # Rename columns like '1750.5' -> '1750'
    rename_map = {}
    for col in methanemip_df.columns:
        col_str = str(col)
        m = re.match(r'^(\d{4})\.5$', col_str)
        if m:
            rename_map[col] = m.group(1)
    if rename_map:
        methanemip_df = methanemip_df.rename(columns=rename_map)
    # Ensure we operate on the variable column (case-insensitive)
    var_col = next((c for c in methanemip_df.columns if c.lower() == 'variable'), None)
    if var_col is None:
        print("Warning: no 'Variable' column found in MethaneMIP data; skipping prefixing.")
    else:
        methanemip_df[var_col] = methanemip_df[var_col].apply(_prefix_emissions)
        # Add "Model" as the first column with the required constant value
        if 'Model' in methanemip_df.columns:
            methanemip_df['Model'] = 'MESSAGE-GLOBIOM-methanemip'
        else:
            methanemip_df.insert(0, 'Model', 'MESSAGE-GLOBIOM-methanemip')

        # Insert "Activity_Id" immediately after "Unit" (or append if "Unit" not present)
        if 'Activity_Id' in methanemip_df.columns:
            methanemip_df['Activity_Id'] = 'not_applicable'
        else:
            if 'Unit' in methanemip_df.columns:
                insert_idx = methanemip_df.columns.get_loc('Unit') + 1
                methanemip_df.insert(insert_idx, 'Activity_Id', 'not_applicable')
            else:
                methanemip_df['Activity_Id'] = 'not_applicable'
    # add CO2 totals combining AFOLU + Energy rows
    var_col_name = next((c for c in methanemip_df.columns if c.lower() == 'variable'), None)
    scenarios = methanemip_df['Scenario'].unique()

    year_cols = [c for c in methanemip_df.columns if re.fullmatch(r"\d{4}", str(c))]
    key_cols = [c for c in ['Unit', 'Region', 'Scenario', 'Model', 'Activity_Id'] if c in methanemip_df.columns]
    new_rows = []
    df_afolu = methanemip_df[methanemip_df[var_col_name] == "Emissions|CO2|AFOLU"].copy()
    df_energy = methanemip_df[methanemip_df[var_col_name] == "Emissions|CO2|Energy and Industrial Processes"].copy()
    for scenario in scenarios:
        df_afolu = methanemip_df[(methanemip_df[var_col_name] == "Emissions|CO2|AFOLU") & (methanemip_df['Scenario'] == scenario)].copy().iloc[0,:]
        df_energy = methanemip_df[(methanemip_df[var_col_name] == "Emissions|CO2|Energy and Industrial Processes") & (methanemip_df['Scenario'] == scenario)].copy().iloc[0,:]
        #print(df_afolu.shape, df_energy.shape)
        # start from AFOLU metadata and replace Variable + year values with sums
        new_row = df_afolu.copy()
        new_row[var_col_name] = "Emissions|CO2"
        new_row.loc[year_cols] = df_afolu.loc[year_cols].values + df_energy.loc[year_cols].values

        for kc in key_cols:
            new_row[kc] = df_afolu.get(kc, np.nan)
        new_rows.append(new_row)

    if new_rows:
        new_df = pd.DataFrame(new_rows)
        # preserve original column order where possible
        methanemip_df = pd.concat([methanemip_df, new_df], ignore_index=True, sort=False)
    else:
        print("Warning: no 'Variable' column found; cannot create aggregated Emissions|CO2 rows.")
    # Restructure as needed (this is a placeholder; actual restructuring logic may vary)
    # For example, renaming columns, adding metadata, etc.
    methanemip_df = cut_old_rcmip_data_and_add_meta(methanemip_df, "non-idealised", 1)
    return methanemip_df

def read_cmip7_lucalbedo_and_irrigation_forcing(scenario):
    """Read CMIP7 land-use albedo and irrigation forcing data from source CSVs.

    This function reads two external CSVs containing CMIP7 land-use albedo
    and irrigation effective radiative forcing data, applies necessary
    metadata, and returns them as DataFrames.

    Returns two DataFrames: albedo_forcing_df and irrigation_forcing_df.
    """
    albedo_path = os.path.join('..', 'data', 'Forcing_AFOLU_CO2.csv')
    irrigation_path = os.path.join('..', 'data', 'Forcing_irrigation_population_scale.csv')
    template_data_dict = {
        'Scenario': [scenario],
        'Unit': ['W/m^2'],
        'Region': ['World'],
        'Model': ['Unspecified'],
        'Activity_Id': ['not_applicable'],
        'Type': ['non-idealised'],
        'Priority': [1],
        'Mip_Era': ['CMIP7'],
    }

    df_full_luc = pd.read_csv(albedo_path)
    df_full_irrigation = pd.read_csv(irrigation_path)

    # Taking historical data from one of the scenarios (1750-2023 is equal across all)
    # VL is a random choice here, could be any of them
    if scenario == "historical":
        year_stop = 2023
        scen_idx = "VL"
    else:
        year_stop = 2501
        scen_idx = scenario.split("scen7-")[-1]
    data_dict_luc = template_data_dict.copy()
    data_dict_luc["Variable"] = "Effective Radiative Forcing|Anthropogenic|Albedo Change|Land Use"
    data_dict_irrigation = template_data_dict.copy()
    data_dict_irrigation["Variable"] = "Effective Radiative Forcing|Anthropogenic|Albedo Change|Irrigation"
    data_dict_tot = template_data_dict.copy()
    data_dict_tot["Variable"] = "Effective Radiative Forcing|Anthropogenic|Albedo Change"
    for year in range(1750, year_stop):
        year_str = str(year)
        data_dict_irrigation[year_str] = df_full_irrigation[scen_idx].values[year-1750]
        data_dict_luc[year_str] = df_full_luc[scen_idx].values[year-1750]
        data_dict_tot[year_str] = data_dict_luc[year_str] + data_dict_irrigation[year_str]
    albedo_df = pd.DataFrame(data_dict_luc)
    irrigation_df = pd.DataFrame(data_dict_irrigation)
    tot_df = pd.DataFrame(data_dict_tot)
    return albedo_df, irrigation_df, tot_df

def add_cmip7_land_use_forcing_to_df(df, scenario):
    """Add CMIP7 land-use albedo and irrigation forcing data to the provided DataFrame.

    This method reads CMIP7 land-use albedo and irrigation forcing data using
    ``read_cmip7_lucalbedo_and_irrigation_forcing`` and appends it to the
    provided DataFrame ``df``.

    Parameters
    ----------
    df : pandas.DataFrame
        The DataFrame to which the land-use forcing data will be appended.
    scenario : str
        The scenario name to use when reading the CMIP7 forcing data.

    Returns
    -------
    pandas.DataFrame
        The updated DataFrame with CMIP7 land-use forcing data appended.
    """
    albedo_df, irrigation_df, luctot_df = read_cmip7_lucalbedo_and_irrigation_forcing(scenario)
    df = pd.concat([df, albedo_df, irrigation_df, luctot_df], ignore_index=True, sort=False)
    return df

def get_scen7_subframe_just_one(
        scen7_inputs: pd.DataFrame, 
        scen_name_full: str, 
        scen_model: str, 
        scenario_name: str, 
        scenario_type: str, 
        scenario_priority: str, 
        scen7_history: pd.DataFrame
        ):
    """
    Take scen7_inputs and cut to a single scenario modle combo for a year, trim and fix
    """
    scen_subframe = scen7_inputs.loc[(scen7_inputs['scenario']==scen_name_full) & (scen7_inputs['region']=='World')]# & (str(scen7_inputs['model']).startswith(scen_model))]
    scen_subframe = scen_subframe.loc[scen7_inputs['model'].str.startswith(scen_model)]
    if scen_subframe.empty:
        print(f"Warning: no data found for scenario7 scenario {scenario_name} with model {scen_model}.")
        sys.exit(4)
    year_cols = [col for col in scen7_inputs.columns if isinstance(col, str) and col.isdigit() and len(col) == 4]
    scen7_history = map_in_cmip7_style_emissions_to_rcmip_format(scen7_history, scen_subframe, year_cols)
    # print("Year columns in scen7_inputs:", year_cols)
    model_full = scen_subframe['model'].values[0]
    em_df = cut_old_rcmip_data_and_add_meta(scen7_history, scenario_type, scenario_priority)
    return em_df, model_full

def get_scenarioinfor_from_dict(scen_name_short, scenario_name):
    scenario_name_dictionary = {
            "M": ['SSP2 - Medium Emissions', "IMAGE"], 
            "ML": ['SSP2 - Medium-Low Emissions', "COFFEE"],
            "HL": ['SSP5 - Medium-Low Emissions_a',"WITCH"], 
            "VL": ['SSP1 - Very Low Emissions', "REMIND"],
            "L": ['SSP2 - Low Emissions', "MESSAGE"],
            "H": ['SSP3 - High Emissions', "GCAM"],
            "LN": ['SSP2 - Low Overshoot_a', "AIM"], 
        }
    if scen_name_short not in scenario_name_dictionary:
        print(f"Warning: scenario7 scenario name {scenario_name} not recognised.")
        sys.exit(4)
    scen_name_full = scenario_name_dictionary[scen_name_short][0]

    scen_model = scenario_name_dictionary[scen_name_short][1]
    return scen_name_full, scen_model


class RcmipInputGenerator:
    """Generator for RCMIP Phase 3 per-scenario input tables.

    This class orchestrates reading the protocol Excel sheet, loading existing
    RCMIP source CSVs, applying historical fills (CMIP7), and producing the
    per-scenario emissions/concentrations/forcing CSV files under
    ``../data/generated_inputs/``.
    """
    def __init__(
            self, 
            excel_file_path, 
            cmip7_scenario_file_path = "../ScenarioMIP/scenariomip_scenarios_for_scms.csv",
            version = "v2.0.0"
            ):
        """Initialise the RCMIP input generator and load source data.

        This constructor reads the RCMIP protocol ``scenario_info`` sheet and
        loads a set of previously-generated RCMIP CSVs required as source
        material. It also prepares historical CMIP7 data and restructures the
        MethaneMIP input table used by some scenario handlers.

        Parameters
        ----------
        excel_file_path : str
            Path to the RCMIP protocol Excel file containing the
            ``scenario_info`` sheet.
        cmip7_scenario_file_path : str, optional
            Path to a CSV file containing CMIP7-style scenario inputs. The
            default points to ``../ScenarioMIP/scenariomip_scenarios_for_scms.csv``.

        Attributes
        ----------
        scenario_df : pandas.DataFrame
            The protocol ``scenario_info`` sheet loaded into a DataFrame.
        emissions_df, concentrations_df, forcing_df : pandas.DataFrame
            Source RCMIP CSV tables read from ``../data``.
        emissions_df_hist, concentrations_df_hist : pandas.DataFrame
            Historical data filled from CMIP7 sources (created by
            ``get_historical_data_for_cmip7``).
        cmip7_scenario_file_path : str
            Path to the CMIP7 scenario CSV used by ``make_scenario7_data``.
        methanemip_df : pandas.DataFrame
            Restructured MethaneMIP table prepared for methanemip scenarios.
        """
        self.version = version
        self.scenario_df = pd.read_excel(excel_file_path, sheet_name="scenario_info", 
                                skiprows=[1, 2], usecols=[1, 2, 5, 6])
        self.read_in_previous_rcmip_data()
        self.get_historical_data_for_cmip7()
        self.cmip7_scenario_file_path = cmip7_scenario_file_path
        self.methanemip_df = read_and_restructure_methanemip_data()
        self.fix_piControl()
        self.fix_piControl(type="concentrations")
        self.fix_old_conc_idealised()
        self.fix_old_em_idealised()


    def fix_old_em_idealised(self):
        """
        Overwrite idealised emissions with 1850 values from historical data.
        
        This method ensures that idealised emissions scenarios (e.g.,
        ``1pctCO2``) have constant values across all year columns by copying the
        1850 value from the historical emissions data for each variable. (i.e.
        same as used in piControl). Then for CO2 that the idealised emissions
        trajectory as defined by the experiment is added to this, such as the 
        bell curve or pulse emissions (esm-bell-1000PgC... esm-pi-CO2pulse, esm-pi-cdr-pulse).
        """
        print("Fixing idealised emissions...")
        experiments = self.emissions_df["Scenario"].unique().tolist()
        to_add = []
        for experiment in experiments:
            if (not experiment.startswith("esm")) or (experiment == "esm-piControl"):
                continue
            piControl_copy = self.emissions_df[self.emissions_df["Scenario"] == "esm-piControl"].copy()
            data_input = self.emissions_df[self.emissions_df["Scenario"] == experiment].copy()

            year_cols_post_1750 = [c for c in piControl_copy.columns 
                            if isinstance(c, str) and re.fullmatch(r"\d{4}", c) and int(c) >= 1750]
            for idx, row in data_input.iterrows():
                variable = row['Variable']
                variable_data = data_input.loc[idx, year_cols_post_1750].values
                if variable_data.sum() == 0 or np.isnan(variable_data).all():
                    continue
                piControl_data = piControl_copy.loc[piControl_copy["Variable"] == variable, year_cols_post_1750].values[0]

                piControl_copy.loc[piControl_copy["Variable"] == variable, year_cols_post_1750] = np.nansum((piControl_data, variable_data), axis = 0)
                piControl_copy["Scenario"] = experiment
            to_add.append(piControl_copy)

            self.emissions_df = self.emissions_df[self.emissions_df["Scenario"] != experiment].copy()
        self.emissions_df = pd.concat([self.emissions_df] + to_add, ignore_index=True)

    def fix_old_conc_idealised(self):
        """
        Overwrite idealised concentrations with 1850 values from historical data.
        
        This method ensures that idealised concentration scenarios (e.g.,
        ``1pctCO2``) have constant values across all year columns by copying the
        1850 value from the historical concentrations data for each variable.
        (i.e. same as used in piControl).
        Then that value is used to define the idealised trajectory from 1850
        onwards, increasing by 1% per year for 1pctCO2 / 1pctCO2-4xt, or doubled,
        quadrupled, or halved for abrupt-2xCO2, abrupt-4xCO2, or abrupt-0.5xCO2 respectively.
        """
        print("Fixing idealised concentrations...")
        experiments = self.concentrations_df["Scenario"].unique().tolist()
        co2_var_name = 'Atmospheric Concentrations|CO2'
        co2_piC = self.concentrations_df.loc[(self.concentrations_df["Scenario"] == "piControl") & (self.concentrations_df["Variable"] == co2_var_name)]["1850"].values[0]

        to_add = []
        for experiment in experiments:
            if not (experiment.startswith("1pct") or experiment.startswith("abrupt")):
                continue
            piControl_copy = self.concentrations_df[self.concentrations_df["Scenario"] == "piControl"].copy()

            year_cols_post_1850 = [c for c in piControl_copy.columns 
                            if isinstance(c, str) and re.fullmatch(r"\d{4}", c) and int(c) >= 1850]
            if experiment.startswith("1pct"):
                last_year = int(np.ceil(np.log(4)/np.log(1.01)))
                values_post_1850 = co2_piC * 1.01 ** (np.array(year_cols_post_1850[:last_year +1]).astype(int) - 1850)
                piControl_copy.loc[piControl_copy["Variable"] == co2_var_name, year_cols_post_1850[:last_year +1]] = values_post_1850
                if experiment == "1pctCO2":
                    piControl_copy.loc[:, year_cols_post_1850[last_year +1:]] = np.nan
                else:
                    piControl_copy.loc[piControl_copy["Variable"] == co2_var_name, year_cols_post_1850[last_year:]] = co2_piC * 4
            else:
                mult_factor = experiment.split("abrupt-")[-1].split("x")[0]
                if mult_factor == "0p5":
                    mult_factor = 0.5
                else:
                    mult_factor = int(mult_factor)
                piControl_copy.loc[piControl_copy["Variable"] == co2_var_name, year_cols_post_1850] = co2_piC * mult_factor
            piControl_copy["Scenario"] = experiment
            to_add.append(piControl_copy)
            self.concentrations_df = self.concentrations_df[self.concentrations_df["Scenario"] != experiment].copy()
        self.concentrations_df = pd.concat([self.concentrations_df] + to_add, ignore_index=True)



    def fix_piControl(self, type = "emissions"):
        """Overwrite esm-piControl emissions with 1750 values from historical data.
        
        This method ensures that esm-piControl emissions remain constant at their
        1750 values across all year columns by copying the 1750 value from the
        historical emissions data for each variable.
        """
        if type == "emissions":
            # Filter to esm-piControl scenario
            picontrol_mask = self.emissions_df["Scenario"] == "esm-piControl"
            to_fix = self.emissions_df
            hist = self.emissions_df_hist
        elif type == "concentrations":
            picontrol_mask = self.concentrations_df["Scenario"] == "piControl"
            to_fix = self.concentrations_df
            hist = self.concentrations_df_hist
        # Identify year columns (4-digit names)
        year_cols = [c for c in to_fix.columns 
                        if isinstance(c, str) and re.fullmatch(r"\d{4}", c)]
        
        if not year_cols or "1850" not in year_cols:
            return
        # This gas doesn't have variables in the historical dataset, and the piControl
        # value we have from before is 0, so it is ok to skip it.
        ok_skip = ["Halon1202"]
        # For each variable in piControl, copy 1750 value across all years
        for idx in to_fix[picontrol_mask].index:
            variable = to_fix.loc[idx, "Variable"]
            vari_hist = variable
            if "CO2" in variable and type == "emissions":
                continue
            if vari_hist not in hist["Variable"].values:
                vari_hist_test = variable.split("|")[-1]
                for test in hist["Variable"].values:
                    if test.endswith(vari_hist_test):
                        vari_hist = test
                        break
            # Find matching variable in historical data
            hist_match = hist[
                hist["Variable"] == vari_hist
            ]
            
            if not hist_match.empty and hist_match.iloc[0]["1850"]>= 0:
                value_1850 = hist_match.iloc[0]["1850"]
                # Fill all year columns with the 1850 value
                to_fix.loc[idx, year_cols] = value_1850
            elif variable.split("|")[-1] in ok_skip:
                print(f"Warning: no matching variable {variable} with 1850 value found in historical data for piControl fix, but skipping due to known exception.")
            else:
                print(hist["Variable"].unique())
                print(f"Warning: no matching variable {variable} with 1850 value found in historical data for piControl fix.")
                sys.exit(4)
            to_fix.loc[idx, "Model"] = "emissions_harmonisation_pipeline"
            to_fix.loc[idx, "Mip_Era"] = "CMIP7"
        

    def read_in_previous_rcmip_data(self):
        """Load previous RCMIP CSVs from ``../data`` and apply initial cleaning.

        The method sets the instance attributes ``emissions_df``,
        ``concentrations_df`` and ``forcing_df`` by reading the canonical RCMIP
        CSV files from ``../data``. It applies a small set of deterministic
        cleaning operations expected by the rest of the pipeline:

        - Rename a couple of MAGICC-labelled CO2 variables to match the
          repository naming convention.
        - Remove rows containing the string 'MAGICC' (case-insensitive) from
          the emissions table.
        - Drop concentration columns prior to 1750 (1700--1749) if present.
        - Keep only rows where ``Region == 'World'`` in all three tables.

        Notes
        -----
        This method mutates the instance by adding/updating the attributes
        listed above and does not return a value.
        """
        self.emissions_df = pd.read_csv(os.path.join("..", "data", "rcmip-emissions-annual-means-v5-1-0.csv"))
        self.concentrations_df = pd.read_csv(os.path.join("..", "data", "rcmip-concentrations-annual-means-v5-1-0.csv"))
        self.forcing_df = pd.read_csv(os.path.join("..", "data", "rcmip-radiative-forcing-annual-means-v5-1-0.csv"))
        self.emissions_df.loc[self.emissions_df['Variable']=='Emissions|CO2|MAGICC AFOLU', 'Variable'] = 'Emissions|CO2|AFOLU'
        self.emissions_df.loc[self.emissions_df['Variable']=='Emissions|CO2|MAGICC Fossil and Industrial', 'Variable'] = 'Emissions|CO2|Energy and Industrial Processes'
        # Drop rows in emissions_df where the Variable contains 'MAGIC' (covers MAGICC/MAGGIC),
        # except keep explicitly allowed MAGICC/MAGGIC rows.
        if 'Variable' in self.emissions_df.columns:
            mask_magic = self.emissions_df['Variable'].str.contains('MAGICC', case=False, na=False)
            # keep rows that are not magic-related, or are exceptions
            self.emissions_df = self.emissions_df[~(mask_magic)]
        # Drop columns in concentrations_df with names from 1700 to 1749
        cols_to_drop = [str(year) for year in range(1700, 1750) if str(year) in self.concentrations_df.columns]
        self.concentrations_df = self.concentrations_df.drop(columns=cols_to_drop)
        self.emissions_df = self.emissions_df[self.emissions_df["Region"] == "World"]
        self.concentrations_df = self.concentrations_df[self.concentrations_df["Region"] == "World"]
        self.forcing_df = self.forcing_df[self.forcing_df["Region"] == "World"]
        return
    
    def get_cmip7_natural_forcing(self):
        """Generate CMIP7 natural forcing data from RCMIP forcing data.

        This function assumes that the provided forcing DataFrame has already
        been filtered to only include the 'historical' scenario and have had
        old data cut and metadata added.

        Returns a DataFrame: forcing_df for CMIP7 natural historical.
        """
        filtered_forc = self.forcing_df[(self.forcing_df['Scenario'] == 'historical') &
                                        (self.forcing_df['Variable'].isin([
                                            'Effective Radiative Forcing|Natural',
                                            'Effective Radiative Forcing|Natural|Solar',
                                            'Effective Radiative Forcing|Natural|Volcanic',
                                        ]))].copy()
        filtered_forc = cut_old_rcmip_data_and_add_meta(filtered_forc, 'non-idealised', 1)
        volcanic_df = pd.read_csv("../data/volcanic_forcing_cmip7_minus20tau.csv")
        solar_df = pd.read_csv("../data/solar_forcing_cmip7_extended.csv")
        filtered_forc.loc[filtered_forc['Variable'] == 'Effective Radiative Forcing|Natural|Volcanic',
                          [str(y) for y in range(1750, 2501)]] = volcanic_df["erf"].values[1:]
        filtered_forc.loc[filtered_forc['Variable'] == 'Effective Radiative Forcing|Natural|Solar',
                          [str(y) for y in range(1850, 2501)]] = solar_df["erf"].values[100:]
        filtered_forc.loc[filtered_forc['Variable'] == 'Effective Radiative Forcing|Natural',
                          [str(y) for y in range(1750, 2501)]] = (
            filtered_forc.loc[filtered_forc['Variable'] == 'Effective Radiative Forcing|Natural|Solar',
                              [str(y) for y in range(1750, 2501)]].values +
            filtered_forc.loc[filtered_forc['Variable'] == 'Effective Radiative Forcing|Natural|Volcanic',
                              [str(y) for y in range(1750, 2501)]].values
        )
        return filtered_forc

    
    def get_historical_data_for_cmip7(self):
        """Generate CMIP7 historical data from RCMIP historical emissions and concentrations.

        This function assumes that the provided emissions and concentrations DataFrames
        have already been filtered to only include the 'historical' scenario and have
        had old data cut and metadata added.

        Returns two DataFrames: emissions_df and concentrations_df for CMIP7 historical.
        """
        filtered_em = self.emissions_df[self.emissions_df['Scenario'] == 'historical']
        filtered_em = cut_old_rcmip_data_and_add_meta(filtered_em, "non-idealised", 1)
        filtered_conc = self.concentrations_df[self.concentrations_df['Scenario'] == 'historical']
        filtered_conc = cut_old_rcmip_data_and_add_meta(filtered_conc, 'non-idealised', 1)
        self.emissions_df_hist, self.concentrations_df_hist = get_cmip7_historical_data(
            filtered_em,
            filtered_conc
        )
        return

    
    
    def loop_over_scenarios_and_generate_inputs(self, outpath = "../data/generated_inputs/"):
        """Loop over scenarios from the protocol and generate per-scenario CSV inputs.

        This method iterates the protocol `scenario_info` sheet, generates
        emissions, concentrations and forcing DataFrames per scenario (by
        calling `process_scenario` and other helpers), concatenates all
        results and writes three CSV files to `outpath`.
        """
        emissions_list = []
        concentrations_list = []
        forcing_list = []
        emissions_list_scen7 = []
        concentrations_list_scen7 = []
        forcing_list_scen7 = []
        for index, scenario in self.scenario_df.iterrows():
            # print(scenario)
            # print(f"Processing scenario {index + 1}:")
            if scenario['Scenario']=="esm-??????":
                # print(f"Skipping scenario {scenario['Scenario']} as they are not well-defined")
                continue
            em_df, conc_df, forc_df = self.process_scenario(scenario)
            if em_df is not None:
                if "scen7" in scenario['Scenario']:
                    emissions_list_scen7.append(em_df)
                else:
                    emissions_list.append(em_df)
                if "1700" in em_df.columns:
                    print(f"Warning: scenario {scenario['Scenario']} emissions data contains data from 1700s.")
                    sys.exit(4)
            if conc_df is not None:
                if "scen7" in scenario['Scenario']:
                    concentrations_list_scen7.append(conc_df)
                else:
                    concentrations_list.append(conc_df)
                if "1700" in conc_df.columns:
                    print(f"Warning: scenario {scenario['Scenario']} concentrations data contains data from 1700s.")
                    sys.exit(4)
            if forc_df is not None:
                if "scen7" in scenario['Scenario']: 
                    forcing_list_scen7.append(forc_df)
                else:
                    forcing_list.append(forc_df)

        print("Done looping over scenarios, now concatenating and writing output files...")
        self.post_process_and_write_output(emissions_list, concentrations_list, forcing_list, ftype = "_")
        if len(emissions_list_scen7) > 0:
            self.post_process_and_write_output(emissions_list_scen7, concentrations_list_scen7, forcing_list_scen7, ftype = "_ScenarioMIP_")
        return

    def post_process_and_write_output(self, em_list, conc_list, forc_list, ftype = "_"):
        """Apply final post-processing steps and write the DataFrame to a CSV file.

        This function applies any necessary final transformations to the DataFrame
        (e.g., fixing column names, sorting) and writes it to the specified output
        path as a CSV file.

        Parameters
        ----------
        df : pandas.DataFrame
            The DataFrame to be post-processed and written.
        output_path : str
            The file path where the CSV should be saved.

        Returns
        -------
        None
            The function writes the DataFrame to a CSV file and does not return a value.
        """
        emissions_df = simplify_fgas_and_montreal_gases(pd.concat(em_list, ignore_index=True))
        concentrations_df = pd.concat(conc_list, ignore_index=True)
        forcing_df = pd.concat(forc_list, ignore_index=True)

        if ftype == "_ScenarioMIP_":
            outdir = os.path.join(rcmip_output_dir, "..",  "ScenarioMIP")
        else:
            outdir = rcmip_output_dir
        os.makedirs(outdir, exist_ok=True)        

        for df, fname in zip(
            [emissions_df, concentrations_df, forcing_df],
            [f"rcmip_phase3_emissions{ftype}{self.version}.csv", 
             f"rcmip_phase3_concentrations{ftype}{self.version}.csv", 
             f"rcmip_phase3_forcing{ftype}{self.version}.csv"]
        ):
            if not df.empty:
                df = fill_in_nans(df, self.scenario_df.copy())
                insert_idx = df.columns.get_loc('Mip_Era') + 1
                df.insert(insert_idx, 'Version', f"RCMIP Phase 3 {self.version}")
                # print(df.head())
                # print(df.columns)
                # sys.exit(4)
                df.to_csv(os.path.join(outdir, fname), index=False)    


    def process_scenario(self, scenario):
        """Process one protocol scenario and produce per-table outputs.

        The method inspects a single row from the protocol ``scenario_info``
        DataFrame and attempts to produce three DataFrames for that scenario:
        emissions, concentrations and radiative forcing. The function handles
        many special-case scenario names (e.g. historical variants, flat
        scenarios, scenario7/scen7 inputs and MethaneMIP variants) and falls
        back to selecting exact matches from the previously-loaded RCMIP
        source tables.

        Parameters
        ----------
        scenario : pandas.Series
            One row from the protocol DataFrame containing at least the
            ``Scenario``, ``Type`` and ``Priority`` entries.

        Returns
        -------
        tuple
            A 3-tuple ``(em_df, conc_df, forc_df)`` where each element is a
            ``pandas.DataFrame`` or ``None`` when not applicable. When a
            scenario is fully handled the method returns the prepared
            DataFrames with metadata inserted (``Type``/``Priority``) and the
            scenario name set in the ``Scenario`` column.

        Raises
        ------
        SystemExit
            The method calls ``sys.exit(4)`` for a variety of unrecoverable
            situations (missing source files or unhandled scenario names) to
            preserve original behaviour.
        """
        # Add your scenario processing logic here
        scenario_name = scenario.get('Scenario', 'N/A')
        scenario_type = scenario.get('Type', 'N/A')
        scenario_priority = scenario.get('Priority', 'N/A')
        try:
            scenario_duration = int(scenario.get('Duration of scenario', 0))
        except ValueError:
            scenario_duration = 0
        # print(f"  Scenario name: {scenario_name}")
        # print(f"  Scenario type: {scenario_type}")
        em_df = None
        conc_df = None
        forc_df = None

        # Map of attribute name to holder variable
        mapping = [
            ('emissions_df', 'em_df'),
            ('concentrations_df', 'conc_df'),
            ('forcing_df', 'forc_df'),
        ]
        missing_data = True
        print(scenario_name)
        # Iterate through tables, filter by scenario, and apply helper
        if scenario_name == "historical":
            em_df = self.emissions_df_hist.copy()
            conc_df = self.concentrations_df_hist.copy()
            forc_df = self.get_cmip7_natural_forcing()
            forc_df = add_cmip7_land_use_forcing_to_df(forc_df, scenario_name)
            return em_df, conc_df, forc_df
        elif scenario_name.startswith("hist-"):
            em_df, conc_df, missing_data = self.make_data_for_hist_partial_scenarios(
                    scenario_name,
                    scenario_type,
                    scenario_priority
                )
            return em_df, conc_df, None
        elif any(sub in scenario_name for sub in ["lowNTCF", "ssp370-lowCH4"]):
            em_df, conc_df, missing_data = self.make_lowntcf_scenarios(
                scenario_name,
                scenario_type,
                scenario_priority
            )
            return em_df, conc_df, None
        elif scenario_name in ["esm-allGHG-ssp534-over-highCH4", "esm-allGHG-ssp585-lowCH4"]:
            em_df, conc_df, missing_data = self.make_different_methane_scenario(
                scenario_name,
                scenario_type,
                scenario_priority
            )
            return em_df, conc_df, None
        if "scen7" in scenario_name:
            if os.path.exists(self.cmip7_scenario_file_path):
                em_df, conc_df, forc_df, missing_data = self.make_scenario7_data(
                    scenario_name,
                    scenario_type,
                    scenario_priority
                )
                return em_df, conc_df, forc_df
        if "methanemip" in scenario_name:
            return self.make_methanemip_scenario_data(
                scenario_name,
                scenario_type,
                scenario_priority
            ), None, None

        for src_attr, _dest in mapping:
            src_df = getattr(self, src_attr)
            # Normal case: exact scenario present in the source table
            if scenario_name in src_df['Scenario'].unique():
                filtered = src_df[src_df['Scenario'] == scenario_name]
                new_df = cut_old_rcmip_data_and_add_meta(filtered, scenario_type, scenario_priority)
                if src_attr == 'emissions_df':
                    em_df = new_df
                    missing_data = False
                elif src_attr == 'concentrations_df':
                    conc_df = new_df
                    if scenario_name == "1pctCO2":
                        years_nan = [str(y) for y in range(1750+scenario_duration,2501)]
                        #print(years_nan)
                        conc_df.loc[:, years_nan] = np.nan
                        #print(conc_df.loc[:, "1980":"1998"].head())
                        #sys.exit(4)
                    missing_data = False
                elif src_attr == 'forcing_df':
                    forc_df = new_df
                    missing_data = False
            elif scenario_name.startswith("esm-allGHG") and "".join(scenario_name.split("-allGHG")) in self.scenario_df['Scenario'].unique():
                missing_data = False
            elif scenario_name.startswith("esm-ssp") and scenario_name.split("esm-")[-1] in self.scenario_df['Scenario'].unique():
                missing_data = False
            elif scenario_name.startswith("esm-allGHG-ssp") and scenario_name.split("esm-allGHG-")[-1] in self.scenario_df['Scenario'].unique():
                missing_data = False
            # Special-case: historical-cmip6 uses the historical table but with renamed Scenario
            elif scenario_name == "historical-cmip6":
                if src_attr == 'emissions_df':
                    filtered = self.emissions_df[self.emissions_df['Scenario'] == 'historical']
                    em_df = cut_old_rcmip_data_and_add_meta(filtered, scenario_type, scenario_priority)
                    em_df["Scenario"] = "historical-cmip6"
                    missing_data = False
                elif src_attr == 'concentrations_df':
                    filtered = self.concentrations_df[self.concentrations_df['Scenario'] == 'historical']
                    conc_df = cut_old_rcmip_data_and_add_meta(filtered, scenario_type, scenario_priority)
                    conc_df["Scenario"] = "historical-cmip6"
                    missing_data = False
                elif src_attr == 'forcing_df':
                    filtered = self.forcing_df[self.forcing_df['Scenario'] == 'historical']
                    forc_df = cut_old_rcmip_data_and_add_meta(filtered, scenario_type, scenario_priority)
                    forc_df["Scenario"] = "historical-cmip6"
                    missing_data = False

        if missing_data and scenario_type == "idealised":
            em_df, conc_df, forc_df, missing_data = self.make_data_for_idealised_scenario(scenario_name)


        data_from_other = ["esm-hist", "esm-hist-allGHG", "esm-hist-cmip6"]
        if missing_data and scenario_name not in data_from_other:
            if "scen7" in scenario_name and not os.path.exists(self.cmip7_scenario_file_path):
                print(f"Warning: no data source file found for scenario {scenario_name} at {self.cmip7_scenario_file_path}.")
                return None, None, None
            print(f"No data found for this scenario: {scenario_name}.")
            sys.exit(4)
        return em_df, conc_df, forc_df
    
    def make_data_for_idealised_scenario(self, scenario):
        """Create emissions/concentrations/forcing for idealised scenarios.

        This handler implements a small set of idealised scenarios required by
        the protocol. The current implementation supports the special-case
        scenario ``1pctCO2-cdr`` (constructed by combining and reversing
        portions of other concentration tables) and a range of simple
        placeholders listed in the code. Flat scenarios are delegated to
        ``make_flat_scenario_emissions``.

        Parameters
        ----------
        scenario : str
            Scenario name (protocol ``Scenario`` value).

        Returns
        -------
        tuple
            Either a 4-tuple ``(em_df, conc_df, forc_df, missing_data)`` where
            each element is a ``pandas.DataFrame`` or ``None`` when not
            applicable, and ``missing_data`` is a boolean flag. For
            unhandled scenarios the function will call ``sys.exit(4)``.
        """
        if scenario == "1pctCO2-cdr":
            conc_values_abrupt = self.concentrations_df[
                (self.concentrations_df['Scenario'] == 'abrupt-4xCO2') &
                (self.concentrations_df['Region'] == 'World')
            ]
            conc_values_1pct4ext = self.concentrations_df[
                (self.concentrations_df['Scenario'] == '1pctCO2-4xext') &
                (self.concentrations_df['Region'] == 'World')
            ]
            conc_values_picontrol = self.concentrations_df[
                (self.concentrations_df['Scenario'] == 'piControl') &
                (self.concentrations_df['Region'] == 'World')
            ]
            conc_data = conc_values_1pct4ext.copy()
            conc_data.loc[:, 'Scenario'] = '1pctCO2-cdr'
            way_up = True
            rev_idx = 0
            co2_reverse = []
            for year in range(1850, 2501):
                co2_conc = conc_data.loc[conc_data['Variable'] == 'Atmospheric Concentrations|CO2'][f'{year}'].values
                co2_4x = conc_values_abrupt.loc[conc_values_abrupt['Variable'] == 'Atmospheric Concentrations|CO2'][f'{year}'].values

                if way_up and co2_conc < co2_4x:
                    # Still in increasing phase
                    co2_reverse.append(co2_conc)
                elif way_up and co2_conc == co2_4x:
                    # At the peak, switch direction next time
                    #co2_reverse.append(co2_conc)
                    way_up = False
                else:
                    # Now reverse towards piControl values
                    rev_idx += 1
                    #co2_reverse[-(rev_idx)]
                    if rev_idx < len(co2_reverse):
                        conc_data.loc[conc_data['Variable'] == 'Atmospheric Concentrations|CO2', f'{year}'] = co2_reverse[-(rev_idx)]
                    else:
                        conc_data.loc[conc_data['Variable'] == 'Atmospheric Concentrations|CO2', f'{year}'] = conc_values_picontrol.loc[conc_values_picontrol['Variable'] == 'Atmospheric Concentrations|CO2'][f'{year}'].values
            #print(conc_data.loc[conc_data['Variable'] == 'Atmospheric Concentrations|CO2', "1850":"2100"])
            #sys.exit(4)
            conc_df = cut_old_rcmip_data_and_add_meta(conc_data, "idealised", 1)
            conc_df["Scenario"] = "1pctCO2-cdr"
            return None, conc_df, None, False
        if scenario in ["1pctCO2-bgc", "1pctCO2-rad", "esm-1pct-brch-1000PgC", "esm-1pct-brch-2000PgC", "esm-1pct-brch-750PgC"]:
            return None, None, None, False
        if "flat" in scenario:
            em_df = self.make_flat_scenario_emissions(scenario)
            return em_df, None, None, False
        print("Idealised scenario handling not supported for scenario:", scenario)
        sys.exit(4)
        return None, None, None, True

    def make_data_for_hist_partial_scenarios(self, scenario, scenario_type, scenario_priority):
        """Construct data for historical-partial scenarios.

        The historical-partial handlers implement variants of the historical
        scenario where subsets of variables are taken from historical data and
        others taken from piControl (or vice versa). Supported sub-scenarios
        include ``hist-aer``, ``hist-CO2`` and ``hist-GHG``.

        Parameters
        ----------
        scenario : str
            The protocol scenario name (e.g. ``hist-aer``, ``hist-CO2``).
        scenario_type : str
            The ``Type`` metadata value to insert via
            ``cut_old_rcmip_data_and_add_meta``.
        scenario_priority : int
            The ``Priority`` metadata value to insert.

        Returns
        -------
        tuple
            ``(em_df, conc_df, missing_data)`` where ``em_df`` and ``conc_df``
            are the constructed emissions and concentrations DataFrames and
            ``missing_data`` is a boolean flag (``False`` on success).

        Raises
        ------
        SystemExit
            For unsupported sub-scenarios the function exits with code 4 to
            preserve original behaviour.
        """
        em_df_start = self.emissions_df_hist.copy()
        conc_df_start = self.concentrations_df_hist.copy()
        em_picontrol = self.emissions_df[
            (self.emissions_df['Scenario'] == 'esm-piControl') &
            (self.emissions_df['Region'] == 'World')
        ].copy()
        conc_picontrol = self.concentrations_df[
            (self.concentrations_df['Scenario'] == 'piControl') &
            (self.concentrations_df['Region'] == 'World')
        ].copy()
        aerosol_vars = [
            'Emissions|BC',
            'Emissions|OC',
            'Emissions|Sulfur',
            'Emissions|NOx',
            'Emissions|NH3'
        ]
        other_non_GHG =[
            'Emissions|CO',
            'Emissions|VOC',
        ]
        years_nan = [str(year) for year in range(2023, 2501)]
        if scenario == "hist-aer":
            # Take all but aerosols from piControl, add aerosols from em_df_start
            em_df = em_picontrol[~em_picontrol['Variable'].isin(aerosol_vars + other_non_GHG)].copy()
            for var in (aerosol_vars + other_non_GHG):
                new_row = em_df_start[em_df_start['Variable'] == var].copy()
                if not new_row.empty:
                    em_df = pd.concat([em_df, new_row], ignore_index=True)
            em_df.loc[:, years_nan] = np.nan
            em_df = cut_old_rcmip_data_and_add_meta(em_df, scenario_type, scenario_priority)
            conc_df = conc_picontrol.copy()
            conc_df.loc[:, years_nan] = np.nan
            conc_df = cut_old_rcmip_data_and_add_meta(conc_df, scenario_type, scenario_priority)

        elif scenario == "hist-CO2":
            # Take all but CO2 from picControl, add CO2 from em_df_start
            co2_variables = ['Emissions|CO2', "Emissions|CO2|AFOLU", "Emissions|CO2|Energy and Industrial Processes"]
            em_df = em_picontrol[~em_picontrol['Variable'].isin(co2_variables)].copy()
            for var in co2_variables:
                new_row = em_df_start[em_df_start['Variable'] == var].copy()
                if not new_row.empty:
                    em_df = pd.concat([em_df, new_row], ignore_index=True)
            em_df.loc[:, years_nan] = np.nan
            em_df = cut_old_rcmip_data_and_add_meta(em_df, scenario_type, scenario_priority)
            #print()
            conc_df = conc_picontrol[~(conc_picontrol['Variable']=='Atmospheric Concentrations|CO2')].copy()
            new_row = conc_df_start[conc_df_start['Variable']=='Atmospheric Concentrations|CO2'].copy()
            conc_df = pd.concat([conc_df, new_row], ignore_index=True)
            conc_df.loc[:, years_nan] = np.nan
            conc_df = cut_old_rcmip_data_and_add_meta(conc_df, scenario_type, scenario_priority)
        elif scenario == "hist-GHG":
            # Take  non-aerosols from 

            em_df = em_df_start[~em_df_start['Variable'].isin(aerosol_vars + other_non_GHG)].copy()
            for var in aerosol_vars + other_non_GHG:
                new_row = em_picontrol[em_picontrol['Variable'] == var].copy()
                if not new_row.empty:
                    em_df = pd.concat([em_df, new_row], ignore_index=True)
            em_df.loc[:, years_nan] = np.nan
            em_df = cut_old_rcmip_data_and_add_meta(em_df, scenario_type, scenario_priority)

            conc_df = conc_df_start.copy()
            conc_df = cut_old_rcmip_data_and_add_meta(conc_df, scenario_type, scenario_priority)
        
        else:
            print("hist-partial handling not supported for scenario:", scenario)
            sys.exit(4)
        conc_df["Scenario"] = scenario
        em_df["Scenario"] = scenario 
        # print(em_df.columns)
        return em_df, conc_df, False

    def make_lowntcf_scenarios(self, scenario_name, scenario_type, scenario_priority):
        """Create datasets for low-NTCF methane and related scenarios.

        The function builds scenario variants where methane (CH4) and
        non-CO2 tropospheric forcing (NTCF) adjustments are applied by mixing
        or replacing rows from existing SSP scenarios. The produced DataFrames
        include the requested ``Type`` and ``Priority`` metadata values.

        Parameters
        ----------
        scenario_name : str
            Protocol scenario name to generate (e.g. ``esm-allGHG-ssp370-lowNTCF``).
        scenario_type : str
            The ``Type`` metadata value to insert via
            ``cut_old_rcmip_data_and_add_meta``.
        scenario_priority : int
            The ``Priority`` metadata value to insert.

        Returns
        -------
        tuple
            ``(final_em_df, final_conc_df, False)`` on success. On failure the
            function will call ``sys.exit(4)``.
        """
        filtered_em = self.emissions_df[self.emissions_df['Scenario'] == 'ssp370-lowNTCF-gidden'].copy()
        new_df_em = cut_old_rcmip_data_and_add_meta(filtered_em, scenario_type, scenario_priority)
        filtered_conc = self.concentrations_df[self.concentrations_df['Scenario'] == 'ssp370-lowNTCF-gidden'].copy()
        new_df_conc = cut_old_rcmip_data_and_add_meta(filtered_conc, scenario_type, scenario_priority)
        if scenario_name == "esm-allGHG-ssp370-lowNTCF":
            new_df_em["Scenario"] = scenario_name
            new_df_conc["Scenario"] = scenario_name
            return new_df_em, new_df_conc, False
        filtered_em_base = self.emissions_df[self.emissions_df['Scenario'] == 'ssp370'].copy()
        base_df_em = cut_old_rcmip_data_and_add_meta(filtered_em_base, scenario_type, scenario_priority)
        filtered_conc_base = self.concentrations_df[self.concentrations_df['Scenario'] == 'ssp370'].copy()
        base_df_conc = cut_old_rcmip_data_and_add_meta(filtered_conc_base, scenario_type, scenario_priority)
        if scenario_name == "esm-allGHG-ssp370-lowCH4":
            # Replace CH4 in base with lowNTCF CH4
            final_em_df = base_df_em[base_df_em['Variable'] != 'Emissions|CH4'].copy()
            ch4_lowntcf = new_df_em[new_df_em['Variable'] == 'Emissions|CH4'].copy()
            final_em_df = pd.concat([final_em_df, ch4_lowntcf], ignore_index=True)
            final_conc_df = base_df_conc[base_df_conc['Variable'] != 'Atmospheric Concentrations|CH4'].copy()
            ch4_lowntcf_conc = new_df_conc[new_df_conc['Variable'] == 'Atmospheric Concentrations|CH4'].copy()
            final_conc_df = pd.concat([final_conc_df, ch4_lowntcf_conc], ignore_index=True)
        elif scenario_name == "esm-allGHG-ssp370-lowNTCF-HighCH4":
            # Replace all but CH4 in base with lowNTCF data
            final_em_df = new_df_em[new_df_em['Variable'] != 'Emissions|CH4'].copy()
            ch4_base = base_df_em[base_df_em['Variable'] == 'Emissions|CH4'].copy()
            final_em_df = pd.concat([final_em_df, ch4_base], ignore_index=True)
            final_conc_df = new_df_conc[new_df_conc['Variable'] != 'Atmospheric Concentrations|CH4'].copy()
            ch4_base_conc = base_df_conc[base_df_conc['Variable'] == 'Atmospheric Concentrations|CH4'].copy()
            final_conc_df = pd.concat([final_conc_df, ch4_base_conc], ignore_index=True)
        else:
            print("LowNTCF handling not supported for scenario:", scenario_name)
            sys.exit(4)
        final_conc_df["Scenario"] = scenario_name
        final_em_df["Scenario"] = scenario_name       
        return final_em_df, final_conc_df, False
    
    def make_different_methane_scenario(self, scenario_name, scenario_type, scenario_priority):
        """Produce emissions and concentrations where methane is swapped.

        This helper creates a scenario where the CH4 component from one SSP is
        swapped into another SSP (for example to create high-CH4/low-CH4
        variants). Other variables are preserved from the base SSP.

        Parameters
        ----------
        scenario_name : str
            Protocol scenario name identifying which methane swap to perform.
        scenario_type : str
            ``Type`` metadata to insert into resulting DataFrames.
        scenario_priority : int
            ``Priority`` metadata to insert.

        Returns
        -------
        tuple
            ``(final_em_df, final_conc_df, False)`` with the constructed
            DataFrames. The function exits with code 4 for unexpected inputs.
        """
        filtered_em_ssp534 = self.emissions_df[self.emissions_df['Scenario'] == 'ssp534-over'].copy()
        filtered_em_ssp534 = cut_old_rcmip_data_and_add_meta(filtered_em_ssp534, scenario_type, scenario_priority)
        filtered_conc_ssp534 = self.concentrations_df[self.concentrations_df['Scenario'] == 'ssp534-over'].copy()
        filtered_conc_ssp534 = cut_old_rcmip_data_and_add_meta(filtered_conc_ssp534, scenario_type, scenario_priority)
        filtered_em_ssp585 = self.emissions_df[self.emissions_df['Scenario'] == 'ssp585'].copy()
        filtered_em_ssp585 = cut_old_rcmip_data_and_add_meta(filtered_em_ssp585, scenario_type, scenario_priority)
        filtered_conc_ssp585 = self.concentrations_df[self.concentrations_df['Scenario'] == 'ssp585'].copy()
        filtered_conc_ssp585 = cut_old_rcmip_data_and_add_meta(filtered_conc_ssp585, scenario_type, scenario_priority)
        if scenario_name == "esm-ssp534-over-highCH4-allGHG":
            # Replace CH4 in ssp534-over with ssp585 CH4
            final_em_df = filtered_em_ssp534[filtered_em_ssp534['Variable'] != 'Emissions|CH4'].copy()
            ch4_ssp585 = filtered_em_ssp585[filtered_em_ssp585['Variable'] == 'Emissions|CH4'].copy()
            final_em_df = pd.concat([final_em_df, ch4_ssp585], ignore_index=True)
            final_conc_df = filtered_conc_ssp534[filtered_conc_ssp534['Variable'] != 'Atmospheric Concentrations|CH4'].copy()
            ch4_ssp585_conc = filtered_conc_ssp585[filtered_conc_ssp585['Variable'] == 'Atmospheric Concentrations|CH4'].copy()
            final_conc_df = pd.concat([final_conc_df, ch4_ssp585_conc], ignore_index=True)
        else:
           # Replace CH4 in ssp585 with ssp534-over CH4
            final_em_df = filtered_em_ssp585[filtered_em_ssp585['Variable'] != 'Emissions|CH4'].copy()
            ch4_ssp534 = filtered_em_ssp534[filtered_em_ssp534['Variable'] == 'Emissions|CH4'].copy()
            final_em_df = pd.concat([final_em_df, ch4_ssp534], ignore_index=True)
            final_conc_df = filtered_conc_ssp585[filtered_conc_ssp585['Variable'] != 'Atmospheric Concentrations|CH4'].copy()
            ch4_ssp534_conc = filtered_conc_ssp534[filtered_conc_ssp534['Variable'] == 'Atmospheric Concentrations|CH4'].copy()
            final_conc_df = pd.concat([final_conc_df, ch4_ssp534_conc], ignore_index=True)            
        final_conc_df["Scenario"] = scenario_name
        final_em_df["Scenario"] = scenario_name
        return final_em_df, final_conc_df, False
    
    def cut_year_columns_beyond_duration(self, df, scenario_name, set_zeros_in_duration= True):
        """Set year columns beyond scenario duration to NaN.

        This helper identifies the year columns in the provided DataFrame and
        sets values to NaN for any years that exceed the duration specified in
        the protocol for the given scenario.

        Parameters
        ----------
        df : pandas.DataFrame
            The DataFrame containing year columns to be processed.
        scenario_name : str
            The name of the scenario to look up duration for.

        Returns
        -------
        pandas.DataFrame
            A copy of the input DataFrame with year columns beyond the scenario
            duration set to NaN.
        """
        duration = self.scenario_df.loc[self.scenario_df['Scenario'] == scenario_name, 'Duration of scenario'].values[0]
        year_cols = [str(year) for year in range(1750 + duration, 2501)]

        df.loc[:, year_cols] = np.nan
        if set_zeros_in_duration:
            year_cols_in_duration = [str(year) for year in range(1750, 1750 + duration)]
            df.loc[:, year_cols_in_duration] = 0.0
        return df
    
    def make_flat_scenario_emissions(self, scenario_name):
        """Generate emissions DataFrame for 'flat' idealised scenarios.

        Flat scenarios are constructed by taking the ``esm-piControl`` CO2
        components (AFOLU and Energy/FFI) and overwriting them with values
        defined in the ``flat10_inputs.csv`` table. The CSV may contain year
        columns with float-like names which are normalised into integer year
        columns by this function.

        Parameters
        ----------
        scenario_name : str
            Protocol scenario name starting with ``esm-flat`` that encodes the
            multiplicative factor and variant (e.g. ``esm-flat10``).

        Returns
        -------
        pandas.DataFrame
            Emissions DataFrame prepared for the scenario with ``Scenario``
            set and CO2 total rows populated.

        Raises
        ------
        SystemExit
            If required mapping rows are missing or an unsupported flat
            scenario variant is encountered.
        """
        filtered_em_base = self.emissions_df[self.emissions_df['Scenario'] == "esm-piControl"].copy()
        filtered_em_base = cut_old_rcmip_data_and_add_meta(filtered_em_base, "idealised", 1)
        filtered_em_base = self.cut_year_columns_beyond_duration(filtered_em_base, scenario_name)
        idxs_keep = []
        for idx,row in filtered_em_base.iterrows():
            if "Emissions|CO2" in row.Variable:
                idxs_keep.append(idx)
        filtered_em_base = filtered_em_base.loc[idxs_keep,:]
        filtered_em_base = filtered_em_base.loc[filtered_em_base['Variable'].str.contains("Emissions|CO2")]
  
        flat_em_df = pd.read_csv(os.path.join("..", "data", "flat10_inputs.csv"))
        # Rename columns that start with four digits to integer column names
        rename_map = {}
        for col in list(flat_em_df.columns):
            m = re.match(r"^(\d{4})", str(col))
            if m:
                rename_map[col] = int(m.group(1))
        if rename_map:
            flat_em_df = flat_em_df.rename(columns=rename_map)

        flat_em_df.loc[flat_em_df['variable']=='Emissions|CO2|MAGICC Fossil and Industrial', 'variable'] = 'Emissions|CO2|Energy and Industrial Processes'
        flat_em_df.loc[flat_em_df['variable']=='Emissions|CO2|MAGICC AFOLU', 'variable'] = 'Emissions|CO2|AFOLU'

        factor = scenario_name.split("esm-flat")[-1].split("-")[0]

        match_scenario = scenario_name[4:].replace("-", "_").replace(factor, "10")
        factor = float(factor)

        flat_em_df_now =flat_em_df[flat_em_df["scenario"]==match_scenario].copy()
        year_cols = [c for c in flat_em_df_now.columns if isinstance(c, int) and c >= 1750]
        year_cols_str = [str(c) for c in year_cols]
        if not flat_em_df_now.empty:
            for idx, row in flat_em_df_now.iterrows():
                var = row['variable']
                if var in filtered_em_base['Variable'].values:
                    filtered_em_base.loc[filtered_em_base['Variable'] == var, year_cols_str] = row.loc[year_cols].values * 1.e3 / MOL_WEIGHTS["C"] * MOL_WEIGHTS["CO2"]*factor/10.0
                else:
                    print(f"Warning: variable {var} not found in esm-piControl data for flat scenario {scenario_name}.")
                    sys.exit(4)

        elif match_scenario == "flat10_nz":
            flat_years = [str(year) for year in range(1850, 1950)]
            linear_reduction_years = [str(year) for year in range(1950, 2000)]
            linear_values = np.linspace(10.0, 0.0, len(linear_reduction_years))
            filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2|Energy and Industrial Processes", flat_years] = 10.0 * 1.e3 / MOL_WEIGHTS["C"] * MOL_WEIGHTS["CO2"]*factor/10
            filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2|Energy and Industrial Processes", linear_reduction_years] = linear_values* 1.e3 / MOL_WEIGHTS["C"] * MOL_WEIGHTS["CO2"]*factor/10
            filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2|Energy and Industrial Processes", [str(year) for year in range(2000, 2150)]] = 0.0
            filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2|AFOLU", [str(year) for year in range(1750, 2150)]] = 0.0
            filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2", [str(year) for year in range(1750, 2150)]] = filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2|Energy and Industrial Processes", [str(year) for year in range(1750, 2150)]]
        elif match_scenario == "flat10_rev":
            flat_years = [str(year) for year in range(1850, 1950)]
            zec_years = [str(year) for year in range(1950, 2050)]
            flat_neg_years = [str(year) for year in range(2050, 2150)]
            filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2|Energy and Industrial Processes", flat_years] = 10.0 * 1.e3 / MOL_WEIGHTS["C"] * MOL_WEIGHTS["CO2"]*factor/10
            filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2|Energy and Industrial Processes", zec_years] = 0.0
            filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2|Energy and Industrial Processes", flat_neg_years] = -10.0 * 1.e3 / MOL_WEIGHTS["C"] * MOL_WEIGHTS["CO2"]*factor/10
            filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2|Energy and Industrial Processes", str(2150)] = 0.0
            filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2|AFOLU", [str(year) for year in range(1750, 2150)]] = 0.0
            filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2", [str(year) for year in range(1750, 2150)]] = filtered_em_base.loc[filtered_em_base['Variable'] == "Emissions|CO2|Energy and Industrial Processes", [str(year) for year in range(1750, 2150)]]
        else:
            print(f"Warning: no data found for flat scenario {scenario_name}.")
            sys.exit(4)

        filtered_em_base["Scenario"] = scenario_name
        filtered_em_base.loc[filtered_em_base['Variable']=='Emissions|CO2', year_cols_str] = filtered_em_base.loc[filtered_em_base['Variable']=='Emissions|CO2|Energy and Industrial Processes', year_cols_str].values
        return filtered_em_base
    
    def make_scenario7_data(self, scenario_name, scenario_type, scenario_priority, scen_file = None):
        """Map Scenario7/SCEN7-style inputs to RCMIP emissions rows.

        The function reads the CMIP7-style scenario CSV indicated by
        ``scen_file`` (or the instance ``cmip7_scenario_file_path``) and
        selects a model-specific subset of rows for the requested short
        scenario name. It then maps the CMIP7 variable names onto the RCMIP
        historical template using ``map_in_cmip7_style_emissions_to_rcmip_format``
        and inserts metadata via ``cut_old_rcmip_data_and_add_meta``.

        Parameters
        ----------
        scenario_name : str
            The protocol scenario name (e.g. ``esm-M``), where the short code
            is mapped to a full scenario/model pair internally.
        scenario_type : str
            ``Type`` metadata value to insert.
        scenario_priority : int
            ``Priority`` metadata value to insert.
        scen_file : str, optional
            Path to the CMIP7-style scenario CSV. Defaults to the instance
            attribute when ``None``.

        Returns
        -------
        tuple
            ``(em_df, None, None, False)`` where ``em_df`` is the prepared
            emissions DataFrame. The function exits with code 4 for
            unrecognised short names or missing data.
        """
        if scenario_name.startswith("esm-") and scenario_name not in ["esm-allGHG-scen7-H-CH4L", "esm-allGHG-scen7-L-CH4H"]:
            return None, None, None, False
        if scen_file is None:
            scen_file = self.cmip7_scenario_file_path
        if not os.path.exists(scen_file):
            return None, None, None, False
        
        scen7_inputs = fix_columns_to_int_string(pd.read_csv(scen_file))
        scen7_history = self.emissions_df_hist.copy()
        if "CH4" not in scenario_name:
            scen_name_short = scenario_name.split("-")[-1].split("C")[0]
            scen_name_full, scen_model = get_scenarioinfor_from_dict(scen_name_short, scenario_name)
            em_df, model_full = get_scen7_subframe_just_one(
                scen7_inputs, 
                scen_name_full, 
                scen_model, 
                scenario_name, 
                scenario_type, 
                scenario_priority, 
                scen7_history
                )
            em_df.loc[:,"Scenario"] = scenario_name.split("C")[0]
            em_df.loc[:,"Model"] = model_full
            em_df.loc[:,"Mip_Era"] = "CMIP7"
            conc_df = self.get_cmip7_scenario_concentrations(scenario_name, scen_name_short)
            if conc_df is not None:
                conc_df.loc[:,"Scenario"] = scenario_name.split("C")[0]
                conc_df.loc[:,"Model"] = model_full
                conc_df.loc[:,"Mip_Era"] = "CMIP7"

            forc_df = self.get_cmip7_natural_forcing()
            forc_df.loc[:,"Scenario"] = scenario_name.split("C")[0]
            forc_df.loc[:,"Mip_Era"] = "CMIP7"
            forc_df = add_cmip7_land_use_forcing_to_df(forc_df, scenario_name.split("C")[0])
            forc_df.loc[:,"Model"] = model_full
            return em_df, conc_df, forc_df, False
        print("Got into the CH4 handling for scenario7")
        scen_short1 = scenario_name.split("-")[-2]
        scen_short2 = scenario_name[-1]
        em_dfs = []
        models_full = []
        for scen_short in [scen_short1, scen_short2]:
            scen_name_full, scen_model = get_scenarioinfor_from_dict(scen_short, scenario_name)
            em_df, model_full = get_scen7_subframe_just_one(
                    scen7_inputs, 
                    scen_name_full, 
                    scen_model, 
                    scenario_name, 
                    scenario_type, 
                    scenario_priority, 
                    scen7_history
                    )
            em_dfs.append(em_df)
            models_full.append(model_full)
        em_df = em_dfs[0][em_dfs[0]["Variable"] != "Emissions|CH4"]
        meth_row = em_dfs[1][em_dfs[1]["Variable"] == "Emissions|CH4"]
        em_df = pd.concat([em_df, meth_row], ignore_index=True)
        em_df.loc[:,"Scenario"] = scenario_name
        em_df.loc[:,"Model"] = f"{models_full[0]}/{models_full[1]}"
        forc_df = self.get_cmip7_natural_forcing()
        forc_df.loc[:,"Scenario"] = scenario_name
        forc_df.loc[:,"Mip_Era"] = "CMIP7"
        forc_df = add_cmip7_land_use_forcing_to_df(forc_df, f"scen7-{scen_short1}")
        forc_df.loc[:,"Scenario"] = scenario_name
        forc_df.loc[:,"Model"] = models_full[0]
        return em_df, None, forc_df, False
        
        

    
    def make_methanemip_scenario_data(self, scenario_name, scenario_type, scenario_priority):
        """Create emissions/concentrations for MethaneMIP-derived scenarios.

        The function looks up a mapping from protocol scenario names to
        MethaneMIP scenario identifiers, extracts the relevant subframe from
        the previously restructured MethaneMIP DataFrame and inserts metadata
        before returning the DataFrame prepared for the scenario.

        Parameters
        ----------
        scenario_name : str
            Protocol scenario name (one of the keys in the internal mapping).
        scenario_type : str
            ``Type`` metadata value to insert.
        scenario_priority : int
            ``Priority`` metadata value to insert.

        Returns
        -------
        pandas.DataFrame
            The prepared MethaneMIP DataFrame for the requested scenario.

        Raises
        ------
        SystemExit
            If no data is found for the requested MethaneMIP mapping.
        """
        methanemip_mapping = {
            "methanemip-TM-allGHG": "ssp245-ModAggr",
            "methanemip-TM+BC-allGHG": "ssp245-Aggr",
        }
        methanemip_subframe = self.methanemip_df.loc[self.methanemip_df['Scenario']==methanemip_mapping[scenario_name]].copy()
        if methanemip_subframe.empty:
            print(f"Warning: no data found for methanemip scenario {scenario_name}.")
            sys.exit(4)
        methanemip_subframe.loc[:, 'Scenario'] = scenario_name
        unique_couples = self.emissions_df_hist[["Variable", "Unit"]].drop_duplicates()

        methanemip_subframe = fix_fgases_methanemip(methanemip_subframe, self.emissions_df_hist["Variable"].unique())
        methanemip_subframe = fix_units_methanemip(methanemip_subframe, unique_couples)
        methanemip_subframe = cut_old_rcmip_data_and_add_meta(methanemip_subframe, scenario_type, scenario_priority)
        return methanemip_subframe

    def get_cmip7_scenario_concentrations(self, scenario_name, scen_name_short):
        """
        Extract concentrations data for CMIP7 scenarios from esgf files.

        Returns
        -------
        pandas.DataFrame
            A concentrations DataFrame prepared for the CMIP7 scenario. Currently
            this is just a copy of the historical concentrations with metadata.
        """
        conc_df_start = self.concentrations_df_hist.copy()
        scen_name_file = scen_name_short.split("-")[-1].split("C")[0].lower()
        scen_file = os.path.join("..", "data", f"{scen_name_file}_output_concentrations.csv")
        if not os.path.exists(scen_file):
            print(f"Warning: no concentrations file found for scenario7 scenario {scenario_name} at {scen_file}.")
            print(scen_file)
            print(scen_name_file)
            sys.exit(4)
            return None
        scen_data = fix_columns_to_int_string(pd.read_csv(scen_file))
        conc_df = merge_cmip7_concentrations_data(conc_df_start, scen_data)
        return cut_old_rcmip_data_and_add_meta(conc_df, "scenario7", 1)

if __name__ == "__main__":
    version = "v2.0.0"
    #read_and_restructure_methanemip_data()
    rcmipdatagen = RcmipInputGenerator(
        excel_file_path=os.path.join(rcmip_output_dir, f"rcmip_phase3_protocol_{version}.xlsx"),
        version=version
    )
    rcmipdatagen.loop_over_scenarios_and_generate_inputs()
    print("Done looping")
    check_experiment_contents_of_generated_data(rcmipdatagen.scenario_df, version= version)