tgCosmoLongsteps.m - cosmo - front and backend for Markov-Chain Monte Carlo inversion of cosmogenic nuclide concentrations
(HTM) git clone git://src.adamsgaard.dk/cosmo
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) 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 maximum likelihood solution
18 mpr0 = fixed_stuff.mpr0; %model at or near maximum likelihood solution
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 ''linearized''')
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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - Start
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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - Start
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 desired
268 Cfloor=floor((2.6e6-tDegla+dtGla*dtIdtG)/(dtGla*(1+dtIdtG)));
269 try
270 tStartGlaRegular = -tDegla+Period_int -(Cfloor:-1:1)*(Period_gla+Period_int);
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 exactly')
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 glaciation starts
292 t_start_int = fixed_stuff.tInt; %load or compute fixed times of glaciation ends
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 during the Holocene
299 %d18Oth should not determine the "start of Holocene": update by tDegla
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 interglacial
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,L_spal);
346 [Css_10Be_nm ] = C_all_intervals2(tsLong,zobs,K_10Be_nm ,tau_10Be,ers2,L_nm );
347 [Css_10Be_fm ] = C_all_intervals2(tsLong,zobs,K_10Be_fm ,tau_10Be,ers2,L_fm );
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,L_spal);
353 [Css_26Al_nm ] = C_all_intervals2(tsLong,zobs,K_26Al_nm ,tau_26Al,ers2,L_nm );
354 [Css_26Al_fm ] = C_all_intervals2(tsLong,zobs,K_26Al_fm ,tau_26Al,ers2,L_fm );
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,L_spal);
360 [Css_21Ne_nm ] = C_all_intervals2(tsLong,zobs,K_21Ne_nm ,tau_21Ne,ers2,L_nm );
361 [Css_21Ne_fm ] = C_all_intervals2(tsLong,zobs,K_21Ne_fm ,tau_21Ne,ers2,L_fm );
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_spal);
367 [Css_14C_nm ] = C_all_intervals2(tsLong,zobs,K_14C_nm ,tau_14C,ers2,L_nm );
368 [Css_14C_fm ] = C_all_intervals2(tsLong,zobs,K_14C_fm ,tau_14C,ers2,L_fm );
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_spal2));
437 %P10_nm = P10_top_nm*(b1*exp(-z*rho/Tau_nm1)+b2*exp(-z*rho/Tau_nm2)+b3*exp(-z*rho/Tau_nm3));
438 %P10_fm = P10_top_fm*(c1*exp(-z*rho/Tau_fm1)+c2*exp(-z*rho/Tau_fm2)+c3*exp(-z*rho/Tau_fm3));
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,ers);
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,ers);
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,ers);
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,ers);
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=',num2str(i)])
610 % C26 = update_profile_dummy(C26,z,S26,erate,dt); %disp(['C26, i=',num2str(i)])
611 % C21 = update_profile_dummy(C21,z,S21,erate,dt); %disp(['C21, i=',num2str(i)])
612 % C14 = update_profile_dummy(C14,z,S14,erate,dt); %disp(['C14, i=',num2str(i)])
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,L_spal);
746 [Css_10Be_nm ,zss,tss] = C_all_intervals(ts,zs,K_10Be_nm ,tau_10Be,ers,L_nm );
747 [Css_10Be_fm ,zss,tss] = C_all_intervals(ts,zs,K_10Be_fm ,tau_10Be,ers,L_fm );
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,L_spal);
753 [Css_26Al_nm ,zss,tss] = C_all_intervals(ts,zs,K_26Al_nm ,tau_26Al,ers,L_nm );
754 [Css_26Al_fm ,zss,tss] = C_all_intervals(ts,zs,K_26Al_fm ,tau_26Al,ers,L_fm );
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,L_spal);
760 [Css_21Ne_nm ,zss,tss] = C_all_intervals(ts,zs,K_21Ne_nm ,tau_21Ne,ers,L_nm );
761 [Css_21Ne_fm ,zss,tss] = C_all_intervals(ts,zs,K_21Ne_fm ,tau_21Ne,ers,L_fm );
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_spal);
767 [Css_14C_nm ,zss,tss] = C_all_intervals(ts,zs,K_14C_nm ,tau_14C,ers,L_nm );
768 [Css_14C_fm ,zss,tss] = C_all_intervals(ts,zs,K_14C_fm ,tau_14C,ers,L_fm );
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(1,2:end)),'.')
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(1,2:end)),'.')
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(1,2:end)),'.')
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:end)),'.')
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(1,2:end))./C10_top,'.')
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(1,2:end))./C26_top,'.')
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(1,2:end))./C10_top,'.')
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:end))./C14_top,'.')
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,2:end))./C14_top,'.g')
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