def METHOD_NAME(solution, sol_old):
    # source term (external pressure gradient
    f = inner(source, test)*dx
    # vertical diffusion (integrated by parts)
    f += -viscosity_v*inner(Dx(solution, 2), Dx(test, 2)) * dx
    # interface term
    diff_flux = viscosity_v*Dx(solution, 2)
    f += (dot(avg(diff_flux), test('+'))*normal[2]('+')
          + dot(avg(diff_flux), test('-'))*normal[2]('-')) * dS_h
    # symmetric interior penalty stabilization
    l = Constant(depth/layers)
    nb_neigh = 2
    o = 1
    d = 3
    sigma = Constant((o + 1)*(o + d)/d * nb_neigh / 2) / l
    gamma = sigma*avg(viscosity_v)
    f += gamma * dot(jump(solution), test('+')*normal[2]('+') + test('-')*normal[2]('-')) * dS_h
    # boundary term
    uv_bot_old = sol_old + Dx(sol_old, 2)*l*0.5
    uv_bot = solution + Dx(solution, 2)*l*0.5  # solver fails
    uv_mag = sqrt(uv_bot_old[0]**2 + uv_bot_old[1]**2) + Constant(1e-12)
    bnd_flux = bottom_drag*uv_mag*uv_bot
    ds_bottom = ds_t
    f += dot(bnd_flux, test)*normal[2] * ds_bottom
    return f