-
Notifications
You must be signed in to change notification settings - Fork 0
/
buildPFD_Chl.m
77 lines (76 loc) · 2.5 KB
/
buildPFD_Chl.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
function [PFdiv,dPFDdk1,dPFDdb1] = buildPFD_Chl(M3d,p,grd);
% this code is used to build Particle Flux Diverngence (PFD)
%
% ________________________
% / A /| POP sinking
% top/_______________________/ | |
% | | | | | -w
% | / V | / |
% bot|/______________________|/ V
%
% PFD = (A*w(top)*POP(top)-A*w(bot)*POP(bot))/V;
% add an exra layer of zeros at the bottom to ensure there is no
% flux in or out of the domain when using periodic shift operators
% for finite differences and averaging
[ny,nx,nz] = size(M3d);
M3D = zeros(ny,nx,nz+1);
M3D(:,:,1:end-1) = M3d;
% add the zw coordinate at the top of the extra layer
ZW3d = grd.ZW3d;
ZW3d = ZW3d(:,:,[1:end,end]);
ZW3d(:,:,end) = grd.ZW3d(:,:,end)+grd.dzt(end);
% areas of the top of the grid box
dAt = (grd.DXT3d.*grd.DYT3d).*M3d;
% volume of the grid boxes
dVt = (grd.DXT3d.*grd.DYT3d.*grd.DZT3d).*M3d;
%
n = nx*ny*(nz+1);
I0 = speye(n);
i0 = zeros(ny,nx,nz+1);
i0(:) = 1:n;
% periodic shifts OK because M3D has a layer of zeros on the bottom
iu = i0(:,:,[nz+1,1:nz]); %use a periodic upward shift
ib = i0(:,:,[2:nz+1,1]); % use a periodic downward shift
IU = I0(iu,:);
IB = I0(ib,:);
% keep only wet boxes
iocn = find(M3D(:));
I0 = I0(iocn,:); I0 = I0(:,iocn);
IU = IU(:,iocn); IU = IU(iocn,:);
IB = IB(:,iocn); IB = IB(iocn,:);
% (averages POP onto the top of the grid boxes)
%AVG = d0((I0+IU)*M3D(iocn))\(I0+IU);
% (compute the divergence in the center of the grid boxes)
DIV = d0(dVt(iocn))\(I0-IB)*d0(dAt(iocn));
% (compute the flux at the top of the grid cells)
% mimics a Martin curve flux attenuation profile
%(see Kriest and Oschelies 2008 in Biogeosciences)
%r = p.kappa_p;
r = p.d1+p.k1;
b = p.bs;
a = r/b;
% particle sinking velocity at the top of the grid cells.
MSK = M3D.*M3D(:,:,[nz+1,1:nz]);
M = MSK.*ZW3d;
w = -a*M;
dadb = -r/(b^2);
dadk1 = 1/b;
dadd1 = 1/b;
dwdb = -dadb*M;
dwdk1 = -dadk1*M;
dwdd1 = -dadd1*M;
d2adb2 = 2*r/(b^3);
d2wdb2 = -d2adb2*M;
%FLUX = d0(w(iocn))*AVG;
FLUX = d0(w(iocn))*IU;
dFLUXdb = d0(dwdb(iocn))*IU;
dFLUXdk1 = d0(dwdk1(iocn))*IU;
dFLUXdd1 = d0(dwdd1(iocn))*IU;
d2FLUXdb2 = d0(d2wdb2(iocn))*IU;
% particle flux divergence operator
PFdiv = DIV*FLUX;
dPFDdb = DIV*dFLUXdb;
dPFDdk1 = DIV*dFLUXdk1;
dPFDdb1 = DIV*dFLUXdd1;
d2PFDdb2 = DIV*d2FLUXdb2;
end