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	gforget	1.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