def run_supy(
        df_forcing: pandas.DataFrame,
        df_state_init: pandas.DataFrame,
        save_state=False,
        n_yr=10,
)->Tuple[pandas.DataFrame, pandas.DataFrame]:
    '''Perform supy simulation.

    Parameters
    ----------
    df_forcing : pandas.DataFrame
        forcing data for all grids in `df_state_init`.
    df_state_init : pandas.DataFrame
        initial model states;
        or a collection of model states with multiple timestamps, whose last temporal record will be used as the initial model states.
    save_state : bool, optional
        flag for saving model states at each time step, which can be useful in diagnosing model runtime performance or performing a restart run.
        (the default is False, which instructs supy not to save runtime model states).
    n_yr : int, optional
        chunk size (`n_yr` years) to split simulation periods so memory usage can be reduced.
        (the default is 10, which implies 10-year forcing chunks used in simulations).

    Returns
    -------
    df_output, df_state_final : Tuple[pandas.DataFrame, pandas.DataFrame]
        - df_output: `output results <df_output_var>`
        - df_state_final: `final model states <df_state_var>`

    Examples
    --------

    >>> df_output, df_state_final = supy.run_supy(df_forcing, df_state_init)


    '''

    # save df_init without changing its original data
    # df.copy() in pandas does work as a standard python deepcopy
    df_init = df_state_init.copy()

    # retrieve the last temporal record as `df_init`
    # if a `datetime` level existing in the index
    if df_init.index.nlevels > 1:
        idx_dt = df_init.index.get_level_values('datetime').unique()
        dt_last = idx_dt.max()
        df_init = df_init.loc[dt_last]
    # add placeholder variables for df_forcing
    # `metforcingdata_grid` and `ts5mindata_ir` are used by AnOHM and ESTM, respectively
    # they are now temporarily disabled in supy
    df_forcing = df_forcing\
        .assign(
            metforcingdata_grid=0,
            ts5mindata_ir=0,
        )\
        .rename(
            # remanae is a workaround to resolve naming inconsistency between
            # suews fortran code interface and input forcing file hearders
            columns={
                '%' + 'iy': 'iy',
                'id': 'id',
                'it': 'it',
                'imin': 'imin',
                'qn': 'qn1_obs',
                'qh': 'qh_obs',
                'qe': 'qe',
                'qs': 'qs_obs',
                'qf': 'qf_obs',
                'U': 'avu1',
                'RH': 'avrh',
                'Tair': 'temp_c',
                'pres': 'press_hpa',
                'rain': 'precip',
                'kdown': 'avkdn',
                'snow': 'snow_obs',
                'ldown': 'ldown_obs',
                'fcld': 'fcld_obs',
                'Wuh': 'wu_m3',
                'xsmd': 'xsmd',
                'lai': 'lai_obs',
                'kdiff': 'kdiff',
                'kdir': 'kdir',
                'wdir': 'wdir',
            }
        )
    # grid list determined by initial states
    list_grid = df_init.index

    # initialise dicts for holding results and model states
    dict_state = {}
    dict_output = {}

    # initial and final tsteps retrieved from forcing data
    tstep_init = df_forcing.index[0]
    tstep_final = df_forcing.index[-1]
    # tstep size retrieved from forcing data
    freq = df_forcing.index.freq

    # dict_state is used to save model states for later use
    dict_state = {
        # (t_start, grid): series_state_init.to_dict()
        (tstep_init, grid): pack_grid_dict(series_state_init)
        for grid, series_state_init
        in df_init.iterrows()
    }

    # remove 'problems.txt'
    if Path('problems.txt').exists():
        os.remove('problems.txt')

    if save_state:
        # use slower more functional single step wrapper

        # convert df to dict with `itertuples` for better performance
        dict_forcing = {row.Index: row._asdict()
                        for row in df_forcing.itertuples()}

        for tstep in df_forcing.index:
            # temporal loop
            # initialise output of tstep:
            # load met_forcing if the same across all grids:
            met_forcing_tstep = dict_forcing[tstep]
            # spatial loop
            for grid in list_grid:
                dict_state_start = dict_state[(tstep, grid)]
                # calculation at one step:
                # series_state_end, series_output_tstep = suews_cal_tstep_df(
                #     series_state_start, met_forcing_tstep)
                dict_state_end, dict_output_tstep = suews_cal_tstep(
                    dict_state_start, met_forcing_tstep)

                # update output & model state at tstep for the current grid
                dict_output.update({(tstep, grid): dict_output_tstep})
                dict_state.update({(tstep + 1*freq, grid): dict_state_end})

        # pack results as easier DataFrames
        df_output = pack_df_output(dict_output).swaplevel(0, 1)
        # drop unnecessary 'datetime' as it is already included in the index
        df_output = df_output.drop(columns=['datetime'], level=0)
        df_state_final = pack_df_state(dict_state).swaplevel(0, 1)

    else:
        # for multi-year run, reduce the whole df_forcing into {n_yr}-year chunks for less memory consumption
        grp_forcing_yr = df_forcing.groupby(df_forcing.index.year // n_yr)
        if len(grp_forcing_yr) > 1:
            df_state_init_yr = df_state_init.copy()
            list_df_output = []
            list_df_state = []
            for grp in grp_forcing_yr.groups:
                # get forcing of a specific year
                df_forcing_yr = grp_forcing_yr.get_group(grp)
                # run supy: actual execution done in the `else` clause below
                df_output_yr, df_state_final_yr = run_supy(
                    df_forcing_yr, df_state_init_yr)
                df_state_init_yr = df_state_final_yr.copy()
                # collect results
                list_df_output.append(df_output_yr)
                list_df_state.append(df_state_final_yr)

            # re-organise results of each year
            df_output = pd.concat(list_df_output).sort_index()
            df_state_final = pd.concat(
                list_df_state).sort_index().drop_duplicates()
            return df_output, df_state_final

        else:
            # for single-chunk run (1 chunk = {n_yr} years), directly put df_forcing into supy_driver for calculation
            # use higher level wrapper that calculate at a `block` level
            # for better performance

            # for grid in list_grid:
            #     dict_state_start_grid = dict_state[(tstep_init, grid)]
            #     dict_state_end, dict_output_array = suews_cal_tstep_multi(
            #         dict_state_start_grid,
            #         df_forcing)
            #     # update output & model state at tstep for the current grid
            #     dict_output.update({grid: dict_output_array})
            #     # model state for the next run
            #     dict_state.update({(tstep_final + freq, grid): dict_state_end})


            # # parallel run of grid_list for better efficiency
            # if os.name == 'nt':
            #     if __name__ == '__main__':
            #         p = Pool(min([len(list_grid), cpu_count()]))
            # else:
            #     p = Pool(min([len(list_grid), cpu_count()]))

            # # construct input list for `Pool.starmap`
            # construct input list for `dask.bag`
            list_input = [
                # (dict_state[(tstep_init, grid)], df_forcing)
                dict_state[(tstep_init, grid)]
                for grid in list_grid
            ]

            # on windows `processes` has issues when importing
            # so set `threads` here
            method_parallel = 'threads' if os.name == 'nt' else 'processes'
            list_res = db.from_sequence(list_input)\
                .map(suews_cal_tstep_multi, df_forcing)\
                .compute(scheduler=method_parallel)
            list_state_end, list_output_array = zip(*list_res)

            # collect output arrays
            dict_output = {
                grid: dict_output_array
                for grid, dict_output_array in zip(list_grid, list_output_array)
            }

            # collect final states
            dict_state_final_tstep = {
                (tstep_final + freq, grid): dict_state_end
                for grid, dict_state_end in zip(list_grid, list_state_end)
            }
            dict_state.update(dict_state_final_tstep)

            # save results as time-aware DataFrame
            df_output0 = pack_df_output_array(dict_output, df_forcing)
            df_output = df_output0.replace(-999., np.nan)
            df_state_final = pack_df_state(dict_state).swaplevel(0, 1)

    # drop ESTM for now as it is not supported yet
    # select only those supported output groups
    df_output = df_output.loc[:, ['SUEWS', 'snow', 'DailyState']]
    # trim multiindex based columns
    df_output.columns = df_output.columns.remove_unused_levels()

    # pack final model states into a proper dataframe
    df_state_final = pack_df_state_final(df_state_final, df_init)

    return df_output, df_state_final
Perform supy simulation.

    Parameters
    ----------
    df_forcing : pandas.DataFrame
        forcing data for all grids in `df_state_init`.
    df_state_init : pandas.DataFrame
        initial model states;
        or a collection of model states with multiple timestamps, whose last temporal record will be used as the initial model states.
    save_state : bool, optional
        flag for saving model states at each time step, which can be useful in diagnosing model runtime performance or performing a restart run.
        (the default is False, which instructs supy not to save runtime model states).
    n_yr : int, optional
        chunk size (`n_yr` years) to split simulation periods so memory usage can be reduced.
        (the default is 10, which implies 10-year forcing chunks used in simulations).

    Returns
    -------
    df_output, df_state_final : Tuple[pandas.DataFrame, pandas.DataFrame]
        - df_output: `output results <df_output_var>`
        - df_state_final: `final model states <df_state_var>`

    Examples
    --------

    >>> df_output, df_state_final = supy.run_supy(df_forcing, df_state_init)