| 1 |
function [psi,mskG,ylat] = calcEulerPsiCube(varargin); |
| 2 |
|
| 3 |
% [psi,mskG,ylat] = calcEulerPsiCube(d,g,flu,rstar,blkFile,[mask]); |
| 4 |
% |
| 5 |
% Input arguments: |
| 6 |
% d [structure] Velocity field (Mass-weighted if rstar=1): |
| 7 |
% UVELMASS,VVELMASS for rstar=1 |
| 8 |
% UVEL,VVEL for rstar=0 |
| 9 |
% g [structure] drF,dxG,dyG,HFacW,HFacS |
| 10 |
% flu [string] 'O' or 'A' for ocean or atmosphere |
| 11 |
% rstar [integer] 1 or 0 if you are using r* coordinates or not |
| 12 |
% blkFile [string] Broken line file ('isoLat_cs32_59.mat') |
| 13 |
% |
| 14 |
% Optional inputs: |
| 15 |
% mask [structure] Optional: Mask field for computation per basin, it assumes that |
| 16 |
% maskW and maskS are provided in a structure |
| 17 |
% |
| 18 |
% Output fields: |
| 19 |
% psi Overturning |
| 20 |
% mskG Land mask |
| 21 |
% ylat Latitude coordinate of psi |
| 22 |
% |
| 23 |
% Description: |
| 24 |
% Calculates overturning streamfunction (psi). For the atmosphere, data |
| 25 |
% is must be in p-coordinates and the output is the mass transport psi |
| 26 |
% [10^9 kg/s]. For the ocean, data should be in z-coordinates and the |
| 27 |
% output is the volume transport psi [10^6 m^3/s = Sv]. If rstar |
| 28 |
% is on (=1), UVELMASS and VVELMASS must be used. If off (rstar=0), |
| 29 |
% the hfacw*.UVEL and hfacs*.VVEL are used (the multiplication being |
| 30 |
% done inside the function). |
| 31 |
% |
| 32 |
% 'psi' is tabulated on broken lines at the interface between cells in |
| 33 |
% the vertical. 'mskG' is for the area between broken lines and between |
| 34 |
% the cell interfaces in the vertical. |
| 35 |
% |
| 36 |
|
| 37 |
% Defaults that can be overriden. |
| 38 |
grav = 9.81; |
| 39 |
masking=0; |
| 40 |
nBas=0; |
| 41 |
|
| 42 |
% Read input parameters. |
| 43 |
d = varargin{1}; |
| 44 |
g = varargin{2}; |
| 45 |
flu = varargin{3}; |
| 46 |
rstar = varargin{4}; |
| 47 |
blkFile = varargin{5}; |
| 48 |
if length(varargin) == 6 |
| 49 |
mask = varargin{6}; |
| 50 |
masking = 1; |
| 51 |
end |
| 52 |
|
| 53 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 54 |
% Prepare / reform incoming data % |
| 55 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 56 |
|
| 57 |
nc = size(g.XC,2); |
| 58 |
nr = length(g.drF); |
| 59 |
|
| 60 |
delM = g.drF; |
| 61 |
dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]); |
| 62 |
dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]); |
| 63 |
if rstar |
| 64 |
nt = size(d.UVELMASS,4); |
| 65 |
hu = reshape(d.UVELMASS(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]); |
| 66 |
hv = reshape(d.VVELMASS(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]); |
| 67 |
else |
| 68 |
nt = size(d.UVEL,4); |
| 69 |
hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 70 |
hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 71 |
hu = reshape(d.UVEL(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]); |
| 72 |
hv = reshape(d.VVEL(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]); |
| 73 |
for it = 1:nt |
| 74 |
hu(:,:,it) = hw.*hu(:,:,it); |
| 75 |
hv(:,:,it) = hs.*hv(:,:,it); |
| 76 |
end |
| 77 |
end |
| 78 |
|
| 79 |
mskWloc = ones(6*nc*nc,1); |
| 80 |
mskSloc = ones(6*nc*nc,1); |
| 81 |
|
| 82 |
if masking == 1 |
| 83 |
mskWloc=reshape(mask.maskW(:,:,1),6*nc*nc,1); |
| 84 |
mskSloc=reshape(mask.maskS(:,:,1),6*nc*nc,1); |
| 85 |
%hu = repmat(reshape(mask.maskW,6*nc*nc,1),[1 nr nt]) .* hu; |
| 86 |
%hv = repmat(reshape(mask.maskS,6*nc*nc,1),[1 nr nt]) .* hv; |
| 87 |
end |
| 88 |
|
| 89 |
% Load broken information. |
| 90 |
load(blkFile); |
| 91 |
ydim = length(bkl_Ylat); |
| 92 |
ylat = [-90,bkl_Ylat,90]; |
| 93 |
|
| 94 |
% Prepare arrays. |
| 95 |
psi = zeros(ydim+2,nr+1,1+nBas,nt); |
| 96 |
mskZ = zeros(ydim+2,nr+1,1+nBas); % Mask of psi |
| 97 |
mskV = zeros(ydim+2,nr,1+nBas); % Mask of the merid. transport |
| 98 |
mskG = zeros(ydim+1,nr,1+nBas); % Mask of the ground |
| 99 |
|
| 100 |
% The variable "bkl_Flg" is -1/1 if edge (on a given broken) has a u point |
| 101 |
% and -2/2 if it has a v point. Positive/negative values contribute |
| 102 |
% positively/negatively to northward heat transport (this depends on the |
| 103 |
% oreientation of the cell). A zero value indicates an end of edges that |
| 104 |
% contribute to a broken line. The u and v information is parced into two |
| 105 |
% seperate fields, ufac and vfac (-2/2 are reduced to -1/1 for vfac). |
| 106 |
ufac = zeros([size(bkl_Flg),1+nBas]); |
| 107 |
vfac = zeros([size(bkl_Flg),1+nBas]); |
| 108 |
ufac(:,:,1) = rem(bkl_Flg,2); |
| 109 |
vfac(:,:,1) = fix(bkl_Flg/2); |
| 110 |
|
| 111 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 112 |
% Compute mass/volume stream function % |
| 113 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 114 |
|
| 115 |
% Compute volume transport through broken lines a hence psi. ut/vt is the |
| 116 |
% velocity times the edge length it is passing through. The sum of this |
| 117 |
% quantity along a broken line (vz) times the cell height is the volume |
| 118 |
% transport through broken line at one layer (delM(k)*vz). psi is then |
| 119 |
% the value of the volume transport through the level above subtracted |
| 120 |
% from the value of psi above. |
| 121 |
for it = 1:nt |
| 122 |
for k = nr:-1:1 |
| 123 |
ut = dyg.*hu(:,k,it).*mskWloc; |
| 124 |
vt = dxg.*hv(:,k,it).*mskSloc; |
| 125 |
for jl = 1:ydim |
| 126 |
ie = bkl_Npts(jl); |
| 127 |
for b = 1:1+nBas |
| 128 |
vz = sum( ufac(1:ie,jl,b).*ut(bkl_IJuv(1:ie,jl)) ... |
| 129 |
+ vfac(1:ie,jl,b).*vt(bkl_IJuv(1:ie,jl)) ); |
| 130 |
psi(jl+1,k,b,it) = psi(jl+1,k+1,b,it) - delM(k)*vz; |
| 131 |
end |
| 132 |
end |
| 133 |
end |
| 134 |
end |
| 135 |
|
| 136 |
psi = squeeze(psi); |
| 137 |
|
| 138 |
%% For Ocean, result in Sv (10^6 m3/s) |
| 139 |
%% For Atmos, results in 10^9 kg/s |
| 140 |
if isequal(flu,'O'), psi = 1e-6*squeeze(psi); end |
| 141 |
if isequal(flu,'A'), psi =-1e-9/grav*squeeze(psi); end |
| 142 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph