 |
tgCosmoLongsteps.m - cosmo - front and backend for Markov-Chain Monte Carlo inv… |
 |
git clone git://src.adamsgaard.dk/cosmo |
 |
Log |
|
Files |
|
Refs |
|
README |
|
LICENSE |
|
--- |
|
tgCosmoLongsteps.m (28426B) |
|
--- |
|
1 function [c10Bes,c26Als,c21Nes,c14Cs,lump] = ... |
|
2 gCosmoLongsteps(ErateInt,ErateGla,tDegla,dtGla,dtIdtG,fixed_stuff); |
|
3 %also saves concentration histories in "lump". |
|
4 %gCosmoLongsteps2: Allows computation with fixed intervals |
|
5 %Now, it is possible to compute the forward solution as a linear |
|
6 %approximation: |
|
7 switch fixed_stuff.testmode |
|
8 case 'fast', |
|
9 dtInt = dtGla*dtIdtG; |
|
10 case 'run_advec', %trim model to steplength in time |
|
11 disp('gCosmoLongstep called with mode run_advec. Takes a long time!') |
|
12 dt_advec = fixed_stuff.dt, |
|
13 tDegla = round(tDegla/dt_advec)*dt_advec, |
|
14 dtGla = round(dtGla/dt_advec)*dt_advec, |
|
15 dtInt = round(dtGla*dtIdtG/dt_advec)*dt_advec, |
|
16 case 'linearized' |
|
17 Gpr0 = fixed_stuff.Gpr0; %matrix of partial derivatives at or near m… |
|
18 mpr0 = fixed_stuff.mpr0; %model at or near maximum likelihood soluti… |
|
19 dpr0 = fixed_stuff.dpr0; %proper g(m0) |
|
20 %>........ Model parameters. Her we may switch parameters on and off |
|
21 % fs.mname{1} = 'ErateInt'; |
|
22 % fs.mname{2} = 'ErateGla'; |
|
23 % fs.mname{3} = 'tDegla'; |
|
24 % fs.mname{4} = 'dtGla'; |
|
25 % fs.mname{5} = 'dtIdtG'; |
|
26 for im=1:length(fixed_stuff.mname) %we pack elements of the m-vector: |
|
27 m(im,1) = eval(fixed_stuff.mname{im}); |
|
28 end |
|
29 mpr = m2mpr(m,fixed_stuff); |
|
30 dpr = dpr0 + Gpr0*(mpr-mpr0); |
|
31 d = dpr2d(dpr,fixed_stuff); |
|
32 c10Bes = NaN*ones(size(fixed_stuff.zobs)); |
|
33 c26Als = c10Bes; c21Nes = c10Bes; c14Cs = c10Bes; |
|
34 Nzobs = length(fixed_stuff.zobs); |
|
35 for iNucl = 1:length(fixed_stuff.Nucleides) |
|
36 switch fixed_stuff.Nucleides{iNucl} |
|
37 case '10Be', c10Bes = d((iNucl-1)*Nzobs + (1:Nzobs));%d=[d;cBes(… |
|
38 case '26Al', c26Als = d((iNucl-1)*Nzobs + (1:Nzobs));%d=[d;cAls(… |
|
39 case '21Ne', c21Nes = d((iNucl-1)*Nzobs + (1:Nzobs));%d=[d;cNes(… |
|
40 case '14C', c14Cs = d((iNucl-1)*Nzobs + (1:Nzobs));%d=[d;cCs(:… |
|
41 end |
|
42 end |
|
43 lump.zss = []; |
|
44 lump.ts = []; |
|
45 lump.ExposureTimeSinceNow = []; |
|
46 |
|
47 return |
|
48 otherwise, |
|
49 error('fixed_stuff,.testmode must be ''fast'' or ''run_advec'' or ''… |
|
50 end |
|
51 |
|
52 |
|
53 Period_gla = dtGla; %Names used i previous implementations and below |
|
54 Period_int = dtInt; %Names used i previous implementations and below |
|
55 erate_gla = ErateGla; |
|
56 erate_int = ErateInt; |
|
57 |
|
58 % Rock density in kg/m3 |
|
59 %rho = 2650; |
|
60 rho = fixed_stuff.rho; |
|
61 |
|
62 %Half lives |
|
63 H10=1.39e6; |
|
64 L10=log(2)/H10; |
|
65 H26=0.705e6; |
|
66 L26=log(2)/H26; |
|
67 H14=5730; |
|
68 L14=log(2)/H14; |
|
69 |
|
70 %Absorption depth scale in kg/m2 |
|
71 Tau_spal=1600; |
|
72 Tau_nm = 9900; |
|
73 Tau_fm = 35000; |
|
74 |
|
75 Tau_spal1=1570; |
|
76 Tau_spal2=58.87; |
|
77 |
|
78 Tau_nm1 = 1600; |
|
79 Tau_nm2 = 10300; |
|
80 Tau_nm3 = 30000; |
|
81 |
|
82 Tau_fm1 = 1000; |
|
83 Tau_fm2 = 15200; |
|
84 Tau_fm3 = 76000; |
|
85 |
|
86 |
|
87 a1 = 1.0747; |
|
88 a2 = -0.0747; |
|
89 b1 = -0.050; |
|
90 b2 = 0.845; |
|
91 b3 = 0.205; |
|
92 c1 = 0.010; |
|
93 c2 = 0.615; |
|
94 c3 = 0.375; |
|
95 |
|
96 % ratios between fast muon capture and negative muons |
|
97 nm_rat10 = 0.1070/(0.1070+0.0940); |
|
98 fm_rat10 = 0.0940/(0.1070+0.0940); |
|
99 nm_rat26 = 0.7/(0.7+0.6); |
|
100 fm_rat26 = (0.6/(0.7+0.6)); |
|
101 |
|
102 % unconfirmed ratios |
|
103 nm_rat14 = 0.4/(0.4+0.35); |
|
104 fm_rat14 = 0.35/(0.4+0.35); |
|
105 nm_rat21 = 2.3/(2.3+2.1); |
|
106 fm_rat21 = 2.1/(2.3+2.1); |
|
107 |
|
108 |
|
109 %>>>BHJ: To be used in analytical expressions |
|
110 % L_spal = Tau_spal/rho; %Decay depth, exp(-z/L) |
|
111 % L_nm = Tau_nm/rho; |
|
112 % L_fm = Tau_fm/rho; |
|
113 %"K" in analytical expressions = P10_top_spal, etc. |
|
114 |
|
115 %10Be production |
|
116 % Input fra Kasper |
|
117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - St… |
|
118 P10_top_spal=5.33e3; %atoms/kg/yr |
|
119 P10_top_nm=0.106e3; %atoms/kg/yr |
|
120 P10_top_fm=0.093e3; %atoms/kg/yr |
|
121 |
|
122 if exist('fixed_stuff.be_prod_spall', 'var') == 1 && ... |
|
123 exist('fixed_stuff.be_prod_muons', 'var') == 1 |
|
124 if ~isempty(fixed_stuff.be_prod_spall) && ... |
|
125 ~isempty(fixed_stuff.be_prod_muons) |
|
126 P10_top_spal = fixed_stuff.be_prod_spall; |
|
127 P10_top_nm = nm_rat10*fixed_stuff.be_prod_muons; |
|
128 P10_top_fm = fm_rat10*fixed_stuff.be_prod_muons; |
|
129 end |
|
130 end |
|
131 |
|
132 %Reference values for Kasper |
|
133 %P10_top_spal=5.33e3; %atoms/kg/yr |
|
134 %P10_top_nm=0.106e3; %atoms/kg/yr |
|
135 %P10_top_fm=0.093e3; %atoms/kg/yr |
|
136 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - End |
|
137 |
|
138 %P10_top_spal=5.33e3; %atoms/kg/yr |
|
139 %P10_top_spal=4.76e3; %atoms/kg/yr |
|
140 |
|
141 %P10_top_nm=0.106e3; %atoms/kg/yr |
|
142 %P10_top_nm=0.959e3; %atoms/kg/yr |
|
143 |
|
144 %P10_top_fm=0.093e3; %atoms/kg/yr |
|
145 %P10_top_fm=0.084e3; %atoms/kg/yr |
|
146 |
|
147 %Prodi rate for Corbett's sample GU110 |
|
148 %P10_top_spal=8.04e3; %atoms/kg/yr |
|
149 %P10_top_nm=0.125e3; %atoms/kg/yr |
|
150 %P10_top_fm=0.114e3; %atoms/kg/yr |
|
151 |
|
152 |
|
153 % P10_spal = P10_top_spal*exp(-z*rho/Tau_spal); |
|
154 % P10_nm = P10_top_nm*exp(-z*rho/Tau_nm); |
|
155 % P10_fm = P10_top_fm*exp(-z*rho/Tau_fm); |
|
156 % |
|
157 % P10_total = (P10_spal + P10_nm + P10_fm); |
|
158 |
|
159 %26Al production |
|
160 |
|
161 % Input fra Kasper |
|
162 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - St… |
|
163 P26_top_spal=31.1e3; %atoms/kg/yr |
|
164 P26_top_nm=0.7e3; %atoms/kg/yr |
|
165 P26_top_fm=0.6e3; %atoms/kg/yr |
|
166 |
|
167 if exist('fixed_stuff.al_prod_spall', 'var') == 1 && ... |
|
168 exist('fixed_stuff.al_prod_muons', 'var') == 1 |
|
169 if ~isempty(fixed_stuff.al_prod_spall) && ... |
|
170 ~isempty(fixed_stuff.al_prod_muons) |
|
171 P26_top_spal = fixed_stuff.al_prod_spall; |
|
172 P26_top_nm = nm_rat26*fixed_stuff.al_prod_muons; |
|
173 P26_top_fm = fm_rat26*fixed_stuff.al_prod_muons; |
|
174 end |
|
175 end |
|
176 |
|
177 %Reference values for Kasper |
|
178 %P26_top_spal=31.1e3; %atoms/kg/yr |
|
179 %P26_top_nm=0.7e3; %atoms/kg/yr |
|
180 %P26_top_fm=0.6e3; %atoms/kg/yr |
|
181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - End |
|
182 |
|
183 %P26_top_spal=31.1e3; %atoms/kg/yr |
|
184 %P26_top_nm=0.7e3; %atoms/kg/yr |
|
185 %P26_top_fm=0.6e3; %atoms/kg/yr |
|
186 |
|
187 %P26_top_spal=54.25e3; %atoms/kg/yr |
|
188 %P26_top_nm=1.074e3; %atoms/kg/yr |
|
189 %P26_top_fm=0.923e3; %atoms/kg/yr |
|
190 |
|
191 % P26_spal = P26_top_spal*exp(-z*rho/Tau_spal); |
|
192 % P26_nm = P26_top_nm*exp(-z*rho/Tau_nm); |
|
193 % P26_fm = P26_top_fm*exp(-z*rho/Tau_fm); |
|
194 % |
|
195 % P26_total = (P26_spal + P26_nm + P26_fm); |
|
196 |
|
197 %21Ne production |
|
198 P21_top_spal=20.8e3; %atoms/kg/yr |
|
199 P21_top_nm=0.4e3; %atoms/kg/yr |
|
200 P21_top_fm=0.35e3; %atoms/kg/yr |
|
201 |
|
202 if exist('fixed_stuff.ne_prod_spall', 'var') == 1 && ... |
|
203 exist('fixed_stuff.ne_prod_muons', 'var') == 1 |
|
204 if ~isempty(fixed_stuff.ne_prod_spall) && ... |
|
205 ~isempty(fixed_stuff.ne_prod_muons) |
|
206 P21_top_spal = fixed_stuff.ne_prod_spall; |
|
207 P21_top_nm = nm_rat21*fixed_stuff.ne_prod_muons; |
|
208 P21_top_fm = fm_rat21*fixed_stuff.ne_prod_muons; |
|
209 end |
|
210 end |
|
211 |
|
212 % P21_spal = P21_top_spal*exp(-z*rho/Tau_spal); |
|
213 % P21_nm = P21_top_nm*exp(-z*rho/Tau_nm); |
|
214 % P21_fm = P21_top_fm*exp(-z*rho/Tau_fm); |
|
215 % |
|
216 % P21_total = (P21_spal + P21_nm + P21_fm); |
|
217 |
|
218 %14C production |
|
219 P14_top_spal=14.6e3; %atoms/kg/yr |
|
220 P14_top_nm=2.3e3; %atoms/kg/yr |
|
221 P14_top_fm=2.1e3; %atoms/kg/yr |
|
222 |
|
223 if exist('fixed_stuff.c_prod_spall', 'var') == 1 && ... |
|
224 exist('fixed_stuff.c_prod_muons', 'var') == 1 |
|
225 if ~isempty(fixed_stuff.c_prod_spall) && ... |
|
226 ~isempty(fixed_stuff.c_prod_muons) |
|
227 P14_top_spal = fixed_stuff.c_prod_spall; |
|
228 P14_top_nm = nm_rat14*fixed_stuff.c_prod_muons; |
|
229 P14_top_fm = fm_rat14*fixed_stuff.c_prod_muons; |
|
230 end |
|
231 end |
|
232 |
|
233 % P14_spal = P14_top_spal*exp(-z*rho/Tau_spal); |
|
234 % P14_nm = P14_top_nm*exp(-z*rho/Tau_nm); |
|
235 % P14_fm = P14_top_fm*exp(-z*rho/Tau_fm); |
|
236 % |
|
237 % P14_total = (P14_spal + P14_nm + P14_fm); |
|
238 %>>>>> BHJ: And now the analytical solution |
|
239 %INPUTS: |
|
240 % tau: Decay time of nucleide |
|
241 % ts: End times of intervals. Normally ts(1)=0=now, and ts(2:end)<0, so |
|
242 % that ts is a decreasing vector, going more and more negative. |
|
243 % zs: Present lamina depths below present surface z=0. |
|
244 % ers: ers(i) is the erosion rate in interval between ts(i) and ts(i+1), |
|
245 % Ks: Ks(i) is the production rate in interval between ts(i) and ts(i+1) |
|
246 % L: Penetration depth of radiation in play. dC/dt(z) = K*exp(-z/L) |
|
247 %OUTPUTS: |
|
248 %Css: Concentrations of nucleide. Css(it,iz) applies to ts(it) and zs(iz) |
|
249 %zss,tss: depths and times so that mesh(zss,tss,Css) works.%The call: |
|
250 %[Css,zss,tss] = C_all_intervals(ts,zs,Ks,tau,ers,L) |
|
251 |
|
252 % ts = -fliplr(time); %First element in ts is the start of the model run |
|
253 % zs = fixed_stuff.z; |
|
254 |
|
255 %Setting up Gla/Int timing: |
|
256 |
|
257 switch fixed_stuff.CycleMode |
|
258 case 'FixedC' |
|
259 if isempty('fixed_stuff.C') |
|
260 C=round((2.6e6-tDegla+dtGla*dtIdtG)/(dtGla*(1+dtIdtG))); |
|
261 else |
|
262 C=fixed_stuff.C; |
|
263 end |
|
264 tstart_glacial = -tDegla+Period_int -(C:-1:1)*(Period_gla+Period_int… |
|
265 tend_glacial = tstart_glacial + Period_gla; |
|
266 case 'FixedQuaternary' |
|
267 tQuat=-fixed_stuff.tQuaternary; %time of first glaciation, adjust as… |
|
268 Cfloor=floor((2.6e6-tDegla+dtGla*dtIdtG)/(dtGla*(1+dtIdtG))); |
|
269 try |
|
270 tStartGlaRegular = -tDegla+Period_int -(Cfloor:-1:1)*(Period_gla+Per… |
|
271 catch |
|
272 keyboard |
|
273 end |
|
274 tStartIntRegular = tStartGlaRegular + Period_gla; |
|
275 dtFirstCycle = tStartGlaRegular(1)-tQuat; |
|
276 if dtFirstCycle<=0 %the regular cycles fitted the Quaternary exactly |
|
277 tstart_glacial = tStartGlaRegular; |
|
278 tend_glacial = tStartIntRegular; |
|
279 C = Cfloor; |
|
280 % disp('FixedQuaternary: the regular cycles fitted the Quaternary … |
|
281 else |
|
282 dtFirstGla = dtFirstCycle / (1+dtIdtG); |
|
283 tstart_glacial = [tQuat,tStartGlaRegular]; |
|
284 tend_glacial = [tQuat+dtFirstGla,tStartIntRegular]; |
|
285 C = Cfloor+1; |
|
286 end |
|
287 case 'FixedTimes' |
|
288 if any(isnan(fixed_stuff.tGal)) |
|
289 error(['fixed_stuff.CycleMode=''FixedTimes'' not implemented yet']) |
|
290 else |
|
291 t_start_gla = fixed_stuff.tGlas; %load or compute fixed times of g… |
|
292 t_start_int = fixed_stuff.tInt; %load or compute fixed times of gl… |
|
293 end |
|
294 C = length(t_start_gla); |
|
295 case 'd18OTimes' |
|
296 tStarts = fixed_stuff.tStarts; |
|
297 relExpos = fixed_stuff.relExpos; |
|
298 relExpos(end+1) = 1; % !!!!!! We hereby impose full exposure … |
|
299 %d18Oth should not determine the "start of Holocene": update by tD… |
|
300 tStarts(end+1) = -tDegla; |
|
301 % C = length(t_start_gla); |
|
302 end |
|
303 if strcmp(fixed_stuff.CycleMode,'d18OTimes') |
|
304 tsLong = [-inf,tStarts(:)',0]; |
|
305 is_ints2 = relExpos; %If relExpos==1 then is_ints2 is 1, i.e. it is an… |
|
306 is_ints2 = [1,is_ints2(:)']; %Exposed before the Quaternary |
|
307 else %regular intervals |
|
308 tsLong = sort([-inf,tstart_glacial,tend_glacial,0]); |
|
309 is_ints2 = 0.5*(1 + -(-1).^(1:2*C+1)); %[1 0 1 0 ... 0 1] |
|
310 end |
|
311 ers2 = erate_int*is_ints2 + erate_gla*(1-is_ints2); |
|
312 |
|
313 tau_10Be = 1/L10; |
|
314 tau_26Al = 1/L26; |
|
315 tau_21Ne = inf; |
|
316 tau_14C = 1/L14; |
|
317 |
|
318 L_spal = Tau_spal/rho; %Decay depth, exp(-z/L) |
|
319 L_nm = Tau_nm/rho; |
|
320 L_fm = Tau_fm/rho; |
|
321 |
|
322 % is_ints(1) = 1; %Ice free before simulation period |
|
323 % ers(1)= erate_int; %same as we used yo initiate the advective solution |
|
324 |
|
325 K_10Be_spal = is_ints2*P10_top_spal; |
|
326 K_10Be_nm = is_ints2*P10_top_nm; |
|
327 K_10Be_fm = is_ints2*P10_top_fm; |
|
328 |
|
329 K_26Al_spal = is_ints2*P26_top_spal; |
|
330 K_26Al_nm = is_ints2*P26_top_nm; |
|
331 K_26Al_fm = is_ints2*P26_top_fm; |
|
332 |
|
333 K_21Ne_spal = is_ints2*P21_top_spal; |
|
334 K_21Ne_nm = is_ints2*P21_top_nm; |
|
335 K_21Ne_fm = is_ints2*P21_top_fm; |
|
336 |
|
337 K_14C_spal = is_ints2*P14_top_spal; |
|
338 K_14C_nm = is_ints2*P14_top_nm; |
|
339 K_14C_fm = is_ints2*P14_top_fm; |
|
340 |
|
341 zobs = fixed_stuff.zobs; |
|
342 % tic |
|
343 %Model 10Be: |
|
344 [Css_10Be_spal,zss,tss,ExposureTimeSinceNow] ... |
|
345 = C_all_intervals2(tsLong,zobs,K_10Be_spal,tau_10Be,ers2… |
|
346 [Css_10Be_nm ] = C_all_intervals2(tsLong,zobs,K_10Be_nm ,tau_10Be,ers2… |
|
347 [Css_10Be_fm ] = C_all_intervals2(tsLong,zobs,K_10Be_fm ,tau_10Be,ers2… |
|
348 c10Bes = Css_10Be_spal(:,end)+Css_10Be_nm(:,end)+Css_10Be_fm(:,end); |
|
349 % cBe_analyt = c10Be_prof(1); |
|
350 |
|
351 %Model 26Al: |
|
352 [Css_26Al_spal] = C_all_intervals2(tsLong,zobs,K_26Al_spal,tau_26Al,ers2… |
|
353 [Css_26Al_nm ] = C_all_intervals2(tsLong,zobs,K_26Al_nm ,tau_26Al,ers2… |
|
354 [Css_26Al_fm ] = C_all_intervals2(tsLong,zobs,K_26Al_fm ,tau_26Al,ers2… |
|
355 c26Als = Css_26Al_spal(:,end)+Css_26Al_nm(:,end)+Css_26Al_fm(:,end); |
|
356 % cAl_analyt = c26Al_prof(1); |
|
357 |
|
358 %Model 21Ne: |
|
359 [Css_21Ne_spal] = C_all_intervals2(tsLong,zobs,K_21Ne_spal,tau_21Ne,ers2… |
|
360 [Css_21Ne_nm ] = C_all_intervals2(tsLong,zobs,K_21Ne_nm ,tau_21Ne,ers2… |
|
361 [Css_21Ne_fm ] = C_all_intervals2(tsLong,zobs,K_21Ne_fm ,tau_21Ne,ers2… |
|
362 c21Nes = Css_21Ne_spal(:,end)+Css_21Ne_nm(:,end)+Css_21Ne_fm(:,end); |
|
363 % cNe_analyt = c21Ne_prof(1); |
|
364 |
|
365 %Model 14C: |
|
366 [Css_14C_spal] = C_all_intervals2(tsLong,zobs,K_14C_spal,tau_14C,ers2,L_… |
|
367 [Css_14C_nm ] = C_all_intervals2(tsLong,zobs,K_14C_nm ,tau_14C,ers2,L_… |
|
368 [Css_14C_fm ] = C_all_intervals2(tsLong,zobs,K_14C_fm ,tau_14C,ers2,L_… |
|
369 c14Cs = Css_14C_spal(:,end)+Css_14C_nm(:,end)+Css_14C_fm(:,end); |
|
370 % cC_analyt = c14C_prof(1); |
|
371 |
|
372 if nargout==5 |
|
373 lump.zss = zss; |
|
374 lump.ts = tss(1,:); |
|
375 lump.ExposureTimeSinceNow =ExposureTimeSinceNow; |
|
376 lump.c14Css = Css_14C_spal+Css_14C_nm+Css_14C_fm; |
|
377 lump.c10Bess = Css_10Be_spal+Css_10Be_nm+Css_10Be_fm; |
|
378 lump.c26Alss = Css_26Al_spal+Css_26Al_nm+Css_26Al_fm; |
|
379 lump.c21Ness = Css_21Ne_spal+Css_21Ne_nm+Css_21Ne_fm; |
|
380 end |
|
381 % disp(['analyt time:']); toc |
|
382 |
|
383 switch fixed_stuff.testmode |
|
384 case 'fast', return |
|
385 case 'run_advec', %just go on |
|
386 otherwise, |
|
387 error('fixed_stuff,.testmode must be ''fast'' or ''run_advec''') |
|
388 end |
|
389 |
|
390 % error('fast did not return') |
|
391 %#####################********************################ |
|
392 % Start of optional advective solution |
|
393 %#####################********************################ |
|
394 |
|
395 D = fixed_stuff.D; |
|
396 z = fixed_stuff.z; |
|
397 |
|
398 % Rock density in kg/m3 |
|
399 rho = 2650; |
|
400 |
|
401 %Half lives |
|
402 H10=1.39e6; |
|
403 L10=log(2)/H10; |
|
404 |
|
405 H26=0.705e6; |
|
406 L26=log(2)/H26; |
|
407 |
|
408 H14=5730; |
|
409 L14=log(2)/H14; |
|
410 |
|
411 %Absorption depth scale in kg/m2 |
|
412 Tau_spal=1600; |
|
413 Tau_nm = 9900; |
|
414 Tau_fm = 35000; |
|
415 |
|
416 %>>>BHJ: To be used in analytical expressions |
|
417 % L_spal = Tau_spal/rho; %Decay depth, exp(-z/L) |
|
418 % L_nm = Tau_nm/rho; |
|
419 % L_fm = Tau_fm/rho; |
|
420 %"K" in analytical expressions = P10_top_spal, etc. |
|
421 |
|
422 %10Be production |
|
423 P10_top_spal=5.33e3; %atoms/kg/yr |
|
424 %P10_top_spal=4.76e3; %atoms/kg/yr |
|
425 |
|
426 P10_top_nm=0.106e3; %atoms/kg/yr |
|
427 %P10_top_nm=0.959e3; %atoms/kg/yr |
|
428 |
|
429 P10_top_fm=0.093e3; %atoms/kg/yr |
|
430 %P10_top_fm=0.084e3; %atoms/kg/yr |
|
431 |
|
432 P10_spal = P10_top_spal*exp(-z*rho/Tau_spal); |
|
433 P10_nm = P10_top_nm*exp(-z*rho/Tau_nm); |
|
434 P10_fm = P10_top_fm*exp(-z*rho/Tau_fm); |
|
435 |
|
436 %P10_spal = P10_top_spal*(a1*exp(-z*rho/Tau_spal1)+a2*exp(-z*rho/Tau_spa… |
|
437 %P10_nm = P10_top_nm*(b1*exp(-z*rho/Tau_nm1)+b2*exp(-z*rho/Tau_nm2)+b3*e… |
|
438 %P10_fm = P10_top_fm*(c1*exp(-z*rho/Tau_fm1)+c2*exp(-z*rho/Tau_fm2)+c3*e… |
|
439 |
|
440 P10_total = (P10_spal + P10_nm + P10_fm); |
|
441 |
|
442 %26Al production |
|
443 P26_top_spal=31.1e3; %atoms/kg/yr |
|
444 P26_top_nm=0.7e3; %atoms/kg/yr |
|
445 P26_top_fm=0.6e3; %atoms/kg/yr |
|
446 |
|
447 P26_spal = P26_top_spal*exp(-z*rho/Tau_spal); |
|
448 P26_nm = P26_top_nm*exp(-z*rho/Tau_nm); |
|
449 P26_fm = P26_top_fm*exp(-z*rho/Tau_fm); |
|
450 |
|
451 P26_total = (P26_spal + P26_nm + P26_fm); |
|
452 |
|
453 %21Ne production |
|
454 P21_top_spal=20.8e3; %atoms/kg/yr |
|
455 P21_top_nm=0.4e3; %atoms/kg/yr |
|
456 P21_top_fm=0.35e3; %atoms/kg/yr |
|
457 |
|
458 P21_spal = P21_top_spal*exp(-z*rho/Tau_spal); |
|
459 P21_nm = P21_top_nm*exp(-z*rho/Tau_nm); |
|
460 P21_fm = P21_top_fm*exp(-z*rho/Tau_fm); |
|
461 |
|
462 P21_total = (P21_spal + P21_nm + P21_fm); |
|
463 |
|
464 %14C production |
|
465 P14_top_spal=14.6e3; %atoms/kg/yr |
|
466 P14_top_nm=2.3e3; %atoms/kg/yr |
|
467 P14_top_fm=2.1e3; %atoms/kg/yr |
|
468 |
|
469 P14_spal = P14_top_spal*exp(-z*rho/Tau_spal); |
|
470 P14_nm = P14_top_nm*exp(-z*rho/Tau_nm); |
|
471 P14_fm = P14_top_fm*exp(-z*rho/Tau_fm); |
|
472 |
|
473 P14_total = (P14_spal + P14_nm + P14_fm); |
|
474 |
|
475 |
|
476 %Starting concentrations, equilibrium since eternity |
|
477 %This is computed as integrated part of the analytical solution |
|
478 switch fixed_stuff.Cstart |
|
479 case 'from dat', |
|
480 C10start = fixed_stuff.C10_start; |
|
481 C26start = fixed_stuff.C26_start; |
|
482 C21start = fixed_stuff.C21_start; |
|
483 C14start = fixed_stuff.C14_start; |
|
484 case 'extend interglacial' |
|
485 %Making the start profile a "fm" steady state profile: |
|
486 ers = erate_int; %assume non-glaciated start conditions |
|
487 |
|
488 [C10_steady_spal] = Lal2002eq8(z,P10_top_spal,L10,rho,Tau_spal,e… |
|
489 [C10_steady_nm] = Lal2002eq8(z,P10_top_nm, L10,rho,Tau_nm, ers… |
|
490 [C10_steady_fm] = Lal2002eq8(z,P10_top_fm, L10,rho,Tau_fm, ers… |
|
491 C10start = C10_steady_spal +C10_steady_nm +C10_steady_fm; |
|
492 |
|
493 [C26_steady_spal] = Lal2002eq8(z,P26_top_spal,L26,rho,Tau_spal,e… |
|
494 [C26_steady_nm] = Lal2002eq8(z,P26_top_nm, L26,rho,Tau_nm, ers… |
|
495 [C26_steady_fm] = Lal2002eq8(z,P26_top_fm, L26,rho,Tau_fm, ers… |
|
496 C26start = C26_steady_spal +C26_steady_nm +C26_steady_fm; |
|
497 |
|
498 L21 = 0; %stable |
|
499 [C21_steady_spal] = Lal2002eq8(z,P21_top_spal,L21,rho,Tau_spal,e… |
|
500 [C21_steady_nm] = Lal2002eq8(z,P21_top_nm, L21,rho,Tau_nm, ers… |
|
501 [C21_steady_fm] = Lal2002eq8(z,P21_top_fm, L21,rho,Tau_fm, ers… |
|
502 C21start = C21_steady_spal +C21_steady_nm +C21_steady_fm; |
|
503 |
|
504 [C14_steady_spal] = Lal2002eq8(z,P14_top_spal,L14,rho,Tau_spal,e… |
|
505 [C14_steady_nm] = Lal2002eq8(z,P14_top_nm, L14,rho,Tau_nm, ers… |
|
506 [C14_steady_fm] = Lal2002eq8(z,P14_top_fm, L14,rho,Tau_fm, ers… |
|
507 C14start = C14_steady_spal +C14_steady_nm +C14_steady_fm; |
|
508 case 'zeros' |
|
509 C10start = zeros(size(z)); |
|
510 C26start = C10start; C21start = C10start; C14start = C10start; |
|
511 end |
|
512 C10 = C10start; |
|
513 C26 = C26start; |
|
514 C21 = C21start; |
|
515 C14 = C14start; |
|
516 |
|
517 % Number of glacial-interglacial cycles |
|
518 % C = 1; |
|
519 C = fixed_stuff.C; |
|
520 |
|
521 for i=1:C; |
|
522 gla_init(i) = (Period_gla*i-Period_gla)+(Period_int*i-Period_int); |
|
523 gla_end(i) = (Period_gla*i)+(Period_int*i-Period_int); |
|
524 |
|
525 gla_hist_for((2*i)-1) = gla_init(i); |
|
526 gla_hist_for(2*i) = gla_end(i); |
|
527 end |
|
528 |
|
529 l = 1; %<----------------- beware of l and 1 |
|
530 gla_hist(:,l) = gla_hist_for; |
|
531 |
|
532 |
|
533 |
|
534 |
|
535 time_gla(l) = Period_gla; |
|
536 time_int(l) = Period_int; |
|
537 |
|
538 dt = fixed_stuff.dt; %time step. |
|
539 %Note that in gcosmo_advec.m glaciations must be full multiples of dt. |
|
540 %Here, in gcosmo_analyt this is not required. |
|
541 |
|
542 %Number of timesteps |
|
543 %Bem?rk, her tilf?jes "Holoc?n", s? det an bekvemt v?re tDegla |
|
544 % ts(l) = (gla_hist(end,l)+Period_int)/dt; |
|
545 ts(l) = (gla_hist(end,l)+tDegla)/dt; |
|
546 |
|
547 time(1) = 0; |
|
548 |
|
549 P10_total_decay = P10_total*dt*exp(-L10*dt); |
|
550 |
|
551 P26_total_decay = P26_total*dt*exp(-L26*dt); |
|
552 |
|
553 P21_total_decay = P21_total*dt; |
|
554 |
|
555 P14_total_decay = P14_total*dt*exp(-L14*dt); |
|
556 |
|
557 tic %advec |
|
558 for i=1:ts(l), |
|
559 |
|
560 D10 = C10*L10*dt; |
|
561 S10_int = P10_total_decay(:) - D10(:); |
|
562 S10_gla = -D10; |
|
563 |
|
564 D26 = C26*L26*dt; |
|
565 S26_int = P26_total_decay(:) - D26(:); |
|
566 S26_gla = -D26; |
|
567 |
|
568 D21 = zeros(size(z))'; |
|
569 S21_int = P21_total_decay(:) - D21(:); |
|
570 S21_gla = -D21; |
|
571 |
|
572 D14 = C14*L14*dt; |
|
573 S14_int = P14_total_decay(:) - D14(:); |
|
574 S14_gla = -D14; |
|
575 |
|
576 time(i+1) = dt*i; |
|
577 |
|
578 |
|
579 for j=1:2:length(gla_hist(:,l)) |
|
580 if time(i) >= gla_hist(j,l); |
|
581 erate = erate_gla(l); |
|
582 S10 = S10_gla; |
|
583 S26 = S26_gla; |
|
584 S21 = S21_gla; |
|
585 S14 = S14_gla; |
|
586 is_int = 0; %<<<< BHJ |
|
587 end |
|
588 |
|
589 if time(i) >= gla_hist(j+1,l); |
|
590 erate = erate_int(l); |
|
591 S10 = S10_int; |
|
592 S26 = S26_int; |
|
593 S21 = S21_int; |
|
594 S14 = S14_int; |
|
595 is_int = 1; %<<<< BHJ |
|
596 end |
|
597 end |
|
598 |
|
599 |
|
600 erosion(i) = erate; |
|
601 is_ints(i+1) = is_int; |
|
602 |
|
603 C10 = update_profile(C10,z,S10,erate,dt); |
|
604 C26 = update_profile(C26,z,S26,erate,dt); |
|
605 C21 = update_profile(C21,z,S21,erate,dt); |
|
606 C14 = update_profile(C14,z,S14,erate,dt); |
|
607 |
|
608 % try |
|
609 % C10 = update_profile_dummy(C10,z,S10,erate,dt); %disp(['C10, i=',num2s… |
|
610 % C26 = update_profile_dummy(C26,z,S26,erate,dt); %disp(['C26, i=',num2s… |
|
611 % C21 = update_profile_dummy(C21,z,S21,erate,dt); %disp(['C21, i=',num2s… |
|
612 % C14 = update_profile_dummy(C14,z,S14,erate,dt); %disp(['C14, i=',num2s… |
|
613 % catch |
|
614 % keyboard |
|
615 % end |
|
616 D10_top(i) = D10(1); |
|
617 S10_top(i) = S10(1); |
|
618 C10_top(i) = C10(1); |
|
619 |
|
620 D26_top(i) = D26(1); |
|
621 S26_top(i) = S26(1); |
|
622 C26_top(i) = C26(1); |
|
623 |
|
624 D21_top(i) = D21(1); |
|
625 S21_top(i) = S21(1); |
|
626 C21_top(i) = C21(1); |
|
627 |
|
628 D14_top(i) = D14(1); |
|
629 S14_top(i) = S14(1); |
|
630 C14_top(i) = C14(1); |
|
631 |
|
632 %line(C10,z,'color','r'); |
|
633 |
|
634 %>>>> BHJ: Pack valus for subsequent analytical run |
|
635 %L_spal, L_nm and L_rm were defined above |
|
636 %{ |
|
637 K_10Be_spal(i+1) = P10_top_spal; |
|
638 K_10Be_nm(i+1) = P10_top_nm; |
|
639 K_10Be_fm(i+1) = P10_top_fm; |
|
640 |
|
641 K_26Al_spal(i+1) = P26_top_spal; |
|
642 K_26Al_nm(i+1) = P26_top_nm; |
|
643 K_26Al_fm(i+1) = P26_top_fm; |
|
644 |
|
645 K_21Ne_spal(i+1) = P21_top_spal; |
|
646 K_21Ne_nm(i+1) = P21_top_nm; |
|
647 K_21Ne_fm(i+1) = P21_top_fm; |
|
648 |
|
649 K_14C_spal(i+1) = P14_top_spal; |
|
650 K_14C_nm(i+1) = P14_top_nm; |
|
651 K_14C_fm(i+1) = P14_top_fm; |
|
652 %} |
|
653 ers(i+1) = erate; |
|
654 |
|
655 end; |
|
656 |
|
657 t_time(l) = time((C*Period_gla+(C-1)*Period_int+tDegla)/dt+1); |
|
658 |
|
659 D10_top = D10_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
660 S10_top = S10_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
661 C10_top = C10_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
662 |
|
663 D26_top = D26_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
664 S26_top = S26_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
665 C26_top = C26_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
666 |
|
667 D21_top = D21_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
668 S21_top = S21_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
669 C21_top = C21_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
670 |
|
671 D14_top = D14_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
672 S14_top = S14_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
673 C14_top = C14_top(1:(C*Period_gla+(C-1)*Period_int+tDegla)/dt); |
|
674 |
|
675 %Concentration in top node at the end of the last interglacia |
|
676 eval(['C10_top_p',num2str(l),'= C10_top(end);']); |
|
677 eval(['C26_top_p',num2str(l),'= C26_top(end);']); |
|
678 eval(['C21_top_p',num2str(l),'= C21_top(end);']); |
|
679 eval(['C14_top_p',num2str(l),'= C14_top(end);']); |
|
680 |
|
681 %The output values cNe,cC,cAl,CBe: |
|
682 cBe = C10_top(end); |
|
683 cAl = C26_top(end); |
|
684 cNe = C21_top(end); |
|
685 cC = C14_top(end); |
|
686 |
|
687 disp(['advec time:']); toc |
|
688 tic %analyt |
|
689 %>>>>> BHJ: And now the analytical solution |
|
690 %INPUTS: |
|
691 % tau: Decay time of nucleide |
|
692 % ts: End times of intervals. Normally ts(1)=0=now, and ts(2:end)<0, so |
|
693 % that ts is a decreasing vector, going more and more negative. |
|
694 % zs: Present lamina depths below present surface z=0. |
|
695 % ers: ers(i) is the erosion rate in interval between ts(i) and ts(i+1), |
|
696 % Ks: Ks(i) is the production rate in interval between ts(i) and ts(i+1) |
|
697 % L: Penetration depth of radiation in play. dC/dt(z) = K*exp(-z/L) |
|
698 %OUTPUTS: |
|
699 %Css: Concentrations of nucleide. Css(it,iz) applies to ts(it) and zs(iz) |
|
700 %zss,tss: depths and times so that mesh(zss,tss,Css) works.%The call: |
|
701 %[Css,zss,tss] = C_all_intervals(ts,zs,Ks,tau,ers,L) |
|
702 |
|
703 ts = -fliplr(time); %First element in ts is the start of the model run |
|
704 zs = z; |
|
705 |
|
706 %Timing of constant intervals: |
|
707 %An infinite "interglacial" followed by C glaciations and C interglacials |
|
708 tstart_glacial = -(C:-1:1)*(Period_gla+Period_int); |
|
709 tend_glacial = tstart_glacial + Period_gla; |
|
710 t_bounds = sort([-inf,tstart_glacial,tend_glacial,0]); |
|
711 is_ints2 = 0.5*(1 + (-1).^(1:2*C+1)); %[1 0 1 0 ... 0 1] |
|
712 |
|
713 ers2 = erate_int*is_ints2 + erate_gla*(1-is_ints2); |
|
714 |
|
715 tau_10Be = 1/L10; |
|
716 tau_26Al = 1/L26; |
|
717 tau_21Ne = inf; |
|
718 tau_14C = 1/L14; |
|
719 |
|
720 L_spal = Tau_spal/rho; %Decay depth, exp(-z/L) |
|
721 L_nm = Tau_nm/rho; |
|
722 L_fm = Tau_fm/rho; |
|
723 |
|
724 is_ints(1) = 1; %Ice free before simulation period |
|
725 ers(1)= erate_int; %same as we used yo initiate the advective solution |
|
726 |
|
727 K_10Be_spal = is_ints*P10_top_spal; |
|
728 K_10Be_nm = is_ints*P10_top_nm; |
|
729 K_10Be_fm = is_ints*P10_top_fm; |
|
730 |
|
731 K_26Al_spal = is_ints*P26_top_spal; |
|
732 K_26Al_nm = is_ints*P26_top_nm; |
|
733 K_26Al_fm = is_ints*P26_top_fm; |
|
734 |
|
735 K_21Ne_spal = is_ints*P21_top_spal; |
|
736 K_21Ne_nm = is_ints*P21_top_nm; |
|
737 K_21Ne_fm = is_ints*P21_top_fm; |
|
738 |
|
739 K_14C_spal = is_ints*P14_top_spal; |
|
740 K_14C_nm = is_ints*P14_top_nm; |
|
741 K_14C_fm = is_ints*P14_top_fm; |
|
742 |
|
743 tic |
|
744 %Model 10Be: |
|
745 [Css_10Be_spal,zss,tss] = C_all_intervals(ts,zs,K_10Be_spal,tau_10Be,ers… |
|
746 [Css_10Be_nm ,zss,tss] = C_all_intervals(ts,zs,K_10Be_nm ,tau_10Be,ers… |
|
747 [Css_10Be_fm ,zss,tss] = C_all_intervals(ts,zs,K_10Be_fm ,tau_10Be,ers… |
|
748 c10Be_prof = Css_10Be_spal(:,end)+Css_10Be_nm(:,end)+Css_10Be_fm(:,end); |
|
749 cBe_analyt = c10Be_prof(1); |
|
750 |
|
751 %Model 26Al: |
|
752 [Css_26Al_spal,zss,tss] = C_all_intervals(ts,zs,K_26Al_spal,tau_26Al,ers… |
|
753 [Css_26Al_nm ,zss,tss] = C_all_intervals(ts,zs,K_26Al_nm ,tau_26Al,ers… |
|
754 [Css_26Al_fm ,zss,tss] = C_all_intervals(ts,zs,K_26Al_fm ,tau_26Al,ers… |
|
755 c26Al_prof = Css_26Al_spal(:,end)+Css_26Al_nm(:,end)+Css_26Al_fm(:,end); |
|
756 cAl_analyt = c26Al_prof(1); |
|
757 |
|
758 %Model 21Ne: |
|
759 [Css_21Ne_spal,zss,tss] = C_all_intervals(ts,zs,K_21Ne_spal,tau_21Ne,ers… |
|
760 [Css_21Ne_nm ,zss,tss] = C_all_intervals(ts,zs,K_21Ne_nm ,tau_21Ne,ers… |
|
761 [Css_21Ne_fm ,zss,tss] = C_all_intervals(ts,zs,K_21Ne_fm ,tau_21Ne,ers… |
|
762 c21Ne_prof = Css_21Ne_spal(:,end)+Css_21Ne_nm(:,end)+Css_21Ne_fm(:,end); |
|
763 cNe_analyt = c21Ne_prof(1); |
|
764 |
|
765 %Model 14C: |
|
766 [Css_14C_spal,zss,tss] = C_all_intervals(ts,zs,K_14C_spal,tau_14C,ers,L_… |
|
767 [Css_14C_nm ,zss,tss] = C_all_intervals(ts,zs,K_14C_nm ,tau_14C,ers,L_… |
|
768 [Css_14C_fm ,zss,tss] = C_all_intervals(ts,zs,K_14C_fm ,tau_14C,ers,L_… |
|
769 c14C_prof = Css_14C_spal(:,end)+Css_14C_nm(:,end)+Css_14C_fm(:,end); |
|
770 cC_analyt = c14C_prof(1); |
|
771 |
|
772 disp(['analyt time:']); toc |
|
773 |
|
774 %Now compare results: |
|
775 [cBe,cBe_analyt;... |
|
776 cAl,cAl_analyt;... |
|
777 cNe,cNe_analyt;... |
|
778 cC,cC_analyt] |
|
779 rel_error=[(cBe-cBe_analyt)/cBe;... |
|
780 (cAl-cAl_analyt)/cAl;... |
|
781 (cNe-cNe_analyt)/cNe;... |
|
782 (cC-cC_analyt)/cC] |
|
783 |
|
784 figure; set(gcf,'name','AdvecSurface vs AnalyTrace') |
|
785 subplot(2,2,1) |
|
786 plot(ts(2:end),(Css_10Be_fm(1,2:end)+Css_10Be_nm(1,2:end)+Css_10Be_spal(… |
|
787 hold on |
|
788 plot(ts(2:end),C10_top,'-') |
|
789 title('10Be') |
|
790 legend('Analy','Advec','Location','Northwest') |
|
791 |
|
792 subplot(2,2,2) |
|
793 plot(ts(2:end),(Css_26Al_fm(1,2:end)+Css_26Al_nm(1,2:end)+Css_26Al_spal(… |
|
794 hold on |
|
795 plot(ts(2:end),C26_top,'-') |
|
796 title('26Al') |
|
797 legend('Analy','Advec','Location','Northwest') |
|
798 |
|
799 subplot(2,2,3) |
|
800 plot(ts(2:end),(Css_21Ne_fm(1,2:end)+Css_21Ne_nm(1,2:end)+Css_21Ne_spal(… |
|
801 hold on |
|
802 plot(ts(2:end),C21_top,'-') |
|
803 title('21Ne') |
|
804 legend('Analy','Advec','Location','Northwest') |
|
805 |
|
806 subplot(2,2,4) |
|
807 plot(ts(2:end),(Css_14C_fm(1,2:end)+Css_14C_nm(1,2:end)+Css_14C_spal(1,2… |
|
808 hold on |
|
809 plot(ts(2:end),C14_top,'-') |
|
810 title('14C') |
|
811 legend('Analy','Advec','Location','Northwest') |
|
812 |
|
813 figure; set(gcf,'name','(Adv-Ana)/Adv') |
|
814 subplot(2,2,1) |
|
815 plot(ts(2:end),(Css_10Be_fm(1,2:end)+Css_10Be_nm(1,2:end)+Css_10Be_spal(… |
|
816 title('10Be') |
|
817 |
|
818 subplot(2,2,2) |
|
819 plot(ts(2:end),(Css_26Al_fm(1,2:end)+Css_26Al_nm(1,2:end)+Css_26Al_spal(… |
|
820 hold on |
|
821 title('26Al') |
|
822 |
|
823 subplot(2,2,3) |
|
824 plot(ts(2:end),(Css_21Ne_fm(1,2:end)+Css_21Ne_nm(1,2:end)+Css_21Ne_spal(… |
|
825 hold on |
|
826 title('21Ne') |
|
827 |
|
828 subplot(2,2,4) |
|
829 plot(ts(2:end),(Css_14C_fm(1,2:end)+Css_14C_nm(1,2:end)+Css_14C_spal(1,2… |
|
830 hold on |
|
831 title('14C') |
|
832 |
|
833 figure; set(gcf,'name','Compare c(z,t=0)') |
|
834 subplot(2,2,1) |
|
835 plot(zs,c10Be_prof,'.g'); hold on %analyt |
|
836 plot(zs,C10,'.m'); %advec |
|
837 plot(zs,C10start,'-b'); %advec |
|
838 title('10Be(z), t=0') |
|
839 xlabel('Depth [m]') |
|
840 |
|
841 subplot(2,2,2) |
|
842 plot(zs,c26Al_prof,'.g'); hold on %analyt |
|
843 plot(zs,C26,'.m'); %advec |
|
844 plot(zs,C26start,'-b'); %advec |
|
845 title('26Al(z), t=0') |
|
846 xlabel('Depth [m]') |
|
847 |
|
848 subplot(2,2,3) |
|
849 plot(zs,c21Ne_prof,'.g'); hold on %analyt |
|
850 plot(zs,C21,'.m') %advec |
|
851 plot(zs,C21start,'-b'); %advec |
|
852 title('21Ne(z), t=0') |
|
853 xlabel('Depth [m]') |
|
854 |
|
855 subplot(2,2,4) |
|
856 plot(zs,c14C_prof,'.g'); hold on %analyt |
|
857 plot(zs,C14,'.m'); %advec |
|
858 plot(zs,C14start,'-b'); %advec |
|
859 title('14C(z), t=0') |
|
860 xlabel('Depth [m]') |
|
861 |
|
862 %And now make check plots of (zobs,cBes) etc. |
|
863 subplot(2,2,1) |
|
864 plot(zobs,c10Bes,'*') |
|
865 subplot(2,2,2) |
|
866 plot(zobs,c26Als,'*') |
|
867 subplot(2,2,3) |
|
868 plot(zobs,c21Nes,'*') |
|
869 subplot(2,2,4) |
|
870 plot(zobs,c14Cs,'*') |
|
871 |
|
872 % figure; set(gcf,'name','(Adv-Ana)/Adv') |
|
873 % subplot(2,2,1) |
|
874 % plot(ts(2:end),(Css_14C_fm(1,2:end)+Css_14C_nm(1,2:end)+Css_14C_spal(1… |
|
875 keyboard |
|
876 %erosion = erosion(1:(C*Period_gla+C*Period_int)/dt); |
|
877 |
|
878 %C10_C26_top = -1*zeros((5*max(P_gla)+5*max(P_int))/100,l); |
|
879 %C10_C26_top(:,l) = C10_top./C26_top; |
|
880 %eval(['C10_C26_top_',num2str(l),'= C10_top./C26_top;']); |
|
881 |
|
882 %C14_C10_top(:,l) = C14_top./C10_top; |
|
883 %eval(['C14_C10_top_',num2str(l),'= C14_top./C10_top;']); |
|
884 |
|
885 %C10_C21_top = -1*zeros((5*max(P_gla)+5*max(P_int))/100,l); |
|
886 %C10_C21_top(:,l) = C10_top./C21_top; |
|
887 %eval(['C10_C21_top_',num2str(l),'= C10_top./C21_top;']); |
|
888 |
|
889 %C14_C26_top = -1*zeros((5*max(P_gla)+5*max(P_int))/100,l); |
|
890 %C14_C26_top(:,l) = C14_top./C26_top; |
|
891 %eval(['C14_C26_top_',num2str(l),'= C14_top./C26_top;']); |
|
892 |
|
893 %C26_C21_top = -1*zeros((5*max(P_gla)+5*max(P_int))/100,l); |
|
894 %C26_C21_top(:,l) = C26_top./C21_top; |
|
895 %eval(['C26_C21_top_',num2str(l),'= C26_top./C21_top;']); |
|
896 |
|
897 %C14_C21_top = -1*zeros((5*max(P_gla)+5*max(P_int))/100,l); |
|
898 %C14_C21_top(:,l) = C14_top./C21_top; |
|
899 %eval(['C14_C21_top_',num2str(l),'= C14_top./C21_top;']); |
|
900 |
|
901 |
|
902 %eval(['erosion_hist_',num2str(l),'= erosion(:);']); |
|
903 return |
|
904 end |
|
905 |
|
906 |
|
907 function C_out = update_profile_dummy(C,z,S,erate,dt); |
|
908 C_out=C+S; |
|
909 |
