%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ice_dynamics_web_version1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

set(0,'defaulttextfontsize',16); set(0,'defaultaxesfontsize',16);

flag_pickup=0;

if flag_pickup
  load('pickup_final.mat');
  p.t_final=15.02e6*p.yr;
  flag=1;
else
  %% set parameters:
  %% ---------------
  p=set_parameters;

  
  %% calculate insolation and surface temperature:
  %% ---------------------------------------------
  [p,W]=calc_daily_insolation(p);
  Wsols=calc_solstice_insolation(p);
  %% calculate the surface temperature assuming that there is no
  %% internal heating 
  Ts=calc_Ts_new(p,zeros(size(p.phi)),W);
%% calculate the initial freezing temperature
  Tf=(0.0901-0.0575*p.mSb)-7.61e-4*(p.g*p.h*p.rhoI/1e4)+273.16;
  Tf=Tf*ones(size(p.theta));

  dh_dtheta=zeros(size(p.theta)); %  first derivative of h
  dh_dtheta2=zeros(size(p.theta)); %  2nd derivative of h
  dq_dtheta=zeros(size(p.theta)); %  first derivative of q
  if p.flag_ocean
    To=Tf; %  first derivative of h
    dTo_dtheta=zeros(size(p.theta)); %  1nd derivative of h
    dTo_dtheta2=zeros(size(p.theta)); %  2nd derivative of h
    if p.flag_salt
      Sb=p.mSb*ones(size(p.theta));
      dSb_dtheta=zeros(size(p.theta)); %  first derivative of Sb
      dSb_dtheta2=zeros(size(p.theta)); %  2nd derivative of Sb
    end
  end

  xi=0:p.dxi:1; % xi is (z-zb)/(zs-zb)
  p.n_xi=length(xi);

  p.n_theta=length(Ts);
  v=zeros(p.n_theta,p.n_xi);
% initialize dv_dz with some typical scale:
  dv_dz=zeros(p.n_theta,p.n_xi)+5e-25;
  dv_dz_p=zeros(p.n_theta,p.n_xi);
  dv_dz_m=zeros(p.n_theta,p.n_xi);
  Qc=p.Qm*ones(p.n_theta,p.n_xi);
  
% term that is used in the source term
  A2=-p.kappaI*p.mu*p.CpI/p.L*(Ts-Tf);
  A3=-p.Q_g/p.L/p.rho0;
  
  A4=cot(p.theta)/p.Radius;
  
% start with uniform ice depth
  hh=ones(size(p.theta))*p.h;
  h_initial=hh;
  
  t_old=0;
  h_old=hh;
  
  flag=1;

  p.flag_first_call_to_calc_T_ice_convection=1;
  sin_theta=sin(p.theta);
  sum_sin_theta=sum(sin(p.theta));
  ice_bottom_heat_flux=zeros(p.n_theta,1);
  kappa_conv_2d=zeros(p.n_theta,p.n_xi-1);
  tidal_heating=zeros(p.n_theta,p.n_xi);
  z_convective_conductive_transition=zeros(p.n_theta,1);
  total_tidal_heating_ice=zeros(p.n_theta,1);
  h_conducting=zeros(p.n_theta,1);
% initialize ice temperature with a linear profile in the vertical:
  T_ice=zeros(p.n_theta,p.n_xi);
  for i=1:p.n_theta
    T_ice(i,:)=Tf(i)+(Ts(i)-Tf(i))*xi;
  end
  t=0;
end

mHconvc=0;

%% integrate in time to steady state:
while t<(p.t_final+p.dt)
  %% fprintf(1,'t=%g kyr\n',t/(1000*p.yr));
  %%fprintf(1,'max hh=%g, min hh=%g, mean hh=%g, time=%g\n',max(hh),min(hh),mean(hh),t);
  dh_dtheta(2:end-1)=0.5*(hh(3:end)-hh(1:end-2))/p.dphi;  % center difference 
  %% dh_dtheta is zero exactly at the boundary which is located at i=0.5. 
  %% For that h_0=h_1. Thus dh_dtheta_1=(h_2-h_0)/2/dtheta=(h_2-h_1)/2/dtheta
  dh_dtheta(1)=0.5*(hh(2)-hh(1))/p.dphi;
  %% For that h_end=h_{end+1}. 
  %% Thus dh_dtheta_end=(h_{end+1}-h_{end-1})/2/dtheta=(h_end-h_{end-1})/2/dtheta
  dh_dtheta(end)=0.5*(hh(end)-hh(end-1))/p.dphi;

  tmp=dh_dtheta.*sin_theta;
  dh_dtheta2(2:end-1)=0.5*(tmp(3:end)-tmp(1:end-2))/p.dphi;  % center diff
  %% tmp is zero exactly at the boundary which is located at i=0.5. 
  %% For that tmp_0=-tmp_1. Thus dh_dtheta2_1=(tmp_2-tmp_0)/2/dtheta=(tmp_2+tmp_1)/2/dtheta
  dh_dtheta2(1)=0.5*(tmp(2)+tmp(1))/p.dphi;
  %% tmp is zero exactly at the boundary which is located at i=end+0.5. 
  %% For that tmp_{end+1}=-tmp_end. Thus dh_dtheta2_end=(tmp_{end+1}-tmp_{end-1})/2/dtheta=-(tmp_end+tmp_{end-1})/2/dtheta
  dh_dtheta2(end)=-0.5*(tmp(end)+tmp(end-1))/p.dphi;
  dh_dtheta2=dh_dtheta2./sin_theta;
  
  if p.flag_salt
    Tf=(0.0901-0.0575*Sb)-7.61e-4*(p.g*hh*p.rhoI/1e4)+273.16; % (-3.7019 C)  
  else
    Tf=(0.0901-0.0575*p.Sb)-7.61e-4*(p.g*hh*p.rhoI/1e4)+273.16; % (-3.7019 C)  
  end
  %%  tmp1=(hh.*abs(dh_dtheta)).^p.n;
  tmp1=(hh.*dh_dtheta).^p.n;
  
  eta_eff=effective_viscosity(p,T_ice,dv_dz,Tf,Qc);
  for i=1:p.n_theta

    %% calculate temperature as function of depth:
    %% -- diffusion only:
    %%[TT]=calc_T_ice_diffusion_only(Ts(i),Tf(i),xi);
    %% -- convection and diffusion:
    %% number of interior temperature points that need to be solved for:
    N=length(xi)-2; 
    %% location of temperature grid points, including top and bottom:
    z=[(xi-p.mu)*hh(i)];
    dz=z(2)-z(1);
    %% kappa points are between temperature points, and there is
    %% one more such point than temperature points:
    z_kappa=[z(2:end-1)-dz/2,z(end-1)+dz/2];
    if 1
      %% calculate the effective viscosity 
      eta=eta_eff(i,:);
      if p.tidal_heating 
        tidal_heating(i,:)=2*p.H_max(i)./(eta./p.eta_max+ ...
                                          p.eta_max./eta);
      else
        tidal_heating(i,:)=0;
      end
      total_tidal_heating_ice(i)=sum(tidal_heating(i,:))*dz;
      if t==0
        heating_bottom=p.Q_g;
      else
        heating_bottom=ice_bottom_heat_flux(i);
      end
      %% depth of ice bottom:
      z_b=-p.mu*hh(i);
      %% depth of ice surface:
      z_s=(1-p.mu)*hh(i);
      %% thickness of conducting layer at top of ice, below which
      %% convection occurs: 
      z_convective_conductive_transition(i)= ...
       interp1(T_ice(i,:),z,Tf(i)-p.delta_T_over_convective_depths_Hussmann,'spline');
      if z_convective_conductive_transition(i)>z_s; 
        z_convective_conductive_transition(i)=z_s; 
      elseif z_convective_conductive_transition(i)<z_b;
        z_convective_conductive_transition(i)=z_b;
      end
      h_conducting(i)=z_s-z_convective_conductive_transition(i);

      kappa_diff_conv= p.kappaI+(p.kappaI_convection-p.kappaI) ...
          .*heavyside_smooth((z_convective_conductive_transition(i)-z_kappa)/p.dz_conv) ...
          .*heavyside_smooth((z_kappa-(z_b+2*p.dz_conv))/(p.dz_conv)); ...
      kappa_conv_2d(i,:)=kappa_diff_conv;
    else
      kappa_diff_conv=p.kappaI+[z*0 0];
      tidal_heating=z*0; 
    end

    %% build matrix equation to solve diffusion equation in ice:
    Diff_matrix=zeros(p.n_xi-2,p.n_xi-2);
    [ice_bottom_heat_flux(i),Temperature_profile,Diff_matrix] ...
      =calc_T_ice_convection(p,N,Tf(i),Ts(i),kappa_diff_conv,tidal_heating(i,:),z,Diff_matrix);
    T_ice(i,2:end-1)=Temperature_profile';
    %% update the top and bottom temperatures which are not
    %% calculated by the diffusion equation:
    T_ice(i,1)=Tf(i); T_ice(i,end)=Ts(i);
    TT=squeeze(T_ice(i,:));

    eta_eff(i,:)=p.eta0*exp(p.Ql/(p.R_gas_const*Tf(i))*(Tf(i)./TT-1));
    dv_dz(i,:)=2*(1-p.mu)*p.rhoI*p.g/p.Radius./eta_eff(i,:).*xi ...
        *hh(i)*dh_dtheta(i);
    v(i,:)=-hh(i)*(sum(dv_dz(i,:),2)-cumsum(dv_dz(i,:),2))*p.dxi;
    v(i,TT<p.T_transition_brittle)=0;
    q(i)=hh(i)*sum(v(i,:),2)*p.dxi;
  
  end
  
  %%calculate the ice flux times sin(p.theta)
  tmp=q.*sin_theta;

  dq_dtheta(2:end-1)=0.5*(tmp(3:end)-tmp(1:end-2))/p.dphi;  % center diff
  
  %% q is zero exactly at the boundary which is located at i=0.5. 
  %% For that q_0=-q_1. Thus dq_dtheta_1=(q_2-q_0)/2/dtheta=(q_2+q_1)/2/dtheta
  dq_dtheta(1)=0.5*(tmp(2)+tmp(1))/p.dphi;
  %% For that q_{end+1}=-q_end. 
  %% Thus dq_dtheta_1=(q_{end+1}-q_{end-1})/2/dtheta=-(q_end+q_{end-1})/2/dtheta
  dq_dtheta(end)=-0.5*(tmp(end)+tmp(end-1))/p.dphi;
  
  if p.flag_ocean
    %% local ocean depth 
    mh=sum(hh.*sin_theta)/sum_sin_theta;
    ho=p.Ho-p.mu*(hh-mh);

    %% first derivative of ocean temperature  
    dTo_dtheta(2:end-1)=0.5*(To(3:end)-To(1:end-2))/p.dphi;  % center difference 
    dTo_dtheta(1)=0.5*(To(2)-To(1))/p.dphi;
    dTo_dtheta(end)=0.5*(To(end)-To(end-1))/p.dphi;

    %% mutiply by sin(theta)    
    dTo_dtheta=ho.*dTo_dtheta.*sin_theta;
    
    %% second derivative of ocean temperature  
    dTo_dtheta2(2:end-1)=0.5*(dTo_dtheta(3:end)-dTo_dtheta(1:end-2))/p.dphi; 
    dTo_dtheta2(1)=0.5*(dTo_dtheta(2)+dTo_dtheta(1))/p.dphi;
    dTo_dtheta2(end)=-0.5*(dTo_dtheta(end)+dTo_dtheta(end-1))/p.dphi;

    %% devide by sin(theta)    
    dTo_dtheta2=dTo_dtheta2./sin_theta;
    
    if p.flag_salt
      %% first derivative of ocean salinity
      dSb_dtheta(2:end-1)=0.5*(Sb(3:end)-Sb(1:end-2))/p.dphi;  % center difference 
      dSb_dtheta(1)=0.5*(Sb(2)-Sb(1))/p.dphi;
      dSb_dtheta(end)=0.5*(Sb(end)-Sb(end-1))/p.dphi;

      %% mutiply by sin(theta)    
      dSb_dtheta=ho.*dSb_dtheta.*sin_theta;
      
      %% second derivative of ocean salinity
      dSb_dtheta2(2:end-1)=0.5*(dSb_dtheta(3:end)-dSb_dtheta(1:end-2))/p.dphi; 
      dSb_dtheta2(1)=0.5*(dSb_dtheta(2)+dSb_dtheta(1))/p.dphi;
      dSb_dtheta2(end)=-0.5*(dSb_dtheta(end)+dSb_dtheta(end-1))/p.dphi;

      %% devide by sin(theta)    
    dSb_dtheta2=dSb_dtheta2./sin_theta;
        
    end

    %% calculate the source term
    S=(ice_bottom_heat_flux'-p.beta*(To-Tf))/(p.rho0*p.L);  % in m/s
  else
    %% calculate the source term: works for no-ocean case:
    S=(ice_bottom_heat_flux'-p.Q_g)/(p.rho0*p.L);
  end


  %% calculate the numerical diffusion term  
  D=0;
  
  dhdt=-dq_dtheta./sin_theta/p.Radius+S/p.mu+D;

  if t>t_old
    t_old=t_old+p.dt_print;
    
    if flag
      figure;
      axis([t/p.yr,p.t_final/p.yr,0.75*p.h/1e3,1.25*p.h/1e3]);
      ylabel('h (km)');
      xlabel('Time (yr)');
      hold on;
      flag=0;
      plot(t/p.yr,min(hh)/1e3,'.b',t/p.yr,max(hh)/1e3,'.r');
      legend('min h','max h')
    else
      plot(t/p.yr,min(hh)/1e3,'.b',t/p.yr,max(hh)/1e3,'.r');
      save('pickup.mat');
    end
    drawnow;
    fprintf(1,'t=%g yrs; <h>=%g m; <h_conv>=%g m; max(abs(dh/dt))=%g\n' ...
            ,t/p.yr,sum(hh.*sin_theta)/sum_sin_theta,mHconvc,max(abs(dhdt)));
    h_old=hh;
    
  end

  hh=hh+p.dt*dhdt;
  if p.flag_ocean
    To=To+p.dt*((p.Q_g-p.beta*(To-Tf))/p.rho0/p.CpO+p.kappa_GM*dTo_dtheta2)./ho; 
    if p.flag_salt
      mSb=sum(Sb.*sin_theta)/sum_sin_theta;
      Sb=Sb+p.dt*(mSb.*S+p.kappa_GM*dSb_dtheta2)./ho; 
    end
  end

  %% update kappa_convective
  if t>p.dt
    mHconvc=sum(sin_theta.*(hh-h_conducting'))/sum_sin_theta;
    %% calculate the Rayleigh number
    Ra=p.const_Rayleigh*mHconvc^3;
    if Ra>p.Rac
       if p.simple_convection
	 p.kappaI_convection=p.kappaI*p.kappaI_conv_cond_factor;
       else
	 p.kappaI_convection=p.kappaI;
       end
    else
      p.kappaI_convection=p.kappaI;
    end
    %% update the surface temperature based on the updated internal
    %% heating 
    Ts=calc_Ts_new(p,(total_tidal_heating_ice+ice_bottom_heat_flux)',W);
  end

  
  if ~isfinite(mean(hh))
    figure;
    plot(p.phi/pi*180,h_old/1000);
    title(['Ice depth at t=' sprintf('%g',(t_old-p.dt_print)/p.yr/1000) ' kyr' ] );
    xlabel('Latitude');
    ylabel('Depth (km)');

    fprintf(1,'*** error: found infinite hh at t=%g yr.\n',t/p.yr);
    return
  end
  
  t=t+p.dt;
end

%% plot results:
plot_results(p,v,hh,h_initial,q,S,xi,dv_dz,dv_dz_p,dv_dz_m,Tf,Ts,T_ice,A2,A3,Qc);
debugging_plots(p,t,xi,W,Wsols,h_conducting,total_tidal_heating_ice, ...
                ice_bottom_heat_flux,hh,dh_dtheta,dh_dtheta2,Ts,Tf, ...
                z_convective_conductive_transition,T_ice,eta_eff, ...
                tidal_heating,kappa_conv_2d);

if p.flag_ocean;
  plot_ocean(p,To,Tf);
  if p.flag_salt;
    plot_salt(p,Sb);
  end
end

Ts_sols=calc_Ts(p,(total_tidal_heating_ice+ice_bottom_heat_flux)',Wsols);
Ts_sols_no_tidal_heating=calc_Ts(p,ice_bottom_heat_flux',Wsols);

%% save results to file:
% Thickness of stagnant ice. 
h_s=-A2/A3;  
% Temperature within the ice as a function of latitude and depth:
TT1=((Ts-Tf)'*xi+Tf'*ones(1,p.n_xi));

if p.flag_ocean
  if p.flag_salt
    save(['output/p.eta0_' sprintf('%g',p.eta0) ...
          '_with_ocean_with_salt_YA.mat'], ... 
         'xi','hh','h_s','Ts','Tf','To','Sb','q','S','v','dv_dz','p', ... 
         'h_conducting','total_tidal_heating_ice', ...
         'ice_bottom_heat_flux','z_convective_conductive_transition', ...
         'T_ice','eta_eff','tidal_heating','kappa_conv_2d'); 
  else
    save(['output/p.eta0_' sprintf('%g',p.eta0) ...
          '_with_ocean_YA.mat'], ... 
         'xi','hh','h_s','Ts','Tf','To','q','S','v','dv_dz','p', ... 
         'h_conducting','total_tidal_heating_ice', ...
         'ice_bottom_heat_flux','z_convective_conductive_transition', ...
         'T_ice','eta_eff','tidal_heating','kappa_conv_2d'); 
  end
else
  save(['output/p.eta0_' sprintf('%g',p.eta0) ...
        '_no_ocean_YA.mat'], ... 
       'xi','hh','h_s','Ts','Tf','q','S','v','dv_dz','p', ... 
       'h_conducting','total_tidal_heating_ice', ...
       'ice_bottom_heat_flux','z_convective_conductive_transition', ...
       'T_ice','eta_eff','tidal_heating','kappa_conv_2d'); 
end
save('pickup_final.mat');

% contourf(p.phi'/pi*180*ones(1,p.n_xi),hh'*xi-p.mu*hh'*ones(1,p.n_xi),((Ts-Tf)'*xi+Tf'*ones(1,p.n_xi)));colorbar


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [W]=calc_insolation(p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W=p.S0/4/pi/(1-p.e^2)^2*(4*cos(p.phi)+p.epsilon^2*sin(p.phi).*tan(p.phi) ... 
    -cos(p.phi)*sin(p.epsilon)^2);
i1=find((p.phi>(pi/2-p.epsilon))|(p.phi<(-pi/2+p.epsilon)));
W(i1)=p.S0/pi/(1-p.e^2)^2* ... 
    (sin(p.epsilon)+0.25/sin(p.epsilon)*cos(p.phi(i1)).^2); %(pi/2-p.phi(i1)).^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p,W]=calc_daily_insolation(p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
omega=0; % precession of Jupiter. Some default number. The
            % observed value is probably different
delta=asin(sin(p.epsilon)*sin(p.lambda));  % declination angle 
p.dt_dlambda=1./(1+p.e*cos(p.lambda-omega-pi)).^2/2/pi; % normalized
                                                        % dt_dlambda
                                                        % (i.e.,
                                                        % dt_dlambda/T) 
tmp1=p.S0/2/pi^2/(1-p.e^2)^2;
W=zeros(length(p.phi),length(p.lambda));
for i=1:length(p.phi)
  H0=acos(-tan(p.phi(i))*tan(delta));  % hour angle at sunrise and sunset 
  i0=find(imag(H0)==0);
  W(i,i0)=2*pi*tmp1*(1+p.e*cos(p.lambda(i0)-omega-pi)).^2.* ... 
         (H0(i0).*sin(p.phi(i)).*sin(delta(i0))+ cos(p.phi(i)).* ...
          cos(delta(i0)).*sin(H0(i0))); 
  i0=find(p.phi(i)>(pi/2-delta)|p.phi(i)<(-pi/2-delta)); %polar day
  W(i,i0)=p.S0*sin(p.phi(i))*sin(delta(i0))/(1-p.e^2)^2.* ... 
          (1+p.e*cos(p.lambda(i0)-omega-pi)).^2;  
  i0=find(p.phi(i)>(pi/2+delta)|p.phi(i)<(-pi/2+delta)); %polar night
  W(i,i0)=0;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [W]=calc_solstice_insolation(p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H0=acos(-tan(p.phi)*tan(p.epsilon));  % hour angle at sunrise and sunset 
W=p.S0/pi*(H0.*sin(p.phi)*sin(p.epsilon)+ ... 
            cos(p.phi).*sin(H0)*cos(p.epsilon));
i0=find(p.phi>(pi/2-p.epsilon));
W(i0)=p.S0*sin(p.phi(i0))*sin(p.epsilon);
% polar winter 
i0=find(p.phi<(-pi/2+p.epsilon));
W(i0)=0;

%% for surface temperature claculation see: http://arxiv.org/abs/1608.07372
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Ts]=calc_Ts(p,Q,W)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ts=((W*(1-p.alpha)*(1-p.p1)+Q)/p.sigma/p.emissivity).^0.25;

%% for surface temperature claculation see: http://arxiv.org/abs/1608.07372
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Ts]=calc_Ts_new(p,Q,W)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(p.phi)
  tmp3=p.dt_dlambda.*((Q(i)+W(i,:)*(1-p.alpha)*(1-p.p1))/p.sigma/ ...
                      p.emissivity).^0.25; 
  Ts(i)=sum(tmp3)*p.dl;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_results(p,v,hh,h_initial,q,S,xi,dv_dz,dv_dz_p,dv_dz_m,Tf,Ts,T_ice,A2,A3,Qc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
plot(p.phi/pi*180,hh/1000,'r')
hold on
plot(p.phi/pi*180,h_initial/1000,'b');
title(['Ice thickness, n=' num2str(p.n) '  \eta_0=' sprintf('%g',p.eta0)] );
xlabel('Latitude');
ylabel('Thickness (km)');
legend('dynamical ice','initial ice thickness')

figure
plot(p.phi/pi*180,q);
title(['q, n=' num2str(p.n) '  \eta_0=' sprintf('%g',p.eta0)] );
xlabel('Latitude');
ylabel('Ice flux (m^2/s)');

figure
plot(p.phi/pi*180,S);
title(['S, n=' num2str(p.n) '  \eta_0=' sprintf('%g',p.eta0)] );
xlabel('Latitude');
ylabel('Melting rate (m/s)');

figure;
imagesc(p.phi/pi*180,xi,v'*365*86400);
ylabel('\Delta h/h');
xlabel('latitude');
title(['v (m/year), n=' num2str(p.n)]);
colorbar;
set(gca,'YDir','normal')

figure;
imagesc(p.phi/pi*180,xi,dv_dz');
ylabel('\Delta h/h');
xlabel('latitude');
title(['dv/dz, n=' num2str(p.n)]);
colorbar;
set(gca,'YDir','normal')

%% calculate and plot the effective viscosity 
eta_eff=effective_viscosity(p,T_ice,dv_dz,Tf,Qc);

figure;
imagesc(p.phi/pi*180,xi,log10(eta_eff)');
ylabel('\Delta h/h');
xlabel('latitude');
title(['log10 of the effective viscosity, n=' num2str(p.n)]);
colorbar;
set(gca,'YDir','normal')

figure;
imagesc(p.phi/pi*180,xi,log10(abs(0.5*eta_eff.*dv_dz))');
ylabel('\Delta h/h');
xlabel('latitude');
title(['log10 of the stress (Pa), n=' num2str(p.n)]);
colorbar;
set(gca,'YDir','normal')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_ocean(p,To,Tf);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
plot(p.phi/pi*180,To-273.16,p.phi/pi*180,Tf-273.16);
title(['Freezing and ocean temperatures']);
xlabel('Latitude');
ylabel('Temperature (C)');
legend('T_o','T_f')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_salt(p,Sb);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
plot(p.phi/pi*180,Sb);
title(['Ocean salinity'] );
xlabel('Latitude');
ylabel('Salinity (ppt)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p=set_parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

set(0,'defaulttextfontsize',18); set(0,'defaultaxesfontsize',18);
p.n=3; % non-Newtonian fluid; using regularity formulation 
p.sigma=5.67e-8; % stephan-boltzman constant
p.emissivity=0.94;  % the emissivity of Europa's ice, based on Spencer (1987)
p.alpha=0.62; % albedo of Europa
p.p1=0.033;  % relative time of the eclipses 
p.S0=51;  % the solar constant of Jupiter W/m^2
p.epsilon=3.7/180*pi;  % obliquity 
p.dl=0.01;
p.lambda=(0.5*p.dl):p.dl:(2*pi-0.5*p.dl);  % longitude of Earth around the Sun
p.dt_dlambda=zeros(size(p.lambda));
p.e=0.048;  % eccentricity of Jupiter

% ice parameters 
         %Tf=-3.7019+273.16;  % freezing temperature (in Kelvin)
p.T_transition_brittle=-100+273.16;  % transition temperature (between viscous to brittle ice)
p.Tstar=263; % T* of Cuffey and Paterson (2010) Eq. 3.35
p.delta_Tstar=10; % difference between the melting temperature and
                  % transition temperature, related to T* of Cuffey and Paterson (2010) Eq. 3.35 
p.Qm=6e4; % J/mol, Q- of Cuffey and Paterson (2010) 3.35
p.Qp=11.5e4; % J/mol, Q+ of Cuffey and Paterson (2010) 3.36
p.Ql=p.Qm; % J/mol linear (diffusion) activation energy of Barr and Showman (2009)
%p.Ql=59.4e3; % J/mol linear (diffusion) activation energy of Barr and
%           % Showman (2009)
p.kappaI=1.54e-6; % heat cond. of ice
p.rhoI=917;  % density of ice
p.CpI=2000;  % heat capacity of ice
p.L=3.34e5; % Latent heat of ice 
p.rho0=1040; % density of ocean water
p.mu=p.rhoI/p.rho0; 
p.Q_g=0.05; % geothermal heat flux from silicate part

p.eta0=1e14; % Bar and Showman (2009) value for large (1 mm) and
            % small (0.1 mm) grain size. p.eta0 is between 1e13 to 1e15

if p.n==3
  p.As=3.5e-25;
  p.Es=5;
  p.eta_b=1/(p.Es*p.As);  % based on E* and A* of  Cuffey and Paterson (2010)
                          % E*=5 but can be as large as 10 or even larger
                          % See Table 3.5 of Cuffey and Paterson (2010)
%    p.eta_b=1/(5*3.5e-25);  % based on Cuffey and Paterson (2010)
end
p.R_gas_const=8.3144598;  % universal gas constant (J/mol/K)

p.Radius=1561e3; % radius of Europa
p.g=1.314; % gravity constant of Europa (m/s^2) 

% ocean parameters
% beta=30; % W/m^2/K ocean sensible heat flux coefficient; value of
         % Pollard and Kasting (2005)
p.beta=400; % W/m^2/K ocean sensible heat flux coefficient; based on
          % MITgcm value which is based on Holland and Jenkins (1999)
p.Sb=75; % ocean salinity of 50 ppt
p.mSb=p.Sb; % ocean salinity of 50 ppt
p.CpO=3994; % heat capacity of the ocean in J/kg/K
p.Ho=1e5; % ocean depth in meters
p.kappa_GM=1e-0/p.Radius^2; % GM diffusion constant; in 1/s ; almost no effect

p.h=9.9e3;  % mean ice depth

% time parameters 
p.dy=86400;
p.yr=p.dy*365;
p.dt=10000*p.yr;
p.dt0=p.dt;
p.dt_print=20*p.dt;
p.t_final=15e3*p.yr;

p.flag_ocean=0;
p.flag_salt=0;

p.kappa_h=0e-4/p.Radius^2;  % diffusion constant of h; for numerical reasons. 

p.dphi=1/180*pi;
p.phi=(-89/180*pi-p.dphi/2):p.dphi:(89/180*pi+p.dphi/2); % latitude
p.theta=p.phi+pi/2;  % co-latitude 

p.dxi=0.25e-2;  % xi is (z-zb)/(zs-zb)

% Calculate 
p.simple_convection=1;  % apply simple convection
p.tidal_heating=0;  % if 1 consider tidal heating if zero ignore
                    % tidal heating 

alpha_Hussmann=1.6e-4;
kappa_Hussmann=1.4e-6;
%% Hussmann specified the top temperature according to viscosity
%% that increase by 10 compare to the base value. This lead to
%% T=T_m/(1+ln(10)/A) where A=26 for pure water (Barr and Showman,
%% 2009). This give dT of approximately 22K: 
p.delta_T_over_convective_depths_Hussmann=22;

eta_b_Hussmann=p.eta0;
p.Rac=1000;  % critical Rayleigh number
p.const_Rayleigh=alpha_Hussmann*p.rhoI*p.g* ... 
    p.delta_T_over_convective_depths_Hussmann/(kappa_Hussmann* ...
                                               eta_b_Hussmann);
Ra=p.const_Rayleigh*p.h^3;
p.kappaI_conv_cond_factor=10;
p.kappaI_convection=p.kappaI*p.kappaI_conv_cond_factor; 

p.dz_conv=500;
fprintf(1,'p.kappaI=%g; p.kappaI_convection=%g\n',p.kappaI, p.kappaI_convection);

p=calc_Hmax_via_tidal_heating_fit(p);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ice_bottom_heat_flux,T,Diff_matrix] ...
  =calc_T_ice_convection(p,N,T_bot,T_top,kappa,tidal_heating,z,Diff_matrix)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Solve d/dz(kappa(z)dT/dz)=-tidal_heating(z) subject to T(z_b)=T_bot and T(z_t)=T_top. 
%% Finite difference: 
%% - kappa(k+1) is specified between T points k and k+1.
%% - for k=2:N-1 (kappa(k+1)*(T(k+1)-T(k))-kappa(k)*(T(k)-T(k-1)))/dz^2=-tidal_heating(k)/(rho*Cp)
%% -     k=1:    (kappa(2)  *(T(2)-T(1))  -kappa(1)*(T(1)-T_bot)) /dz^2=-tidal_heating(1)/(rho*Cp)
%% -     k=N:    (kappa(N+1)*(T_top-T(N)) -kappa(N)*(T(N)-T(N-1)))/dz^2=-tidal_heating(N)/(rho*Cp)
%% -     k=0:    T(0)=T_bot
%% -     k=N+1:  T(N+1)=T_top
%% - Need to write this as a matrix equation Diff_matrix*v=RHS for
%%   the vector v=(T(1),...,T(N)). 
%% - The coefficients of T(k), k=2,...,N-1 are:
%% T(1<k<N): 0...0,kappa(k)/dz^2,-(kappa(k)+kappa(k+1))/dz^2,kappa(k+1)/dz^2,0...0
%%             RHS(k)=-tidal_heating(k)/(Cp*rhoI)
%% T(1):       -(kappa(1)+kappa(2))/dz^2,kappa(2)/dz^2,0...0
%%             RHS(2)=-tidal_heating(2)/(Cp*rhoI)-kappa(1)*T_bot/dz^2
%% T(N):     0...0,kappa(N)/dz^2,-(kappa(N)+kappa(N+1))/dz^2
%%             RHS(N)=-tidal_heating(N)/(Cp*rhoI)-kappa(N)*T_top/dz^2

dz=z(2)-z(1);

%% Construct the matrix:
diagonal_above(1:N-1)=kappa(2:N)/dz^2;        % above diagonal
diagonal(1:N)=-(kappa(1:N)+kappa(2:N+1))/dz^2;  % diagonal
diagonal_below(1:N-1)=kappa(2:N)/dz^2;      % below diagonal
%% Indices for creating sparse diffusion matrix: 
if p.flag_first_call_to_calc_T_ice_convection
  p.flag_first_call_to_calc_T_ice_convection=0;
  %% row: [below, diagonal, above]:
  I=[2:N 1:N 1:N-1];
  %% column: [below, diagonal, above]:
  J=[1:N-1 1:N 2:N];
end
entries=[diagonal_below(1:N-1) diagonal(1:N) diagonal_above(1:N-1)];
Diff_matrix=sparse(I,J,entries,N,N);

%% Construct the RHS:
X1=(p.CpI*p.rhoI);
X2=dz^2;
RHS(1)=    -tidal_heating(1)    /X1 -kappa(1)  *T_bot/X2;
RHS(2:N-1)=-tidal_heating(2:N-1)/X1                     ;
RHS(N)=    -tidal_heating(N)    /X1 -kappa(N+1)*T_top/X2;

%% Solve:
T=Diff_matrix\RHS';

%% Calculate bottom flux (positive upward):
ice_bottom_heat_flux=-p.CpI*p.rhoI*kappa(1)*(T(1)-T_bot)/dz;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [T]=calc_T_ice_diffusion_only(Ts,Tf,xi)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Solve for ice temperature assuming diffusion only and therefore
%% a linear profile.

T=xi*(Ts-Tf)+Tf;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function eta_eff=effective_viscosity(p,T_ice,dv_dz,Tf,Qc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate effective viscosity at all latitudes and depths
eta_eff=zeros(p.n_theta,p.n_xi);
for i=1:length(p.theta)
  TT=T_ice(i,:);
  eta_eff(i,:)=p.eta0*exp((p.Qm/(p.R_gas_const*Tf(i)))*(Tf(i)./TT-1 ));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function heavy=heavyside_smooth(x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
heavy=0.5*(tanh(x)+1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p=calc_Hmax_via_tidal_heating_fit(p);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% fit to Nimmo et al 2007:

%% load Nimmo data:
strain=load('strain-rate-Nimmo-et-al-2007.txt');
strain_mean=mean(strain(:,2));
strain_amplitude=max(strain(:,2))-mean(strain(:,2));
A=1.25*strain_mean;
B=1.2*strain_amplitude;
C=0.75;
strain_fit=A-B*(1+cos(2*p.phi)).^C;

omega=2*pi/306806;
mu_E=3.3e9; % Tobie et al 2003
p.eta_max=mu_E/omega;
%%p.H_max=(mu_E/omega)*strain_fit.^2;
p.H_max=(mu_E/omega)*strain_fit.^2; %% XX divide by a factor 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function debugging_plots(p,t,xi,W,Wsols,h_conducting,total_tidal_heating_ice, ...
                    ice_bottom_heat_flux,hh,dh_dtheta,dh_dtheta2,Ts,Tf, ...
                    z_convective_conductive_transition,T_ice,eta_eff, ...
                    tidal_heating,kappa_conv_2d)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ts_sols=calc_Ts(p,(total_tidal_heating_ice+ice_bottom_heat_flux)',Wsols);
Ts_sols_no_tidal_heating=calc_Ts(p,ice_bottom_heat_flux',Wsols);
figure; plot(p.phi*180/pi,Ts, p.phi*180/pi,calc_Ts_new(p,zeros(size(p.phi)),W), ...
             p.phi*180/pi,Ts_sols, p.phi*180/pi,Ts_sols_no_tidal_heating);
legend('Ts','Ts, Q=0','Ts at solstice',['Ts at solstice without tidal ' ...
                    'heating']);title('Surface temperature');   
figure; plot(p.phi*180/pi,Tf),title('Tf');
figure; plot(p.phi*180/pi,h_conducting/1000); title(sprintf('h conducting (km), for t=%g yr',t/p.yr));
figure; plot(p.phi*180/pi, total_tidal_heating_ice, p.phi*180/pi, ...
             ice_bottom_heat_flux, p.phi*180/pi, ...
             total_tidal_heating_ice+ice_bottom_heat_flux); 
title(sprintf('heat fluxes for t=%g yr',t/p.yr)); ...
    xlabel('Latitude');ylabel('Flux W/m^2');
legend('ice tidal heating','ice bottom heat flux', 'total heating');
%figure; plot(p.phi*180/pi,total_tidal_heating_ice); title(sprintf('total tidal heating, for t=%g yr',t/p.yr));
%figure; plot(p.phi*180/pi,ice_bottom_heat_flux); title(sprintf('ice bottom heat flux, for t=%g yr',t/p.yr));
figure; plot(p.phi*180/pi,hh/1000); title(sprintf('ice depth (km), for t=%g yr',t/p.yr));
figure; plot(p.phi*180/pi,dh_dtheta); title(sprintf('dh_dtheta, for t=%g yr',t/p.yr));
figure; plot(p.phi*180/pi,dh_dtheta2); title(sprintf('dh_dtheta2, for t=%g yr',t/p.yr));
figure; plot(p.phi*180/pi,z_convective_conductive_transition/1000); title(sprintf('z_{diff} (km), for t=%g yr',t/p.yr));

figure;
i1=round(p.n_theta/4); plot(T_ice(i1,:),xi,'r');
hold on
i2=round(p.n_theta/2); plot(T_ice(i2,:),xi,'g');
i3=round(3*p.n_theta/4); plot(T_ice(i3,:),xi,'b');
i4=1; plot(T_ice(i4,:),xi,'k');
ylabel('xi');
xlabel('Ice temperature at 4 latitudes');
title(sprintf('T_{ice}, t=%g kyr',t/(p.yr*1000)));
legend(sprintf('lat=%g',p.phi(i1)*180/pi), ...
       sprintf('lat=%g',p.phi(i2)*180/pi),...
       sprintf('lat=%g',p.phi(i3)*180/pi),...
       sprintf('lat=%g',p.phi(i4)*180/pi));

figure; imagesc(p.phi*180/pi,xi,T_ice'); colorbar; title(sprintf('Tice, for t=%g yr',t/p.yr));set(gca,'YDir','normal')
figure; imagesc(p.phi*180/pi,xi,log10(eta_eff)'); colorbar; title(sprintf('log10(eta_{eff}), for t=%g yr',t/p.yr));set(gca,'YDir','normal')
figure; imagesc(p.phi*180/pi,xi,tidal_heating'); colorbar; title(sprintf('tidal heating, for t=%g yr',t/p.yr));set(gca,'YDir','normal')

figure;
i1=round(p.n_theta/4); plot(tidal_heating(i1,:),xi,'r');
hold on
i2=round(p.n_theta/2); plot(tidal_heating(i2,:),xi,'g');
i3=round(3*p.n_theta/4); plot(tidal_heating(i3,:),xi,'b');
ylabel('xi');
xlabel('Tidal heating at 3 latitudes');
title(sprintf('tidal heating, t=%g kyr',t/(p.yr*1000)));
legend(sprintf('lat=%g',p.phi(i1)*180/pi), ...
       sprintf('lat=%g',p.phi(i2)*180/pi),...
       sprintf('lat=%g',p.phi(i3)*180/pi));

figure; imagesc(kappa_conv_2d'); colorbar; title(sprintf('kappa conv, for t=%g yr',t/p.yr));set(gca,'YDir','normal');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fun=ploynomial_for_dvdz(x,polynomial_coefficients) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fun=polyval(polynomial_coefficients,x);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [h,Ts,Tf]=calc_initial_h_Ts_Tf(p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
set(0,'defaulttextfontsize',16); set(0,'defaultaxesfontsize',16);

[p,W]=calc_daily_insolation(p);

Ts=calc_Ts_new(p,p.Q_g*ones(size(p.phi)),W);

% based on Tobie et al. (2003)

xi=0:p.dxi:1; % xi is (z-zb)/(zs-zb)
p.n_xi=length(xi);

Hmax=p.H_max'*ones(1,p.n_xi);

% initial depth guess
h=p.h*ones(size(p.phi));
% Assume 10 km of ice below the sea level (thus don't multiply by the ratio between density of ice and density of water)
Pb=p.g*h*p.rhoI/1e4; % in dBar
% compute the approximate freezing temperature from the MITgcm manual
Tf=273.16+(0.0901-0.0575*p.Sb)-7.61e-4*Pb; % (-3.7019 C)

% compute initial depth 
h=p.rhoI*p.CpI*p.kappaI/p.Q_g*(Tf-Ts);

% perform iteration to find the thickness when the tidal heating is taken
% into account 
for i=1:30
  ggamma=p.Qm/p.R_gas_const./Tf.^2; % constant of the viscosity exponent
  A=p.eta0*exp(-(ggamma.*(Ts-Tf))'*xi);
  H=2*Hmax./(A/p.eta_max+p.eta_max./A);  % tidal heating 
  Qtides=trapz(xi,H,2)'.*h;
  h=p.rhoI*p.CpI*p.kappaI./(p.Q_g+Qtides).*(Tf-Ts);
  Ts=calc_Ts_new(p,p.Q_g+Qtides,W);
  Pb=p.g*h*p.rhoI/1e4; % in dBar
  Tf=273.16+(0.0901-0.0575*p.Sb)-7.61e-4*Pb; % (-3.7019 C)
end
