-
Notifications
You must be signed in to change notification settings - Fork 0
/
revscavMn.m
76 lines (60 loc) · 2.43 KB
/
revscavMn.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
function [A,b,output] = revscavMn(A,b,do);
%REVSCAVPOC scavenges your element downwards by reversible scavenging
% This function adds reversible scavenging of your element. It is based
% on the assumption that your element scavenges like Th from the Hulten model.
fprintf('%s','revscavMn...')
% load the grid and POC data (units of mmole m-3)
load ([do.highestpath '/data/ao.mat'])
load pMn
% unpack the scavenging rate constant
K = do.revscavMn.K;
sedremin = do.revscavMn.sedremin;
% get the Mn concentrations from each grid cell
pMn = pMn(ao.iocn);
% multiply by 10^9 to get nM
pMn = pMn*10^9;
% the amount of element E which sinks out depends on the equilibrium
% scavenging constant, multiplied by the POC concentration, times the
% sinking rate divided by the height of the grid cell
sinkout = K*(pMn./ao.Height);
% find the equation position (positions in the A matrix) of the grid cells
% which lie below each cell
EQNPOSBELOW = cat(3,ao.EQNPOS(:,:,2:24),zeros(91,180,1));
% define the equation positions that particles are sinking from, as well as
% the volumes and heights of those grid cells
frompos = ao.EQNPOS(EQNPOSBELOW~=0);
fromvol = ao.Vol(frompos);
fromheight = ao.Height(frompos);
% define the equation positions that particles are sinking to, as well as
% the volumes and heights of those grid cells
topos = EQNPOSBELOW(EQNPOSBELOW~=0);
tovol = ao.Vol(topos);
% create the sinkout A matrix, and fill in the diagonal with the magnitude
% of the sinking flux out
sinkoutA = speye(ao.nocn,ao.nocn);
sinkoutA(sinkoutA==1)=sinkout;
% calculate the amount of element transferred into each grid cell by
% sinking with K, the sinking rate devided by the grid cell height from
% which sinking occurs, and a correction for volume
sinkin = sinkout(frompos).*(fromvol./tovol);
% create the A matrix for sinking in
sinkinA = sparse(topos,frompos,sinkin,ao.nocn,ao.nocn);
% find the equation positions of cells which lie on the bottom, and the
% amount of reminerlization there is equal to the amount which sinks out of
% that grid cell
btmeqnpos = ao.EQNPOS(ao.ibtm);
sedreminA = sparse(btmeqnpos,btmeqnpos,sinkout(btmeqnpos),ao.nocn,ao.nocn);
% add the A matrix with the sinking matrices
A = A - sinkoutA + sinkinA;
% add sedimentary remineralization if switched on
if sedremin
A = A + sedreminA;
end
% package outputs
output.K=K;
output.sinkoutA=sinkoutA;
output.sinkinA=sinkinA;
if sedremin
output.sedreminA=sedreminA;
end
output.citations=cell(1,1);