def METHOD_NAME(solution):
    # vertical diffusion operator integrated by parts
    # f = -diffusivity_v*inner(Dx(solution, 2), Dx(test, 2)) * dx
    # interface term
    # diffFlux = diffusivity_v*Dx(solution, 2)
    # f += (dot(avg(diffFlux), test('+'))*normal[2]('+') +
    #      dot(avg(diffFlux), test('-'))*normal[2]('-')) * dS_h
    # symmetric interior penalty stabilization
    # L = Constant(depth/layers)
    # nbNeigh = 2
    # d = 3
    # sigma = Constant((deg + 1)*(deg + d)/d * nbNeigh / 2) / L
    # gamma = sigma*avg(diffusivity_v)
    # jump_test = test('+')*normal[2]('+') + test('-')*normal[2]('-')
    # f += gamma * dot(jump(solution), jump_test) * dS_h
    def grad_v(a):
        return as_vector((0, 0, Dx(a, 2)))
    n = as_vector((0, 0, normal[2]))
    grad_test = grad_v(test)
    diff_flux = diffusivity_v*grad_v(solution)
    diff_flux_jump = diffusivity_v*jump(solution, n)
    f = -inner(grad_test, diff_flux)*dx
    # Interior penalty as in horizontal case
    degree_h, degree_v = fs.ufl_element().degree()
    dim = 3.0
    sigma = (degree_v + 1.0)*(degree_v + dim)/dim/v_elem_size
    # sigma = 1.0/v_elem_size
    alpha = avg(sigma)
    ds_interior = (dS_h)
    f += -alpha*inner(jump(test, n), diff_flux_jump)*ds_interior
    f += +inner(avg(grad_test), diff_flux_jump)*ds_interior
    f += +inner(jump(test, n), avg(diff_flux))*ds_interior
    # symmetric boundary terms
    f += inner(test, dot(diff_flux, n))*(ds_t + ds_b)
    return f