This is a very nice question, but to answer that, we should first clarify the OP.
-- First scenario: Abrupt change
The OP considers a scenario in which a ferromagnetic material prepared in an ordered state (eg., with $ \langle m \rangle = \uparrow $) is suddenly subjected to a finite “opposite” external magnetic field (ie., antiparallel to the initial magnetisation) while the temperature is kept fixed below the critical temperature. The question is about the latent heat involved in this purported “phase transition of 1st order”.
Firstly, I think this scenario cannot be called a phase transition in the common sense. Such a procedure will induce an abrupt change in the properties of the system (esp. its ground-state). Free energy will be discontinuous itself, so it is a “zeroth-order/discontinuous phase transition”, if you wish. Phase transitions are usually defined as processes in which a slow variation of a thermodynamic variable (usually, temperature) leads to drastic changes in the thermodynamic phase of a system. Abrupt changes are trivial in this regard: we know a priori that they induce qualitative changes in the system. In our case, if one changes the external magnetic field abruptly, and waits long enough to let the system relax to its equilibrium, one will see that the magnetisation of the system is reversed to be parallel to the final external field; the system adapts itself to the applied field if we allow it to relax (see diagram below).
*************
apply H ** complicated **
initial phase ———————————> ** relaxation ** ———> final phase
** process **
*************
Notice that the devices of equilibrium statistical mechanics cannot deal with the intervening inherently non-equilibrium relaxation process. They can only explain the equilibrium intial and final (ordered) phases. In our scenario, we implicitly assume that we wait long enough for the system to relax.
So, the scenario explained by OP is an abrupt change in the system, not a continuous or discontinuous phase transition — in its proper sense.
-- Mean-field analysis
To see the above explanations more quantitatively, let's model a ferromagnetic material with a nearest-neighbour Ising model ($ S = \frac{1}{2} $), and “solve” it within the mean-field approximation. The derivations of the mean-field Hamiltonian and partition function are given in modern statistical physics textbooks and will not be repeated here (see eg., Schwabl, F. “Statistical Mechanics” (2010) [WCat]). The dimensionless mean-field free energy density reads
$$ f(h, T; m) := \frac{1}{T_c} \frac{F}{N} = \frac{1}{2} m^2 - \frac{T}{T_c} \ln \big( 2\cosh( \mathcal{M}_h ) \big) $$
where $ \mathcal{M}_h := \frac{T_c}{T} ( h + m ) $, $ h $ is the rescaled external magnetic field, $ h := h_{ext}/T_c $, $ m $ is the mean-field, $ T_c $ is the critical temperature for spontaneous magnetisation, and $ T $ is the temperature. Note that we use natural units where $ k_B = 1 = \hbar $.
To determine the mean-field $ m $, we minimize the free energy and obtain the self-consistent mean-field equation,
$$ m = \tanh( \mathcal{M}_h ) ~. $$
The nature of the solutions for $ T < T_c $ differs from those for $ T > T_c $: ie., in absence of external fields, $ h \rightarrow 0 $, we have a finite spontaneous magnetisation ($ m \neq 0 $) when $ T < T_c $, but no spontaneous magnetisation ($ m = 0 $) when $ T > T_c $.
From the free energy, one can readily obtain the entropy*,
$$ s := \frac{S}{N} = - \frac{1}{N} \frac{\partial F}{\partial T} = m \, \mathcal{M}_h + \ln \big( 2 \cosh(\mathcal{M}_h) \big) ~, $$
and specific heat at constant volume, $ c_V $,
$$ c_V := \frac{C_V}{N} = \frac{1}{N} T \frac{\partial S}{\partial T} \Big\vert_V = -\frac{1}{N} T \frac{\partial^2 F}{\partial T^2} \Big\vert_V = \mathcal{M}_h \big( \frac{3 m + h}{1 - \frac{t}{1 - m^2}} - 2 m \big) ~. $$
-- THE FOLLOWING PART IS NOT IN FINAL FORM YET
From these relations, we can obtain the difference between the thermodynamic quantities in the scenario above: we prepare the system in a preferred magnetised state (say, $ \uparrow $) with a tiny external field $ h_\uparrow $ in a temperature $ T < T_c $, and then suddenly apply a finite external magnetic field $ h_\downarrow $ in the opposite direction. Then we calculate the free energy and entropy in the initial and final phases and changes thereof:
$$ \Delta f := f(h_\downarrow, T) - f(h_\uparrow, T) \\ \Delta s := s(h_\downarrow, T) - s(h_\uparrow, T) $$
from this we obtain the exchanged heat in the process,
$$ \Delta q = T \Delta s ~. $$
For concreteness, let's suppose the following
$$ \begin{align*} \frac{T}{T_c} &= \frac{1}{2} \quad \text{: deep in the ordered phase} \\ h_\uparrow &= +10^{-3} \quad \text{: small positive field for initial ordering} \\ h_\downarrow &= -1 \quad \text{: large opposite field applied suddenly} \end{align*} $$ where all the quantities are dimensionless. Then, we obtain**
$$ \begin{align*} \Delta f &= -0.989 \\ \Delta s &= 4.220 \\ \Delta q &= 8.440 ~. \end{align*} $$
Notice the finite jump in free energy. This process is actually similar to magnetic cooling.
One can also consider a variation of the original scenario where, at a fixed $ T < T_c $, the strength of the magnetic field is slowly (quasi-statically) increased from a positive value to some negative value. We observe that in this case, $ m $ has a jump at $ h = 0 $, therefore, a discontinuous (1st-order) phase transition has happened.
We can easily see that in this case, entropy is a smooth function of $ h $. This smoothness can be demonstrated analytically as follows:
Near the transition point, $ \frac{T}{T_c} \rightarrow 1 $, and for small $ m $ and $ h $, we can expand the free energy as
$$ f \approx \frac{1}{2} r \, m^2 + \frac{1}{12} m^4 - h \, m ~, $$
where $ r = \frac{T - T_c}{T_c} $, and only the leading linear term in $ h $ is kept (“linear response” approximation). The mean-field equation will then be
$$ \frac{\partial f}{\partial m} = r \, m + \frac{1}{3} m^3 - h = 0 \Rightarrow m \sim \frac{h}{r} ~, $$
within the linear response approximation. Hence, entropy reads
$$ s = - \frac{\partial f}{\partial T} \sim \frac{1}{2} m + m^2 - \frac{1}{3 h} m^5 ~. $$
If we replace $ m $ in the above equation by its approximate value, $m \sim \frac{h}{r} $, we obtain an expansion of the entropy in terms of the applied field $ h $,
$$ s \sim \frac{h}{r} - (\frac{h}{r})^2 + \frac{h^2}{r^3} ~, $$
which shows that entropy is an analytic (smooth) function of $ h $ near the transition point.
${}^{\ast}$ We have used the mean-field equation, and $ \frac{\partial m}{\partial T} = \frac{q}{ 1 - \frac{T}{1 - m^2} } $ derived from it.
${}^{\ast\ast}$ The Python code to compute and visualize all of these is accessible here.