toy_models/ice_thick_distrib/ice_remap.m

function [V,C,H]=ice_remap(v,c,h);
%(Lipscomb's approach using linear distributions; simplified.)
%
%object: remap ice categories after thermodynamical ice growth/melt
%       (v increments) may result in ice categories outgrowing
%       their limits ( v/c>h(n) or v/c<h(n-1) ). Here we assume a
%   homogeneous distribution (linear more generally; see Lipscomb 01)
%   around the updated h1=v/c and compute the overlapping integral 
%   within the prescribed neighbor category limits. The result (dC&dV) 
%   is then moved to that neighboring category (thicker or thinner).
%inputs:        v is the ice volume in each cat
%               c is the ice concentration in each cat
%               h is the cat thickness limits
%outputs:       V,C,H are the recasted vol.,conc., and thick.

V=v; C=c; H=V./c;
nh=size(V,1);
nx=size(V,2);

%0) invert formula for new distribution
%======================================

dh=0.25;

c1=c(2:end-1,:);%category concentration
v1=v(2:end-1,:);%category volume
h1=v1./c1;%moved category center
HL1=h1-dh/2;%moved lower limit
HR1=h1+dh/2;%moved upper limit

HL0=h(1:end-2,:);%fixed lower limit
HR0=h(2:end-1,:);%fixed upper limit
h0=(HL0+HR0)/2;%fixed category center

etaR=HR1-HL1;
etaN=h1-HL1;
g1=12*c1./(etaR.^3).*(etaN-etaR/2);
g0=6*c1./(etaR.^2).*(2/3*etaR-etaN);
%new distribution is 0 outside [HL1 HR1] and otherwise is
%g(h)=g0+g1*(h-HL1);
%
%note: as compared with Lipscomb 01, I simplified the 
%computation of lower/upper boundaries, which implies that
%the distributions are homogeneous -- so g(h)=g0 in fact.
%after simplification:
g1=0*ones(size(etaR)); g0=c1/dh;

%1) recast to thinner categories
%===============================

%1.1) compute integrals between HR0 and HR1, which
%   will be transfered to the next category
dC0=(g0-g1.*HL1).*(HL0-HL1)+1/2*g1.*(HL0.^2-HL1.^2);
dV0=1/2*(g0-g1.*HL1).*(HL0.^2-HL1.^2)+1/3*g1.*(HL0.^3-HL1.^3);

%avoid leaving exponentially small numbers : full recasting
tmp1=V(2:end-1,:); tmp2=C(2:end-1,:);
ii=find(tmp1<1e-4|tmp2<1e-4);
dV0(ii)=tmp1(ii); dC0(ii)=tmp2(ii);

%avoid category jump (HL1<HL0-dh) that leads to 
%dC0>C or dV0>V and ultimately to negative values
tmp1=V(2:end-1,:); tmp2=C(2:end-1,:);
ii=find(HL1<HL0-dh);
dV0(ii)=tmp1(ii); dC0(ii)=tmp2(ii);

%1.2) mask out points where no remapping
ii=find(HL0==0|HL1>=HL0|c1==0|v1==0); dC0(ii)=0; dV0(ii)=0;

%1.3) time step V and C:
V(2:end-1,:)=V(2:end-1,:)-dV0;
V(1:end-2,:)=V(1:end-2,:)+dV0;

C(2:end-1,:)=C(2:end-1,:)-dC0;
C(1:end-2,:)=C(1:end-2,:)+dC0;

%2) recast to thicker categories
%===============================

%2.1) compute integrals between HR0 and HR1, which
%   will be transfered to the next category
dC1=(g0-g1.*HL1).*(HR1-HR0)+1/2*g1.*(HR1.^2-HR0.^2);
dV1=1/2*(g0-g1.*HL1).*(HR1.^2-HR0.^2)+1/3*g1.*(HR1.^3-HR0.^3);

%avoid leaving exponentially small numbers : full recasting
tmp1=V(2:end-1,:); tmp2=C(2:end-1,:);
ii=find(tmp1<1e-4|tmp2<1e-4);
dV1(ii)=tmp1(ii); dC1(ii)=tmp2(ii);

%avoid category jump (HR1>HR0-dh) that leads to 
%dC0>C or dV0>V and ultimately to negative values
tmp1=V(2:end-1,:); tmp2=C(2:end-1,:);
ii=find(HR1>HR0+dh);
dV0(ii)=tmp1(ii); dC0(ii)=tmp2(ii);

%2.2) mask out points where no remapping
ii=find(HR1<=HR0|c1==0|v1==0); dC1(ii)=0; dV1(ii)=0;

%2.3) time step V and C:
V(2:end-1,:)=V(2:end-1,:)-dV1;
V(3:end,:)=V(3:end,:)+dV1;

C(2:end-1,:)=C(2:end-1,:)-dC1;
C(3:end,:)=C(3:end,:)+dC1;
    
1	function [V,C,H]=ice_remap(v,c,h);
2	%(Lipscomb's approach using linear distributions; simplified.)
3	%
4	%object: remap ice categories after thermodynamical ice growth/melt
5	% (v increments) may result in ice categories outgrowing
6	% their limits ( v/c>h(n) or v/c<h(n-1) ). Here we assume a
7	% homogeneous distribution (linear more generally; see Lipscomb 01)
8	% around the updated h1=v/c and compute the overlapping integral
9	% within the prescribed neighbor category limits. The result (dC&dV)
10	% is then moved to that neighboring category (thicker or thinner).
11	%inputs: v is the ice volume in each cat
12	% c is the ice concentration in each cat
13	% h is the cat thickness limits
14	%outputs: V,C,H are the recasted vol.,conc., and thick.
15
16	V=v; C=c; H=V./c;
17	nh=size(V,1);
18	nx=size(V,2);
19
20	%0) invert formula for new distribution
21	%======================================
22
23	dh=0.25;
24
25	c1=c(2:end-1,:);%category concentration
26	v1=v(2:end-1,:);%category volume
27	h1=v1./c1;%moved category center
28	HL1=h1-dh/2;%moved lower limit
29	HR1=h1+dh/2;%moved upper limit
30
31	HL0=h(1:end-2,:);%fixed lower limit
32	HR0=h(2:end-1,:);%fixed upper limit
33	h0=(HL0+HR0)/2;%fixed category center
34
35	etaR=HR1-HL1;
36	etaN=h1-HL1;
37	g1=12c1./(etaR.^3).(etaN-etaR/2);
38	g0=6c1./(etaR.^2).(2/3*etaR-etaN);
39	%new distribution is 0 outside [HL1 HR1] and otherwise is
40	%g(h)=g0+g1*(h-HL1);
41	%
42	%note: as compared with Lipscomb 01, I simplified the
43	%computation of lower/upper boundaries, which implies that
44	%the distributions are homogeneous -- so g(h)=g0 in fact.
45	%after simplification:
46	g1=0*ones(size(etaR)); g0=c1/dh;
47
48	%1) recast to thinner categories
49	%===============================
50
51	%1.1) compute integrals between HR0 and HR1, which
52	% will be transfered to the next category
53	dC0=(g0-g1.HL1).(HL0-HL1)+1/2g1.(HL0.^2-HL1.^2);
54	dV0=1/2(g0-g1.HL1).(HL0.^2-HL1.^2)+1/3g1.*(HL0.^3-HL1.^3);
55
56	%avoid leaving exponentially small numbers : full recasting
57	tmp1=V(2:end-1,:); tmp2=C(2:end-1,:);
58	ii=find(tmp1<1e-4\|tmp2<1e-4);
59	dV0(ii)=tmp1(ii); dC0(ii)=tmp2(ii);
60
61	%avoid category jump (HL1<HL0-dh) that leads to
62	%dC0>C or dV0>V and ultimately to negative values
63	tmp1=V(2:end-1,:); tmp2=C(2:end-1,:);
64	ii=find(HL1<HL0-dh);
65	dV0(ii)=tmp1(ii); dC0(ii)=tmp2(ii);
66
67	%1.2) mask out points where no remapping
68	ii=find(HL0==0\|HL1>=HL0\|c1==0\|v1==0); dC0(ii)=0; dV0(ii)=0;
69
70	%1.3) time step V and C:
71	V(2:end-1,:)=V(2:end-1,:)-dV0;
72	V(1:end-2,:)=V(1:end-2,:)+dV0;
73
74	C(2:end-1,:)=C(2:end-1,:)-dC0;
75	C(1:end-2,:)=C(1:end-2,:)+dC0;
76
77	%2) recast to thicker categories
78	%===============================
79
80	%2.1) compute integrals between HR0 and HR1, which
81	% will be transfered to the next category
82	dC1=(g0-g1.HL1).(HR1-HR0)+1/2g1.(HR1.^2-HR0.^2);
83	dV1=1/2(g0-g1.HL1).(HR1.^2-HR0.^2)+1/3g1.*(HR1.^3-HR0.^3);
84
85	%avoid leaving exponentially small numbers : full recasting
86	tmp1=V(2:end-1,:); tmp2=C(2:end-1,:);
87	ii=find(tmp1<1e-4\|tmp2<1e-4);
88	dV1(ii)=tmp1(ii); dC1(ii)=tmp2(ii);
89
90	%avoid category jump (HR1>HR0-dh) that leads to
91	%dC0>C or dV0>V and ultimately to negative values
92	tmp1=V(2:end-1,:); tmp2=C(2:end-1,:);
93	ii=find(HR1>HR0+dh);
94	dV0(ii)=tmp1(ii); dC0(ii)=tmp2(ii);
95
96	%2.2) mask out points where no remapping
97	ii=find(HR1<=HR0\|c1==0\|v1==0); dC1(ii)=0; dV1(ii)=0;
98
99	%2.3) time step V and C:
100	V(2:end-1,:)=V(2:end-1,:)-dV1;
101	V(3:end,:)=V(3:end,:)+dV1;
102
103	C(2:end-1,:)=C(2:end-1,:)-dC1;
104	C(3:end,:)=C(3:end,:)+dC1;
105