import numpy as np

from veros import VerosSetup, veros_routine
from veros.variables import allocate, Variable
from veros.core.operators import numpy as npx, update, at

from veros.pyom_compat import load_pyom, setup_pyom

from test_base import compare_state


yt_start = -39.0
yt_end = 43
yu_start = -40.0
yu_end = 42


def set_parameter_pyom(pyom_obj):
    m = pyom_obj.main_module

    (m.nx, m.ny, m.nz) = (30, 42, 15)
    m.dt_mom = 4800
    m.dt_tracer = 86400 / 2.0
    m.runlen = 86400 * 365

    m.coord_degree = 1
    m.enable_cyclic_x = 1

    m.congr_epsilon = 1e-8
    m.congr_max_iterations = 10_000

    i = pyom_obj.isoneutral_module
    i.enable_neutral_diffusion = 1
    i.k_iso_0 = 1000.0
    i.k_iso_steep = 500.0
    i.iso_dslope = 0.005
    i.iso_slopec = 0.01
    i.enable_skew_diffusion = 1

    m.enable_hor_friction = 1
    m.a_h = 2.2e5
    m.enable_hor_friction_cos_scaling = 1
    m.hor_friction_cospower = 1

    m.enable_bottom_friction = 1
    m.r_bot = 1e-5

    m.enable_implicit_vert_friction = 1
    t = pyom_obj.tke_module
    t.enable_tke = 1
    t.c_k = 0.1
    t.c_eps = 0.7
    t.alpha_tke = 30.0
    t.mxl_min = 1e-8
    t.tke_mxl_choice = 2
    t.kappam_min = 2e-4

    i.k_gm_0 = 1000.0
    e = pyom_obj.eke_module
    e.enable_eke = 1
    e.eke_k_max = 1e4
    e.eke_c_k = 0.4
    e.eke_c_eps = 0.5
    e.eke_cross = 2.0
    e.eke_crhin = 1.0
    e.eke_lmin = 100.0
    e.enable_eke_superbee_advection = 1
    e.enable_eke_isopycnal_diffusion = 1

    i = pyom_obj.idemix_module
    i.enable_idemix = 1
    i.enable_idemix_hor_diffusion = 1
    i.enable_eke_diss_surfbot = 1
    i.eke_diss_surfbot_frac = 0.2
    i.enable_idemix_superbee_advection = 1
    i.tau_v = 86400.0
    i.jstar = 10.0
    i.mu0 = 4.0 / 3.0

    m.eq_of_state_type = 3


def set_grid_pyom(pyom_obj):
    m = pyom_obj.main_module
    ddz = [50.0, 70.0, 100.0, 140.0, 190.0, 240.0, 290.0, 340.0, 390.0, 440.0, 490.0, 540.0, 590.0, 640.0, 690.0]
    m.dxt[:] = 2.0
    m.dyt[:] = 2.0
    m.x_origin = 0.0
    m.y_origin = -40.0
    m.dzt[:] = ddz[::-1]
    m.dzt[:] *= 1 / 2.5


def set_coriolis_pyom(pyom_obj):
    m = pyom_obj.main_module
    m.coriolis_t[:, :] = 2 * m.omega * np.sin(m.yt[None, :] / 180.0 * np.pi)


def set_topography_pyom(pyom_obj):
    m = pyom_obj.main_module
    (X, Y) = np.meshgrid(m.xt, m.yt)
    X = X.transpose()
    Y = Y.transpose()
    m.kbot[...] = (X > 1.0) | (Y < -20)


def set_initial_conditions_pyom(pyom_obj):
    m = pyom_obj.main_module

    # initial conditions
    m.temp[:, :, :, :] = ((1 - m.zt[None, None, :] / m.zw[0]) * 15 * m.maskt)[..., None]
    m.salt[:, :, :, :] = 35.0 * m.maskt[..., None]

    # wind stress forcing
    taux = np.zeros(m.ny + 1)
    yt = m.yt[2 : m.ny + 3]
    taux = (0.1e-3 * np.sin(np.pi * (m.yu[2 : m.ny + 3] - yu_start) / (-20.0 - yt_start))) * (yt < -20) + (
        0.1e-3 * (1 - np.cos(2 * np.pi * (m.yu[2 : m.ny + 3] - 10.0) / (yu_end - 10.0)))
    ) * (yt > 10)
    m.surface_taux[:, 2 : m.ny + 3] = taux * m.masku[:, 2 : m.ny + 3, -1]

    t = pyom_obj.tke_module
    t.forc_tke_surface[2:-2, 2:-2] = (
        np.sqrt(
            (0.5 * (m.surface_taux[2:-2, 2:-2] + m.surface_taux[1:-3, 2:-2])) ** 2
            + (0.5 * (m.surface_tauy[2:-2, 2:-2] + m.surface_tauy[2:-2, 1:-3])) ** 2
        )
        ** 1.5
    )


def set_forcing_pyom(pyom_obj):
    m = pyom_obj.main_module
    t_star = (
        15 * np.invert((m.yt < -20) | (m.yt > 20))
        + 15 * (m.yt - yt_start) / (-20 - yt_start) * (m.yt < -20)
        + 15 * (1 - (m.yt - 20) / (yt_end - 20)) * (m.yt > 20.0)
    )
    t_rest = m.dzt[None, -1] / (30.0 * 86400.0) * m.maskt[:, :, -1]
    m.forc_temp_surface = t_rest * (t_star - m.temp[:, :, -1, m.tau - 1])


class ACCSetup(VerosSetup):
    @veros_routine
    def set_parameter(self, state):
        settings = state.settings
        settings.identifier = "acc"

        settings.nx, settings.ny, settings.nz = 30, 42, 15
        settings.dt_mom = 4800
        settings.dt_tracer = 86400 / 2.0
        settings.runlen = 86400 * 365

        settings.x_origin = 0.0
        settings.y_origin = -40.0

        settings.coord_degree = True
        settings.enable_cyclic_x = True

        settings.enable_neutral_diffusion = True
        settings.K_iso_0 = 1000.0
        settings.K_iso_steep = 500.0
        settings.iso_dslope = 0.005
        settings.iso_slopec = 0.01
        settings.enable_skew_diffusion = True

        settings.enable_hor_friction = True
        settings.A_h = 2.2e5
        settings.enable_hor_friction_cos_scaling = True
        settings.hor_friction_cosPower = 1

        settings.enable_bottom_friction = True
        settings.r_bot = 1e-5

        settings.enable_implicit_vert_friction = True

        settings.enable_tke = True
        settings.c_k = 0.1
        settings.c_eps = 0.7
        settings.alpha_tke = 30.0
        settings.mxl_min = 1e-8
        settings.tke_mxl_choice = 2
        settings.kappaM_min = 2e-4

        settings.K_gm_0 = 1000.0
        settings.enable_eke = True
        settings.eke_k_max = 1e4
        settings.eke_c_k = 0.4
        settings.eke_c_eps = 0.5
        settings.eke_cross = 2.0
        settings.eke_crhin = 1.0
        settings.eke_lmin = 100.0
        settings.enable_eke_superbee_advection = True
        settings.enable_eke_isopycnal_diffusion = True

        settings.enable_idemix = 1
        settings.enable_idemix_hor_diffusion = 1
        settings.enable_eke_diss_surfbot = 1
        settings.eke_diss_surfbot_frac = 0.2
        settings.enable_idemix_superbee_advection = 1
        settings.tau_v = 86400.0
        settings.jstar = 10.0
        settings.mu0 = 4.0 / 3.0

        settings.eq_of_state_type = 3

        var_meta = state.var_meta
        var_meta.update(
            t_star=Variable("t_star", ("yt",), "deg C", "Reference surface temperature"),