function sea_ice()
% solves sea ice model
%   dH/dt = - Fo - Q/(1+H)  Q<=0
%           - Fo - Q        Q>0
%   for H>=0, OR H = 0
% uses a single ODE for both cases H>=0 and H<0.  Numerical errors
% causes H to go slightly below 0; dHdt is then set to zero until it
% becomes positive allowing ice to start to grow again

%% Parameters
Fo = 0;
Q0 = 0;
omega = 2*pi*.1;
Q = @(t) Q0 + cos(omega*t);

%% Solve ODE
opts = odeset('reltol',1e-8);
[t,H] = ode45(@dHdt_fun,[0 2*2*pi/omega],0,opts);
% ODE
function dHdt = dHdt_fun(t,H)
    dHdt = ( -Fo-Q(t)./(1+H) ).*( Q(t)<=0 ).*( H>=0 ) + ...
            ( -Fo-Q(t) ).*( Q(t)>0 ).*( H>=0 ) + ...
            max(0,-Fo-Q(t)./(1+H)).*( H<0 );
end

%% Plot solution
figure(1); clf; width = 2; fontsize = 14;
set(gcf,'DefaultAxesFontSize',fontsize,'DefaultTextFontSize',fontsize);
set(gcf,'units','centimeters','Paperpositionmode','auto','position',[2 2 20 10]);
    plot(omega/2/pi*t,H,'linewidth',width);
%     hold on; plot(omega/2/pi*t,dHdt_fun(t,H)); % derivative
    hold on; plot(omega/2/pi*t,max(0,-1-Q(t)/Fo),'r--'); % quasi-steady approximation (small omega)
    xlabel('$\omega t / 2\pi$','interpreter','latex');
    ylabel('$H / [H]$','interpreter','latex');

% %% Save figure  
% text(0.02,0.95,['$\omega/2\pi = $',num2str(omega/2/pi)],'interpreter','latex','units','normalized','horizontalalignment','left','verticalalignment','top');
% print(gcf,'-depsc2',[mfilename,'_omega',num2str(omega/2/pi),'.eps'],'-loose');

end