-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathLoscardCAPerm1BU.m
4409 lines (3928 loc) · 156 KB
/
LoscardCAPerm1BU.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
%--------------------------------------
%12/22/2016:
% Thing to do:
% 1) Rewrite GEOCARB so that it takes org C carbon as a driving parameter
% instead of d13C. So you first run it the way it is originally set up to
% get the d13c curve. Next, you take the org C value and rewrite the model
% so that d13C is now predicted. If everything done correctly, this should
% give us a d13C curve identical to the original one
% 2) Now choose a starting year (60 or 50Ma) and set the fluxes in LOSCAR
% to be exactly the same to the ones in the given year in GEOCARB so that
% both models have the same steady-states.
% 3) Run GEOCARB for a million years and supply the data to LOSCAR. LOSCAR
% will calculate a new org C burial, feed this burial to GEOCARB and
% run it for a million year period and calculate all other fluxes and
% supply them to LOSCAR. Keep repeating this for the duration of simulation
% 4) Because the run is going to be over a long time period need to include
% the change in Ca and Mg over this time period because the ratio of the
% two will have effect on dissociation constants. Even though Ca is
% predicted in the model right the predicted Ca doesn't match the
% observation. So need to get both Ca and Mg data over this time period and
% prescribed them instead
% file: Loscar.m
%
% LOSCAR Model: Long-term Ocean-atmosphere-Sediment
% CArbon cycle Reservoir Model
%
% ocean box model (+1 for atmosphere)
%
%
% 0
%
% updates:
% 05/19/10 44Ca tracer added but only for CAflag=0 &
% fsed && ftys && fdox =1
%
% 07/06/08 A few comments added
%
% 02/06/08 ffflag results are slightly different
% for MM-O2 from those saved Oct 2007.
% Set KMMOX = 0.0.
%
% 10/27/07 Michaelis-Menton for Oxygen included.
%
% 10/20/07 Indices of initial f0 (sediments, load)
% corrected in Y0 when oxygen is included
% (caused f>1 during erosion)
%
% 10/20/07 TDflag: derivs with dYflag = 1 now called with
% TCvt because TCv can't be updated (odeXX not called).
% Shouldn't be an issue for TDflag = 0 as
% TCv IS updated.
%
% 05/13/07 constraints on CBl and d13CBl input
% ccdA(1)-min(ccdA)
% kk=150; min(d13c(kk,1))-d13c(1,1)
%
% 01/29/07 Tethys mv=3.5, rrain=
% 01/28/07 Oxygen included. Lots of changes:
% hs, TH, TT, rrain
% 01/22/07 [CO3=] grad: co3tv(jt,9)/co3tv(jt,7)
%
% 12/00/06 sedrate/phi (WRONG!) in old sed model
% new is OK!
% 07/13/06 Temp for co3sat corrected
% 07/05/06 rcak for Ca/Cam corrected (Tethys)
% 07/01/06 Version B runs. Nearly same results.
% 06/26/06 New sediment model (Version B)
% 05/27/06 New erosion
% 05/03/06 Pac CCD shallowed (x in THmfun)
% 03/23/06 Effect of Ca/Mg on K's included
% 03/12/06 Adjust dissolution to Ca. nc 2.4
% 03/11/06 shelf/deep rain C13 missing. done
% 03/08/06 mix/biopump changed
% 03/04/06 Millero sat, Ca = 20, 1000 ppmv
% 03/01/06 H-Lat mix & 13Cin changed for basin-d13C
% 01/18/06 switch TH SO-NP-SO
% 12/29/05 Tethys ocean+sed complete
% 12/26/05 P-Error in dafunPE corrected
% 12/24/05 Tethys
% 12/10/05 Tuning good with water column diss
% 12/08/05 Tuning
% 11/26/05 10-box: rhos, phi = phi(fc)
% 11/22/05 03-box: DEQ in df/dt
% 11/05/05 porosity, phi = phi(fc)
% 10/27/05 error fixed (rsed could become < 0)
% 10/25/05 C13 10-box ocean + sediment
% 10/23/05 C13 10-box ocean
% 10/21/05 C13 3-box sediment
% 10/18/05 resumed
% 05/26/05 10-box (c,a,p) + sed complete
% 05/11/05 new file
%
%
%--------------------------------------
solflag = 1; % 0: skip solver/load
if (solflag == 1)
clear all all;
solflag = 1;
end;
logax = 0; % plot: log axes on/off
axx = [-0.5e5 2.0e5]; % x-axes limits (time)
%axx = [000 1000];
global myflag kasflag Fem20 kt tst Dtst yst dYst stflag ...
Cam Ca Mgm Mg y2s;
kt = 1; % counter time step
Cam = 10.3e-3; % 10.3 (mol/kg) Calcium modern
Mgm = 53.0e-3; % 53.0 (mol/kg) Magnesium modern
Ca = Cam;
Mg = Mgm;
y2s = 3600.*24.*365.; % year to seconds
myflag = 01;
% 01: 10-box model + N sediment + Tethys
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 10 BOX + sediment + Tethys
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if myflag == 01;
if solflag == 01;
%========================================
% global variables:
% known to Loscar10 and LoscarDif10
%========================================
global Aoc rho Nb V A TH0 TH TS mv mhd kasv kkv TCv Sv Pv gp tA tI ...
phflag k1k2flag EPH fEPL rrain ep ept ec ect REDPC REDNC REDO2C eI ...
Ns asvA asvI asvP zv m2kg m2kgca rhos frrf frrfca phic hs FiN ...
phivtA phivtI phivtP phi0 phi1 gam ...
phiiA phiiI phiiP phiiAca phiiIca phiiPca phiiT phiiTca dsv ...
Fpr frain rsedv rburvtA rburvtI rburvtP dissvtA ...
dissvtI dissvtP FprtA FprtI FprtP it co3satv ...
Kd KS cst dsflag nc dYflag nu fsed ...
Fpr13tA Fpr13tI Fpr13tP Fpr44tA Fpr44tI Fpr44tP Fint Fin13t Fin44t Fpr44tT FSit FSi13t FSi44t...
fc0A fc0I fc0P fca0A fca0I fca0P f13c0A f13c0I f13c0P ...
dissv13tA dissv13tI dissv13tP ...
dissv44tA dissv44tI dissv44tP...
co3s0 as zs0 klid nli Rst epsp epspca Rin FiN13 ...
BlFlag CBl RBl kb FVC FVC13 Rvc pCSi nSi nCC Focb0 Fkg13 ...
ftys nOC TT asvT kliT fc0T fca0T f13c0T rsedvtT dissvtT dissv13tT dissv44tT ...
FprtT Fpr13tT Fpr44tT phivtT fcon swcon ...
TCv0 TCvt ntL ntH fsh fshI fshP fdpv fshT nshT ...
ffflag tem em RlsCtv ...
fdox KMMOX vask DTS DTS2 ts3 DTS3 DTS4 ...
k1st tfinal TDflag omegCSvt omegASvt ...
THt mv0 mhd0 oxA CAvflag FSichck Finchck kspCHCK...
CHECK1 CHECK2 CHECK3 CHKFin CHKFSi CHCKbioL CHCKbioI CHCKbioD ...
CHKcarB1 CHKcarB2 CHKcarB3 CHKcar1 CHKcar2 CHKcar3 CHKcar4 CAv g Cinpc Cadv slj avCALC AJDE Voc...
RinCA RcasAcheck RsAcheck f13cvAchck f44cavAchck fcvAchck CHKcarC1 CHKcarC2 CHKcarC3...
CHKcarD1 CHKcarD2 CHKcarD3 CHKcarD4 Rbchck Rbcachck...
rburvtAca rburvtIca rburvtPca rsedvtTca phivtAca phivtIca phivtPca...
FprtAca FprtIca FprtPca dissvtTca phivtTca FprtTca EPLv EXLv ECLv...
fwcv fwsv...
fwsv1 fwgv fmgv pco2gca tgc kgc fbgv fmcv FSi0 fbcv pco20 Focw0 fbch fbbv...
frkc fekc fdkc flakc FERT ACT RUN gamma fGG FiN0 Rkg LTflag acv ybbv ybv ycv epspv wcvtA wcvtI wcvtP wcvtT...
dissvtAca dissvtIca dissvtPca epspcaV inorgF epspF biopF Fpw0...
fEPLv EPHv Rstca dincaV epsSensF runtime Fopb0 Ffep0 Fcap0 oxicf0 anoxf0 Pfeed chck1st Focbv...
EALvv ECALvv ECAHv EAHv dissCavv EPLvv oIv oIpv PPLvv Ffepv Fcapv counter PPHv Pscenario O0...
Floegel rREGv REDPCv meanSpco2v;
counter = 0;
% Loscar-GEOCARB long-term stuff
LTflag = 0; % Flag to switch Long-Term, LOSCAR-GEOCARB run
chck1st=-1;
kgc=0;
if(LTflag)
pco20=300; %default 300 ppm
FERT=0.40;
ACT=0.05;
[ybbv,ybv,ycv,aav,bbv,ccv,fGG,fgkc,frkc,fekc,fdkc,flakc,fbch,fbbv,fwcv,fwsv1,fwsv2,fbcv,fbgv,rco2,pco2gca,fwgv,fmcv,fmgv,ggc,cgc,dg,dc]=myGeoCarbMODULE1(1);
end
% Initial year in Ma
tgc=59;
if(LTflag)
% The length of the run in years
runtime=55e6;
else
runtime=4.e5; %4.e5 end-Permian, 2e5 PETM
end
%%%%Payne data set for d44Ca and d13C end-Permian
pd = csvread('PayneData.csv');
stel = pd(:,1); % stratigraphic elevation in meters
d44cap = pd(:,2);
d13cp = pd(:,4);
depth = [0 160]; % used for calculating time
for i=1:length(stel)
time(i) = interp1(depth, [0 1] , stel(i));
end
% +++++ Edit this line to add correct path to solver ! +++++ %
% *****LINES CHANGED TO GET PETM set-up with Ca= 10.3; 292 540 659 and save
addpath('myode/');
% addpath('/home/komar/Dropbox/geocarb/Loscar d44Ca/myode');
CBl = 3000.e15; % Blast end=Perm 13200, 43200, 10000
d13CBl = -50.; % Blast d13C -55 -60 -34 n-50 -7. -5 for Permian
% set DTS - duration of the initial carbon input pulse
% this for BlFlag = 2;
DTS = 5e3; % 6e3 years for PETM,1e5 and 0.5e5 end-Permian
nSi0 = 0.2; % 0.2
nC0 = 0.4; % 0.4
nSi = 0.2; % 0.2 0.3
nCC = 0.4; % 0.4 0.3 1.0 0.5
plotflag = 1; % plot results
CAvflag =2; % 2:calcium is a tracer and changes,1:calcium is a tracer but stays constant
% 0: Ca is not a tracer
Pfeed = 1; % 0 no Phosphate cycle, 1 phosphate cycle incuded
% note that this variable is also used as a multilying
% constant in TA calculations so that when there is no
% PO4 cycle included and there is no decrease in
% remineralization included, TA is not affected by PO4
Pscenario = 1; % 0 - default Phosphate weathering scenario, po4bf 0.5%
% 1 - fraction of P buried twice, 1%;
% 2 - fraction of P buried halved, 0.25%;
% 3 - fraction of P 0.5%, ocbf 0.5%;
% 4 - fraction of P 1.0%, ocbf 0.5%;
% 5 - fraction of P 0.25%, ocbf 0.5%;
% 6 - fraction of P 0.5%, ocbf 2.0%;
% 7 - fraction of P 1.0%, ocbf 2.0%;
% 8 - fraction of P 0.25%, ocbf 2.0%;
% 9 - fraction of P 0.5%, ocbf 4.0%;
% 10 - fraction of P 1.0%, ocbf 4.0%;
% 11 - fraction of P 0.25%, ocbf 4.0%;
% 12 - fraction of P 2.0%, ocbf 2.0%;
% 13 - the same as 7 but fEPL =0.7 and
% rrain = 7.0 in order to get a deeper
% CCD of 4km in Atl and CaCO3 constant
%
% POx0 - default Pscenarion but oxA = 0;
% POx1 - oxA = 1;
% Pfloegel0 - default Fluegel with oxA =0.4
% Pfloegel1 - default Fluegel with oxA =0.0
% Pfloegel2 - default Fluegel with oxA =1
Floegel = 0; % 0: Slomp parameterizatiom, >0 Floegel-Wallmann param % Pfloegel0
% 1: po4bf0=0.02, ocbf0 = 0.01; % Pfloegel3
% 2: po4bf0=0.02, ocbf0 = 0.02; % Pfloegel4
% 3: po4bf0=0.02, ocbf0 = 0.005;% Pfloegel5
% 4: po4bf0=0.04, ocbf0 = 0.01; % Pfloegel6
% 5: po4bf0=0.01, ocbf0 = 0.01; % Pfloegel7
if(Pscenario)
Floegel = 0;
end
% Specifying which initial file to load and where to save
if(Floegel)
Pscenario = 0;
filepath = ['dat/PETM44CaP/'];
folder = ['Floegel'];
fnum = num2str(Floegel);
else
filepath = ['dat/PETM44CaP/'];
folder = ['Psc'];
fnum = num2str(Pscenario);
end
saveOutput = 0;
savf = 0; % save end state
loadf = 1; % load initial steady-state
BlFlag = 0; % Blast 1: shot0, 2: cont release
inorgF = 0; % inorganic CaCO3 precip. should be 1 when runninig Permian
% for any other time period should be 0
epspF = 0; %Ca fractionation depends on CO32-;
biopF = 0; %increases intial bio pump
epsSensF = 0; % flag for d44Ca fractionations sensitivuty runs
% loads different steady states in order to get different Ca
% residence time in the ocean; 0-standard, 1 doubles(b) the residence
% time by halving the riverine inputs, 2 - half the residence (c) time
% by doubling riverine input
if(epsSensF==0)
Sfvc = 1.0;
Sfin = 1.0;
epspF=1;
elseif(epsSensF==1)
Sfvc = 1.0;%0.5
Sfin = 1.0;%0.5
epspF=0;
elseif(epsSensF==2)
Sfvc = 1.0;%2.0
Sfin = 1.0;%2.0
epspF=0;
end
fcon = 3; % 1: NADW, 2: NPDW, 3: SO
swcon = 0; % 1: SO -> NP -> SO switch
bath = 2; % 1,2,3 bathymetry; 4-made up bathym for end perm
dsflag = 3; % 1,2,3 dissolution parameter
ftys = 1; % Tethys
fsed = 1; % include sediments
parflag = 0; % write parameter to file
ffflag = 0; % Anthropogenic CO2 (fossil fuel)
fdox = 1; % include dissolved oxygen
TDflag = 0; % Temp sens to doubling CO2
ccdrun = 0; % parameter run: CBl, oxA
oxA = 0.4; % fraction released in deep Atl 0.4
disp(' ');
disp('@==================== RUN Loscar ====================@');
disp(' ');
load_Blast_fcon_bathym_diss_TDflag = ...
sprintf(' %d %d %d %d %d %d',...
loadf,BlFlag,fcon,bath,dsflag,TDflag);
load_Blast_fcon_bathym_diss_TDflag
% Tracer -> DIC ALK PO4 O2 DIC-13 Catm Catm-13 sediments ...
%
% Boxes:
%
% Low, Interm, Deep, High
% Atlantic, Pacific, Indic
%
% AL 1
% IL 2
% PL 3
% AI 4
% II 5
% PI 6
% AD 7
% ID 8
% PD 9
% H 10
Voc = 1.2918235e18; % (m3) volume ocean
Aoc = 3.49e14; % (m2) area ocean
Hav = Voc/Aoc; % (m) average depth
rho = 1.025e3; % kg/m3 (1.025: Toggweiler)
rhos = 2.50e3; % kg/m3 sed. density 2.50
m2kg = 100/1e3; % mol C -> kg CaCO3
m2kgca = 100/1e3; % mol Ca -> kg CaCO3
REDPC = 1/130; % 130 Redfield P:C
REDNC = 15/130; % 15/130 Redfield N:C
REDO2C = 165/130; % 165/130 Redfield O2:C 169
% flags CO2 system
phflag = 0;
k1k2flag = 1;
% Number of oceans
if(ftys)
nOC = 4;
else
nOC = 3;
end;
on3 = ones(1,3);
%--------------------- Ocean Boxes
% A I P
if(ftys)
fA3 = [.15 .14 .52]; %03/31/06 % Area fraction AIP
%fA3 = [.17 .18 .46]; % Area fraction AIP
fH = 0.10; % Area fraction H box
fT = 0.09; % Area fraction Tethys
fA = [fA3 fA3 fA3 fH fT*on3]; %
A = fA*Aoc;
HLI = [100. 900.]; % (m) height L I boxes
DTM = sum(HLI); % (m) depth thermocline
HH = 250.; % (m) depth H box
%HDT = [1000.]; % (m) height Deep Tethys
HDT = [ 200.]; % (m) height Deep Tethys
Vres = Voc-(DTM*(1-fH)+HH*fH+HDT*fT)*Aoc;
%VD = Vres*fA3/(sum(fA3)); % (m3) Vol Deep AIP
VD = Vres*[16.0 16.0 68.]/100;% (m3) Vol Deep AIP
% 04/02/06
HD = VD./(fA3*Aoc); % (m) H Deep AIP
HLID = [HLI HD];
H = [HLID(1)*on3 HLID(2)*on3 HD ...
HH HLI HDT]; % (m) height of boxes
V = A.*H;
else
%--------------------- Ocean Boxes
% A I P
fA3 = [.26 .18 .46 ]; % Area fraction
fA = [fA3 fA3 fA3 .10]; % 0.05
A = fA*Aoc;
HLI = [100. 900.]; % (m) height L I boxes
HD = Hav-sum(HLI);
HLID = [HLI HD];
H = [HLID(1)*on3 HLID(2)*on3 ...
HLID(3)*on3 250.]; % (m) height of boxes
V = fA.*H*Aoc; % (m) Volume of boxes
% volume below H box = A(10)*(Hav-H(10))
% add to deep boxes
V(7:9) = V(7:9) + A(10)*(Hav-H(10))/3;
H = V./A;
end;
% Number of ocean boxes
Nb = length(V);
onV = ones(1,length(V));
% volume of basins
for i=1:3
VO(i) = sum(V([i i+3 i+6]));
end;
VO(4) = V(10);
if (ftys)
VO(5) = sum(V(11:13));
end
% Temperature and feedback
%ntL = 0.4; % 0.3 Low Lat sensitivity
%ntH = 0.5; % 0.4 HighLat sensitivity
TC3 = [20. 10. 2.]; % (degC) temp. of boxes
TCv0 = [TC3(1)*on3 TC3(2)*on3 TC3(3)*on3 2.0];
CA3 = [10.3 10.3 10.3]*1e-3; % mol/kg Ca of boxes
CAv0 = [CA3(1)*on3 CA3(2)*on3 CA3(3)*on3 10.3];
if(ftys)
TC3 = [25. 16. 12.]; % (degC) temp. of boxes
TCv0 = [TC3(1)*on3 TC3(2)*on3 TC3(3)*on3 12.0+0];
TCT = [18. 14. 12.+0]; % 18/25 16/14 12
TCv0 = [TCv0 TCT]+0;
%CA3 = [20 20 20]*1e-3; % mol/kg Ca of boxes
if(inorgF)
if(epsSensF==0)
CA3 = [10 10 10]*1e-3; % mol/kg Ca of boxes
elseif(epsSensF==1)
CA3 = [20 20 20]*1e-3; % mol/kg Ca of boxes
elseif(epsSensF==2)
CA3 = [5 5 5]*1e-3; % mol/kg Ca of boxes
end
else
CA3 = [20 20 20]*1e-3; % mol/kg Ca of boxes
end
CAv0 = [CA3(1)*on3 CA3(2)*on3 CA3(3)*on3 CA3];
end;
TCv = TCv0;
Soc = 34.72; % Sal whole ocean
Sv = onV*Soc; % Salinity vector
CAv= CAv0;
% Pressure vector. Note: H(k+6) all different
for k=1:3
Hv3(k,:)= [H(k)/2 H(k)+H(k+3)/2 H(k)+H(k+3)+H(k+6)/2];
end;
Pv = [Hv3(:,1)' Hv3(:,2)' Hv3(:,3)' H(10)/2]/10;
if(ftys)
k = 11;
HTv = [H(k)/2 H(k)+H(k+1)/2 H(k)+H(k+1)+H(k+2)/2];
Pv = [Pv HTv/10];
end;
% Overturning
TH0 = 20.e6*3600*24*365; % (m3/y) 25 20 Sv conveyor transport
if(ftys)
TH0 = 25.e6*3600*24*365; % (m3/y) 25 20 Sv conveyor transport
end;
TT = 02.e6*3600*24*365; % (m3/y) 03 02 Sv conveyor transport
TH = TH0;
TS = 0.0;
% TH branches
tA = 0.20; % 0.20 upwelled into intermdt Atl 0.27 0.15
tI = 0.20; % 0.20 upwelled into intermdt Ind 0.29 0.30
% mixing A I P TLI TII % 3.5 3.5 8.5
mv0 = [5.5 4.5 6.5 2.5 2.]*1e6; % Sv 5.5 4.5 6.5 2.5 2
% high-deep
mhd0 = [03. 02. 8.0 1.0]*1e6; % Sv 3 2 8 | 4 4 6
if(ftys)
mv0 = [3.5 3.5 7.0 3.2 2.]*1e6; % Sv 5.5 4.5 8.5 3.0 2
mhd0 = [04. 04. 6.0 0.7]*1e6; % Sv 3 2 8 | 4 4 6
end;
mv0 = 3.8*mv0 *365*24*3600; % (m3/y) 3.8 4.0
mhd0 = 1.3*mhd0*365*24*3600; % (m3/y) 1.3
mv = mv0;
mhd = mhd0;
% air-sea CO2/O2
kasv = NaN*onV;
vask = NaN*onV;
if(ftys)
kkv = [1 2 3 10 11];
else
kkv = [1 2 3 10];
end;
xkh = 1.*0.06; % Wally's CO2 exchange coeff.
kasv(kkv) = xkh*A(kkv); % (mol/uatm/y) air sea exch coeff Llat
pv = 3.*365; % (m/day) -> (m/y) piston velocity
vask(kkv) = pv*A(kkv); % m3/y
%============== Biological Pump ============%
%
if(biopF)
EPH = 2*1.8*A(10); % (mol/y) 1.8 1.6 H Export, mol C/m2/y*A = mol/y
fEPL = 0.95; % 0.80 0.9 LL utilization
else
EPH = 1.8*A(10); % (mol/y) 1.8 1.6 H Export, mol C/m2/y*A = mol/y
fEPL = 0.8; % 0.80 0.9 LL utilization
if(Pscenario==13)
fEPL = 0.8;
end
end
PPH = EPH*REDPC;
rrain = 6.1; % 6.1 6.2 6.7 export rain ratio (Corg/CaCO3)
% 5.9(?) 6.1(2,3)
if(ftys)
if(inorgF)
rrain = 6.7*1e20; % 6.7 7 4.2 export rain ratio (Corg/CaCO3)
else
rrain = 6.7;
if(Pscenario==13)
rrain = 8.0;
end
end
end; % 8.0(1,1) 6.3(1,3) 6.6/6.2/6.0?(2,3)
nu = 0.31; % 0.31 water column dissolution
eI = 0.78; % 0.78 0.8 fraction EPL, remineralized in I boxes
if(~fsed)
nu = 0.;
end;
% fraction EPH, remineralized in deep A,I,P boxes
gp = 0.*ones(1,Nb);
gp(7:9) = [.3 .3 .4]; % .3 .3 .4
%gp(7:9) = [1 1 1]/3; % .7 .3 0
%============== silicate weathering: volc degass
pRef = 280.; % uatm, weathering ref 280
pCSi = 280.; % uatm, std-stt atm pCO2 280
if(ftys)
if(~inorgF)
pRef = 0500.*1.; % uatm, weathering ref 500 574 350 750 /2
pCSi = 1000.*1; % uatm, std-stt atm pCO2 560 1000 700 1000 /2 887.8625
else
if(biopF)
pRef = 0157.25*1.; % uatm, weathering ref 370 500 574 350 750 /2
pCSi = 850.*1; % uatm, std-stt atm pCO2 560 1000 700 1000 /2 887.8625
else
pRef = 0157.25*1.; % uatm, weathering ref 500 574 350 750 /2
pCSi = 850.*1; % uatm, std-stt atm pCO2 2000
end
end
if(LTflag)
pCSi = pco2gca(tgc);
end;
end;
if(inorgF)
FVC = Sfvc*3.e12/Aoc; % 5.e12 mol C, degassing /m2/y @280 uatm
else
FVC = 1*5.e12/Aoc; % 5.e12 mol C, degassing /m2/y @280 uatm
end
%nSi = 1.5; % 0.2 0.3
FVC = FVC*(pCSi/pRef)^nSi0; % initial 9.2071
if(LTflag)
FVC = (fmcv(tgc)*1.e12)/Aoc;
FSi0 = ((fbcv(tgc)-fwcv(tgc))*1.e12)/Aoc;
end;
%============= phosphate
if(Pfeed)
Fpw0 = (3.6e10+2.416153e9)/Aoc; % mol P /m2/y modern phosphate weathering flux
else
Fpw0 = 3.6e10/Aoc;
end
oxicf0 = 0.86; % surface water total oxic fraction
anoxf0 = 1-oxicf0; % surface water total anoxic fraction
Fopb0 = 1.5e10/Aoc; % 2.0e10 organic P burial
Ffep0 = 0.6e10/Aoc; % 0.6e10 iron-sorbed P burial
Fcap0 = 1.5e10/Aoc; % 1.5e10 calcium-bound P burial
%============== kerogen oxidation
Focb0 = 1*09.e12/Aoc; % kerogen burial mol C /m2/y 09
Focw0 = 1*09.e12/Aoc; % kerogen weathering
if(ftys)
Focb0 = 04.581e12/Aoc;%-3.95374e+012/Aoc; % mol C /m2/y 05 3.5425
if(Pfeed)
Focw0 = (3.899999889985892e+012+6.282e11)/Aoc;%1*05.e12/Aoc;
else
Focw0 = 1*04.581e12/Aoc;
end
if(LTflag)
Focb0 = (fbgv(tgc)*1.e12)/Aoc; % mol C /m2/y 05
Focw0 =((fwgv(tgc)+fmgv(tgc))*1.e12)/Aoc;
end
end
%============== CaCO3 in-flux ===============%
%
if(inorgF)
FiN = Sfin*5.e12/Aoc; % 12e12 mol C /m2/y riverine flux 1.3
else
FiN = 1.0*12.e12/Aoc; % 12e12 mol C /m2/y riverine flux 1.3
end
% nCC = 0.40; % 0.4 0.3 1.0 0.5
FiN = FiN*(pCSi/pRef)^nC0; % 12.7732
Fpr = 3.*FiN; % mol C /m2/y production 3.6
if(LTflag)
FiN = (fwcv(tgc)*1.e12)/Aoc; % 12.7732
FiN0 = (fwcv(tgc)*1.e12)/Aoc; % 12.7732
end
% rain of 'remainder'
%frrf = 1.5*1.15*0.180; % g/cm2/ky remainder 1.5*1.15
frrf = 0.35; % g/cm2/ky remainder .311
frrf = frrf*1e4/1e3/1e3; % -> kg/ m2/ y
frain = Fpr*m2kg/(Fpr*m2kg+frrf);
frrfca = 0.35; % g/cm2/ky remainder .311
frrfca = frrfca*1e4/1e3/1e3; % -> kg/ m2/ y
%======= Carbon-13
Rst = 0.011; % 13C: R standard (value irrelevant)
epsp = -27.; % -27 fractionation Corg
if(ftys && ~LTflag)
epsp = -33.; % -33 fractionation Corg
end;
if(ftys && LTflag)
epsp = -acv(tgc); % fractionation Corg
end;
d13Cin = +2.0; % d13C of riverine flux 3.0 2.0 2.6
Rin = Rst*(d13Cin/1e3+1);
FiN13 = Rin*FiN; % mol C /m2/y riverine flux
% silicate weathering: volc degass
d13Cvc = -3.0; % d13C -5 +0.3 -0.7 +2.0
if(ftys)
d13Cvc = -5.0; % d13C -5 +0.3 -0.7 +2.0
end;
Rvc = Rst*(d13Cvc/1e3+1);
FVC13 = Rvc*FVC; % mol C /m2/y
% kerogen oxidation
d13Ckg = -22.3; % d13C -22.3 -28.3
Rkg = Rst*(d13Ckg/1e3+1);
% if(LTflag)
Fkg13 = Rkg*Focw0; % mol C /m2/y
% else
% Fkg13 = Rkg*Focb0;
% end
%======= Calcium-44
Rstca = 0.0208; % 13C: R standard (value irrelevant)
epspca = -1.4; % -1.3 fractionation between seawater and carbonate minerals
if(ftys)
if(epspF)
if(biopF)
epspca = -0.9249;
else
epspca = -0.9983; %1.0111, 1.0305
end
else
if(epsSensF==0)
epspca = -0.9249;%-1.4; % -1.3 fractionation between seawater and carbonate minerals
elseif(epsSensF==1)
epspca = -0.9937;%-1.4;
elseif(epsSensF==2)
epspca = -0.5773;%-1.4; -0.9419
end
end
end;
epspca = -1.4;
d44CAin = -0.6; % d44Ca of riverine flux 3.0 2.0 2.6
RinCA = Rstca*(d44CAin/1e3+1);
FiN44ca = RinCA*FiN; % mol C /m2/y riverine flux
% silicate weathering: volc degass
%d44CAvc = -3.0; % d13C -5 +0.3 -0.7 +2.0
if(ftys)
%d44CAvc = -5.0; % d13C -5 +0.3 -0.7 +2.0
end;
% RvcCA = Rstca*(d44CAvc/1e3+1);
% FVC13ca = RvcCA*FVC; % mol C /m2/y
% kerogen oxidation
% d44CAkg = -22.3; % d13C -22.3 -28.3
% Rkgca = Rstca*(d44CAkg/1e3+1);
% Fkg13ca = Rkgca*Focb0; % mol C /m2/y
if(fsed) %========= sediments ============== fsed
% dissolution parameter
%
% ~fc*(1-om)
if (dsflag == 1)
Kd = 1.*365/100 % 1/d -> 1/y
nc = 4.5; % 4.50 calc dissolution order
elseif(dsflag == 2)
% ~fc^0.5*(1-om)
nc = 2.35; % 2.35 calc dissolution order
Kd = 0.12;
elseif(dsflag == 3)
% ~fc^0.5*(cs-c) % Kd defnd in DEQ
nc = 2.40; % 2.40 2.35 calc dissolution order
cst = 100.e-6; %
KS = 20.36e10; % mol/m2/y
end;
%====== shelf/deep rain
fsh = 1.00; % increase shelf rain
fshI = 1;
fshP = 1;
nshT = 1.00; %
if(ftys)
if(inorgF)
fsh = 1; % 8.9
fshI = 1; % 17.2
fshP = 1; % 33.3
nshT = 0.27; % 0.6(1,1) 0.35/0.6(2,3).3
else
fsh = 4.5; % 8.9
fshI = 4.5; % 17.2
fshP = 4.5; % 33.3
nshT = 0.4; % 0.6(1,1) 0.35/0.6(2,3).3
end
end;
%----------------- sediment boxes (bathymetry) ------------%
%
if (bath == 1)
dsv = [ .1 .6 1.5 2.5 3.5 4.5 5.5 6.5]*1000;
asvA = [ 7.0297 5.1729 4.2988 8.5975 19.3421 32.4777 22.3425 0.7388]/100.;
asvI = [ 3.5710 2.6844 3.5792 10.0293 25.2598 36.6442 16.9915 1.2407]/100.;
asvP = [ 1.6358 2.5901 3.2590 6.8744 21.8550 35.0822 26.9567 1.7468]/100.;
if(ftys)
asvT = [16.3934 16.3934 16.3934 20.4918 20.4918 8.1967 1.6393 0.0001]/100.;
end;
elseif(bath == 2)
if(ftys) % 03/31/06, 2x2\deg Bice
dsv = [.1 .6 1 1.5 2. 2.5 3. 3.5 4. 4.5 5. 5.5 6.5]*1000; % depth (m)
% area fraction A I P
asvA = [1.1407 10.0216 8.7160 6.3724 4.7915 3.4973 12.2935 ...
10.7438 8.4983 11.4134 11.0948 11.4167 1e-6]/100;
asvI = [0.2501 5.5345 5.5145 8.9550 4.6283 6.5361 11.7221 ...
12.6050 14.6295 13.8384 8.0807 7.7057 1e-6]/100;
asvP = [0.1673 2.8333 3.0599 2.5389 1.4218 5.0153 9.8023 ...
14.0117 10.1975 20.0019 13.7155 17.2346 1e-6]/100;
asvT = [7.0534 46.5363 22.4068 7.1501 2.4261 3.9946 2.7063 ...
0.8532 0.9814 3.0802 2.3583 0.4533 1e-6]/100;
else %#! old 03/22/06
dsv = [.1 .6 1 1.5 2. 2.5 3. 3.5 4. 4.5 5. 5.5 6.5]*1000; % depth (m)
% area fraction A I P
asvA = [7.0297 5.1729 1.9106 2.3882 4.2988 4.2988 9.6711 ...
9.6711 16.2389 16.2389 11.1712 11.1712 0.7388]/100;
asvI = [3.5710 2.6844 1.5907 1.9884 5.0146 5.0146 12.6299 ...
12.6299 18.3221 18.3221 8.4957 8.4957 1.2407]/100;
asvP = [1.6358 2.5901 1.4484 1.8105 3.4372 3.4372 10.9275 ...
10.9275 17.5411 17.5411 13.4784 13.4784 1.7468]/100;
%if(ftys)
%asvT = [16.3934 16.3934 7.2860 9.1075 10.2459 10.2459 10.2459 ...
% 10.2459 4.0984 4.0984 0.8197 0.8197 0.0001]/100;
%end;
end; % ftys old/new
elseif(bath == 3)
dsv = [.1 .6 1 1.5 2 2.5 2.75 3 3.25 3.5 3.75 4 ...
4.25 4.5 4.75 5. 5.25 5.5 6.0 6.5]*1000; % depth (m)
% area fraction A I P
asvA = [07.0297 5.1729 1.9106 2.3882 4.2988 4.2988 4.8355 4.8355 4.8355 4.8355 8.1194 8.1194 8.1194 8.1194 5.5856 5.5856 5.5856 5.5856 0.3694 0.3694]/100;
asvI = [03.5710 2.6844 1.5907 1.9884 5.0146 5.0146 6.3149 6.3149 6.3149 6.3149 9.1610 9.1610 9.1610 9.1610 4.2479 4.2479 4.2479 4.2479 0.6204 0.6204]/100;
asvP = [01.6358 2.5901 1.4484 1.8105 3.4372 3.4372 5.4638 5.4638 5.4638 5.4638 8.7705 8.7705 8.7705 8.7705 6.7392 6.7392 6.7392 6.7392 0.8734 0.8734]/100;
if(ftys)
asvT = [16.3934 16.3934 7.2860 9.1075 10.2459 10.2459 5.1229 5.1229 5.1229 5.1229 2.0492 2.0492 2.0492 2.0492 0.4098 0.4098 0.4098 0.4098 0.0001 0.0001]/100;
end;
elseif(bath == 4)
% 03/31/06, 2x2\deg Bice
dsv = [.1 .6 1 1.5 2. 2.5 3. 3.5 4. 4.5 5. 5.5 6.5]*1000; % depth (m)
% area fraction A I P
asvA = [1.1407*2 10.0216*2 8.7160 6.3724 4.7915 3.4973 12.2935 ...
10.7438 8.4983 11.4134 11.0948-1.1407 11.4167-10.0216 1e-6]/100;
asvI = [0.2501*2 5.5345*2 5.5145 8.9550 4.6283 6.5361 11.7221 ...
12.6050 14.6295 13.8384 8.0807-0.2501 7.7057-5.5345 1e-6]/100;
asvP = [0.1673*2 2.8333*2 3.0599 2.5389 1.4218 5.0153 9.8023 ...
14.0117 10.1975 20.0019 13.7155-0.1673 17.2346-2.8333 1e-6]/100;
asvT = [7.0534*3 46.5363 22.4068-2*7.0534 7.1501 2.4261 3.9946 2.7063 ...
0.8532 0.9814 3.0802 2.3583 0.4533 1e-6]/100;
end;
% Number of sediment boxes
Ns = length(dsv);
onNs = ones(1,Ns);
% kLID: assign sediment to ocean boxes
% Low, Interm, or Deep
kl = find(dsv <= HLI(1));
ki = find(dsv > HLI(1) & dsv <= (HLI(1)+HLI(2)));
kd = find( dsv > (HLI(1)+HLI(2)));
klid(kl) = 01; klid(ki) = 04; klid(kd) = 07;
if(ftys)
kliT(kl) = 11; kliT(ki) = 12; kliT(kd) = 13;
end
nli = [length(kl) (length(kl)+length(ki))];
zv = [000.:10:6000.]; % z, continuous (1 or 10 m)
lzv = length(zv);
% calculate calcite saturation at depth of sediment boxes
satflg = 2; % 1/2, 1:Wally 2:Millero
zsatv = dsv;
if (satflg == 1)
as = 0.189/1.e3;
zs0 = 3.82e3; % 3.82 m
co3s0 = 88.7e-6; % 88.7 mol/kgn
co3satv = co3s0*exp(as*(zsatv-zs0));
elseif(satflg == 2)
% 2. Millero
Cam = 10.3e-3; % 10.3 (mol/kg) modern
Mgm = 53.0e-3; % 53.0 (mol/kg)
Ca = Cam;
Mg = Mgm;
if(ftys)
if(inorgF)
if(epsSensF==0)
Ca = 10.00e-3; % Modern 10.3 PETM 20.0
elseif(epsSensF==1)
Ca = 20.00e-3; % Modern 10.3 PETM 20.0
elseif(epsSensF==2)
Ca = 5.0e-3; % Modern 10.3 PETM 20.0
end
else
Ca = 20.00e-3; % Modern 10.3 PETM 20.0
end
Mg = 30.0e-3; % 53.0 30.0
end;
if(CAvflag == 0)
Tdv = TCv(klid);
Sdv = Sv(klid);
for k=1:Ns
[kspc(k),x] = ...
kspfun(Tdv(k),Sdv(k),zsatv(k)/10.,Ca,Mg);
end;
co3satv = kspc/Ca;
elseif(CAvflag == 1)
Tdv = TCv(klid);
Sdv = Sv(klid);
for k=1:Ns
[kspc(k),x] = ...
kspfun(Tdv(k),Sdv(k),zsatv(k)/10.,Ca,Mg);
end;
co3satv = kspc/Ca;
else
% Temp for co3sat corrected 07/13/06
Tdv = TCv(klid);
Sdv = Sv(klid);
Cadv = CAv(klid);
CALCv=Ca;
clear kspc;
for k=1:Ns
[kspc(k),x] = ...
kspfun(Tdv(k),Sdv(k),zsatv(k)/10.,CALCv,Mg);
end;
co3satv = kspc./(CALCv*1e-3);
end
end;% satflag
kspc
%-------------- Porosity --------------------%
%phic = 0.78; % porosity 0.75 0.78
if(phic) % do NOT use exist! phic does exist (global!)
phiiA = ones(1,Ns)*phic;
phiiI = ones(1,Ns)*phic;
phiiP = ones(1,Ns)*phic;
phiiAca = ones(1,Ns)*phic;
phiiIca = ones(1,Ns)*phic;
phiiPca = ones(1,Ns)*phic;
if(ftys)
phiiT = ones(1,Ns)*phic;
phiiTca = ones(1,Ns)*phic;
end;
else % phic
phi0 = 0.85; % porosity max 0.85 0.88
gam = 0.23; % 0.23 0.28
phi1 = phi0-gam; % porosity min
end;
hs = 0.08; % (m ) bioturbated layer 0.08 0.1
VsvA = asvA*A(1)*hs; % (m3) Volume sediment boxes A
VsvI = asvI*A(2)*hs; % (m3) Volume sediment boxes I
VsvP = asvP*A(3)*hs; % (m3) Volume sediment boxes P
VsA = sum(VsvA);
VsI = sum(VsvI);
VsP = sum(VsvP);
if(ftys)
VsvT = asvT*A(11)*hs;% (m3) Volume sediment boxes T
VsT = sum(VsvT);
end;
% set initial calcite fraction
fc0A = 0.46*ones(1,Ns);
fc0I = 0.46*ones(1,Ns);
fc0P = 0.46*ones(1,Ns);
fca0A = 0.46*ones(1,Ns);
fca0I = 0.46*ones(1,Ns);
fca0P = 0.46*ones(1,Ns);
if(ftys)
fc0T = 0.46*ones(1,Ns);
fca0T = 0.46*ones(1,Ns);
end;
% calc initial phi
if(isempty(phic)) % phi = phi(fc)
FF = (phi1-phi0)/(1-phi1);
phiiA = (phi0+FF*fc0A)./(1+FF*fc0A);
phiiI = (phi0+FF*fc0I)./(1+FF*fc0I);
phiiP = (phi0+FF*fc0P)./(1+FF*fc0P);
FFca = (phi1-phi0)/(1-phi1);
phiiAca = (phi0+FFca*fca0A)./(1+FFca*fca0A);
phiiIca = (phi0+FFca*fca0I)./(1+FFca*fca0I);
phiiPca = (phi0+FFca*fca0P)./(1+FFca*fca0P);
if(ftys)
phiiT = (phi0+FF*fc0T)./(1+FF*fc0T);
phiiTca = (phi0+FFca*fca0T)./(1+FFca*fca0T);
end;
end;
% calc initial calcite mass
mc0vA = (fc0A.*rhos.*(1-phiiA));
mc0vI = (fc0I.*rhos.*(1-phiiI));
mc0vP = (fc0P.*rhos.*(1-phiiP));
mca0vA = (fca0A.*rhos.*(1-phiiAca));
mca0vI = (fca0I.*rhos.*(1-phiiIca));
mca0vP = (fca0P.*rhos.*(1-phiiPca));
if(ftys)
mc0vT = (fc0T.*rhos.*(1-phiiT));
mca0vT = (fca0T.*rhos.*(1-phiiTca));
end;
%====== Carbon-13
m13c0vA = Rin*mc0vA; % -> kg CaCO3/m3 Rin
m13c0vI = Rin*mc0vI; % -> kg CaCO3/m3 Rin
m13c0vP = Rin*mc0vP; % -> kg CaCO3/m3 Rin
f13c0A = fc0A.*m13c0vA./mc0vA;
f13c0I = fc0I.*m13c0vI./mc0vI;
f13c0P = fc0P.*m13c0vP./mc0vP;
if(ftys)
m13c0vT = Rin*mc0vT; % -> kg CaCO3/m3 Rin
f13c0T = fc0T.*m13c0vT./mc0vT;
end;
%====== Calcium-44
m44ca0vA = RinCA*mca0vA; % -> kg CaCO3/m3 Rin
m44ca0vI = RinCA*mca0vI; % -> kg CaCO3/m3 Rin
m44ca0vP = RinCA*mca0vP; % -> kg CaCO3/m3 Rin
f44ca0A = fca0A.*m44ca0vA./mca0vA;
f44ca0I = fca0I.*m44ca0vI./mca0vI;
f44ca0P = fca0P.*m44ca0vP./mca0vP;
if(ftys)
m44ca0vT = RinCA*mca0vT; % -> kg CaCO3/m3 Rin
f44ca0T = fca0T.*m44ca0vT./mca0vT;