hydrotrans.m
by
Webmaster Legos
—
last modified
Jul 05, 2012 02:18 PM
hydrotrans.m
—
Objective-C source code,
15 kB (15,853 bytes)
File contents
%function
% KEY:
% USAGE :
%
%
%
%
% DESCRIPTION :
%
%
% INPUT: IPdir boxname modelid (parameters)
% <boxname>_equn.mat : contains equation matrix and parameters
% <boxname>_bhat.mat : contains solution bhat, P from inversion
%
%
% OUTPUT: (section isb, layer ilay, property iprop)
% Tr{isb}(ilay,iprop) : relative transport (including Ekman)
% Ta{isb}(ilay,iprop) : absolute transport from inversion
% dTa{isb}(ilay,iprop): uncertainty (1std. dev)
%
% Wtrans(ilayint,iprop): vertical transport through interface ilayint
% ilayint=1:Nlay-2
% Dwtrans(ilayint,iprop): uncertainty
%
% Kztrans{ibox}% diffusive transport
% Dkztrans{ibox}% diffusive transport uncertainty
%
% Resr(ilay,iprop): residuals in the box before inversion
% Resa(ilay,iprop): residuals after inversion
% Dres(ilay,iprop): uncertainty on residuals
%
% AUTHOR : A.Ganachaud (ganacho@gulf.mit.edu) , Jul 97
%
% SIDE EFFECTS :
%
% SEE ALSO : mkequats
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALLER:
% CALLEE:
clear Wtrans Dwtrans boxsec Kztrans Dkztrans
clear Tbi Tr cTr Ta cTa dTa tekk Resr Resa Dres Dresr
if ~exist('p_res_units')
p_res_units=0; %convert residual units in mol/m2/yr
end
if ~exist('p_calculate_Ekman')
p_calculate_Ekman=0;
end
if ~exist('bhat')
str=['load ' IPdir 'bhat_' Modelid '_' invid '.mat'];
disp(str)
eval(str)
if 0
IPdir='/data1/ganacho/Boxmod/Natl/';
boxname='natl_XVII';
modelid='Alison';
IPdir='/data1/ganacho/Boxmod/Natl/';
Modelid='NewNatl';
invid='GM1.0';
end
else
end
str=[ 'load ' IPdir Modelid '_equn.mat '];
disp(str)
eval(str)
if ~exist('boxsec')
boxsec=ones(1,length(gsecs.name));
gisb=1:length(gsecs.name);
%boxsec=[1 1 1 2]; %Which section belogs to which box
%gisb =[1 2 3 2]; %relative indice in the box
end
%(For Ekman transport calculation)
Nconseq=size(boxi.liavg,2);
Nlay=boxi.nlay;
Nsec=length(gsecs.name);
if any(strcmp(fieldnames(pm),'gifwcol')) %contains multiple boxes
gifwcol=pm.gifwcol;
gilwcol=pm.gilwcol;
multbox=1;
if exist('p_compw')& ~p_compw
error('indices will be wrong if multibox and different layers')
end
else
gifwcol=pm.ifwcol;
gilwcol=pm.ilwcol;
pm.ibfrow=1;
multbox=0;
end
if any(strcmp(fieldnames(pm),'gifKzcol')) %contains multiple boxes
gifkzcol=pm.gifKzcol;
gilkzcol=pm.gilKzcol;
elseif any(strcmp(fieldnames(pm),'ifKzcol'))
gifkzcol=pm.ifKzcol;
gilkzcol=pm.ilKzcol;
end
Nbox=length(gifwcol);
for iprop=1:Nconseq
disp(sprintf('CALCULATING %s TRANSPORTS ...',upper(propnm{iprop})))
%COMPUTE NET TRANSPORT FOR EACH SECTION/LAYER
for isa=1:Nsec
girow=pm.ibfrow(boxsec(isa))-1+[(iprop-1)*Nlay+1:iprop*Nlay];
gicol=pm.gifcol(isa):pm.gilcol(isa);
if exist('p_combinesec')&~isempty(p_combinesec{isa})&...
any(p_combinesec{isa})
for isa2=p_combinesec{isa}
if isa2==isa
error('self combination : p_combinesec(isa)==isa !')
end
disp(['COMBINING ' gsecs.name{isa} ' WITH ' ...
gsecs.name{isa2}])
gicol=[gicol,pm.gifcol(isa2):pm.gilcol(isa2)];
if multbox
gsecs.Ekt{boxsec(isa)}(:,gisb(isa),iprop)=...
gsecs.Ekt{boxsec(isa)}(:,gisb(isa),iprop)+...
gsecs.Ekt{boxsec(isa2)}(:,gisb(isa2),...
iprop);
else
gsecs.Ekt(:,gisb(isa),iprop)=gsecs.Ekt(:,gisb(isa),iprop)+...
gsecs.Ekt(:,gisb(isa2),iprop);
end
end %for isa2=p_combinesec{isa}
if isempty(findstr(gsecs.name{isa},'+'))
gsecs.name{isa}=[gsecs.name{isa} '+'];
end
end
if multbox
tekk=gsecs.Ekt{boxsec(isa)}(:,gisb(isa),iprop);
inboxdirr=gsecs.inboxdir{boxsec(isa)}(gisb(isa));
else
tekk=gsecs.Ekt(:,gisb(isa),iprop);
inboxdirr=gsecs.inboxdir(gisb(isa));
end
%Re-CREATE ORIGINAL A MATRIX FOR NET TRANSPORT COMPUTATION IF ANOMALY
%WERE USED IN CONSERVATION EQUATIONS.
if any(strcmp(fieldnames(boxi.conseq),'anom'))&...
length(boxi.conseq.anom)>=iprop &...
boxi.conseq.anom(iprop)~=0 & iprop~=1
if iprop ==2
disp('NET TRANSPORTS ARE ABSOLUTE (NOT ANOMALEOUS)')
end
C0=boxi.lavg(:,iprop)*1e9/1000/lscale(iprop);
Amatnet=Amat(girow,:)+C0*ones(1,size(Amat,2)).*Amat(1:Nlay,:);
Gunsgnnet=Gunsgn(girow,:)+C0*ones(1,size(Amat,2)).*Gunsgn(1:Nlay,:);
if exist('p_t2bsalianom') & p_t2bsalianom & ...
strcmp(propnm{iprop},'sali')
disp('KEEP TOP-TO-BOTTOM ANOMALY FOR SALT')
C0=ones(Nlay,1)*boxi.lavg(Nlay,iprop)*1e9/1000/lscale(iprop);
Amatnet=Amatnet-...
C0*ones(1,size(Amat,2)).*Amat(1:Nlay,:);
Gunsgnnet=Gunsgnnet-...
C0*ones(1,size(Amat,2)).*Gunsgn(1:Nlay,:);
tekk=tekk-C0.*gsecs.Ekt(:,gisb(isa),1);
if any(strcmp(fieldnames(gsecs),'wbc'))&...
(length(gsecs.wbc)>=isa) &...
~isempty(gsecs.wbc{isa})
Wbccor([1,boxi.nlay])=-C0([1,boxi.nlay])*gsecs.wbc{isa}(1);
%!! C0 is used again below, for baroclinic contrib.
Wbccor=Wbccor(:);
else
Wbccor=0;
end
end
else
Amatnet=Amat(girow,:);
Gunsgnnet=Gunsgn(girow,:);
Wbccor=0;
end
if exist('p_save_Amat')&p_save_Amat %save Amat rows for
Amatnets{isa}(iprop,:)=-inboxdirr*...
full(Amatnet(length(girow),gicol));
end
%ADD IN WESTERN BOUNDARY CURRENT IF NECESSARY IN THE FIRST LAYER
%ALL IN FIRST LAYER
if any(strcmp(fieldnames(gsecs),'wbc'))&(length(gsecs.wbc)>=isa) &...
~isempty(gsecs.wbc{isa})
disp(['adding bcurrent across section ' gsecs.name{isa}])
Twbc([1,boxi.nlay])=gsecs.wbc{isa}(iprop);
Twbc=Twbc(:)+Wbccor;
%Wbc correction if salinity anomaly transport
else
Twbc=0;
end
if ~exist('p_pltstream')|p_pltstream==0
Tr{isa}(:,iprop)=sum(Gunsgnnet(:,gicol),2)+tekk+...
-inboxdirr*Amatnet(:,gicol)*Binit(gicol)+Twbc;
Ta{isa}(:,iprop)=sum(Gunsgnnet(:,gicol),2)+tekk+...
-inboxdirr*Amatnet(:,gicol)*bhat(gicol)+Twbc;
if exist('p_savealltrans')&p_savealltrans
h_savealltrans
end
dTa{isa}(:,iprop)=...
full(sqrt(diag(Amatnet(:,gicol)*P(gicol,gicol)*Amatnet(:,gicol)')));
%Cumulative transports from West or North
disp('Does not include Ekman transport for Cumulative transport')
np=gsecs.npair(isa)-any(strcmp(fieldnames(pm),'giEkcol'));
cTr{isa,iprop}=cumsum(Gunsgnnet(:,gicol(1:np)),2)+...
-inboxdirr*cumsum(Amatnet(:,gicol(1:np)).*(ones(Nlay,1)*...
Binit(gicol(1:np))'),2);
cTa{isa,iprop}=cumsum(Gunsgnnet(:,gicol(1:np)),2)+...
-inboxdirr*cumsum(Amatnet(:,gicol(1:np)).*(ones(Nlay,1)*...
bhat(gicol(1:np))'),2);
%Get barotropic and baroclinic transports versus horiz.
%Barotropic (net T X thetabar) is last layer; BI is sum
%1:Nlay-1
if iprop > 1
ddivd=Lays{isa,1}.lavgprop; ddivd(~ddivd)=Inf;
vzonal=Ta{isa}(:,1)./ddivd';
vertlayratio=Lays{isa,iprop}.lverarea(1:Nlay-1)'/...
Lays{isa,iprop}.lverarea(Nlay);
vzonal=vzonal-[vzonal(Nlay)*vertlayratio;0];
rhoTzonal=Lays{isa,iprop}.lavgprop';
-[Lays{isa,iprop}.lavgprop(Nlay)*ones(Nlay-1,1);0];
if iprop==2
Tbi{isa}(:,iprop)=rhoTzonal.*vzonal*lscale(1)/lscale(2);
else
Tbi{isa}(:,iprop)=rhoTzonal.*vzonal; %kmol/s
end
if exist('p_t2bsalianom') & p_t2bsalianom & ...
strcmp(propnm{iprop},'sali')
%in that case removes barotropic component as it is
% zero by definition. BI is OK.
Tbi{isa}(Nlay,iprop)=0;
end
end
else
%Aug 98:
disp('Get cumulative transports from bottom (Streamfunction)')
cccc=sum(Gunsgnnet(:,gicol),2)+tekk+...
-inboxdirr*Amatnet(:,gicol)*bhat(gicol);
ctest=flipud(cumsum(flipud(cccc(1:Nlay-1)),1));
%define "cumulator" matrix
cmat=triu(ones(Nlay,Nlay));
cmat(1:Nlay-1,Nlay)=0; %remove last layer (top to bottom)
Tr{isa}(:,iprop)=cmat*( sum(Gunsgnnet(:,gicol),2)+tekk+...
-inboxdirr*Amatnet(:,gicol)*Binit(gicol) )+Twbc;
Ta{isa}(:,iprop)=cmat*(sum(Gunsgnnet(:,gicol),2)+tekk+...
-inboxdirr*Amatnet(:,gicol)*bhat(gicol) )+Twbc;
%plot(Ta{isa}(:,iprop)-[ctest;cccc(Nlay)]);ppause
if any((Ta{isa}(:,iprop)-[ctest;cccc(Nlay)]-Twbc)>1e-6)
error('(Ta{isa}(:,iprop)-[ctest;cccc(Nlay)]-Twbc)>1e-6')
end
dTa{isa}(:,iprop)=...
full(sqrt(diag(...
cmat*Amatnet(:,gicol)*P(gicol,gicol)*Amatnet(:,gicol)'*cmat')));
end %if ~exist('pltstream')|pltstream==0
end %isa
%VERTICAL TRANSPORTS (Only if original equations!)
if ~exist('p_compw') | p_compw
for ibox=1:Nbox
for ilayint=1:Nlay-2
girow=pm.ibfrow(ibox)-1+[(iprop-1)*Nlay+1:iprop*Nlay];
% iwkzrow=1+girow(ilayint);
iwkzrow=1+ilayint;
iwcol=gifwcol(ibox)+ilayint-1;
Wtrans{ibox}(ilayint,iprop)=full(Amatnet(iwkzrow,iwcol)*bhat(iwcol));
Dwtrans{ibox}(ilayint,iprop)=sqrt(full(...
Amatnet(iwkzrow,iwcol)*P(iwcol,iwcol)*Amatnet(iwkzrow,iwcol)'));
if any(strcmp(fieldnames(boxi.conseq),'Kzstd'))
ikzcol=gifkzcol(ibox)+ilayint-1;
Kztrans{ibox}(ilayint,iprop)=full(Amatnet(iwkzrow,ikzcol)*bhat(ikzcol));
Dkztrans{ibox}(ilayint,iprop)=sqrt(full(...
Amatnet(iwkzrow,ikzcol)*P(ikzcol,ikzcol)*Amatnet(iwkzrow,ikzcol)'));
end
end %ilayint
if iwcol~=gilwcol(ibox) |...
(any(strcmp(fieldnames(boxi.conseq),'Kzstd')) & ikzcol~=gilkzcol(ibox))
error('PROBLEM WITH W or Kz INDICES')
end
if exist('p_savewkztrans') & p_savewkztrans
h_savewkztrans
end
%RESIDUALS FOR CONSERVATION EQUATION
if any(boxi.conseq.freshw)&strcmp(propnm{iprop},'mass')
disp('INCLUDING FRESHWATER IN RESIDUAL COMPUTATION')
fw2subs([1 Nlay])=boxi.conseq.freshw;fw2subs=fw2subs(:);
else
fw2subs=0;
end
if any(strcmp(fieldnames(boxi.conseq),'rhs'))
disp(['Substracting initial rhs to ' boxi.name])
thisrhs=boxi.conseq.rhs(:,iprop);
else
thisrhs=0;
end
Resr{ibox}(:,iprop)=G(girow)+Ekman(girow)-fw2subs+...
Amat(girow,:)*Binit-thisrhs;
Dresr{ibox}(:,iprop)=1./(Rwght(girow)+~Rwght(girow));
Dresr{ibox}(~Rwght(girow))=NaN;
Resa{ibox}(:,iprop)=G(girow)+Ekman(girow)-fw2subs+...
Amat(girow,:)*bhat-thisrhs;
Dres{ibox}(:,iprop)=sqrt(full(diag(Amat(girow,:)*P*Amat(girow,:)')));
if p_res_units & (boxi.harea(1) ~= 0)
switch propnm{iprop}
case {'oxyg','phos','sili','nita'}
disp(['Scaling ' propnm{iprop} ' residual by average layer area'])
scalefac0=1./boxi.harea(1:Nlay-1)';
scalefac0(isinf(scalefac0))=0;
%kmol/s -> mol/m2/yr
scalefac=1e9/1e6*365*24*3600*[scalefac0;1./boxi.harea(1)];
resunit{iprop}='mol yr^{-1}m^{-2}';
if strcmp(propnm{iprop},'oxyg')
resname{iprop}='OUR';
elseif strcmp(propnm{iprop},'heat')
resname{iprop}='Heat flux';
resunit{iprop}='W m^{-2}';
scalefac=1e15*[scalefac0;1./boxi.harea(1)];
else
resname{iprop}=[propnm{iprop} ' utilization'];
end
Resa{ibox}(:,iprop)=-scalefac.*Resa{ibox}(:,iprop);
Resr{ibox}(:,iprop)=-scalefac.*Resr{ibox}(:,iprop);
Dres{ibox}(:,iprop)= scalefac.*Dres{ibox}(:,iprop);
otherwise
resname{iprop}=[propnm{iprop} ' residuals'];
resunit{iprop}=lunits{iprop};
end %switch propnm{iprop}
else
resname{iprop}=[propnm{iprop} ' residuals'];
resunit{iprop}=lunits{iprop};
end %p_res_units
end %ibox
else %if p_compw
Wtrans=NaN;Dwtrans=NaN;
Kztrans=NaN;Dkztrans=NaN;
Resr=NaN;Dresr=NaN;
Resa=NaN;Dres=NaN;
end %if p_compw
end %iprop
%FRESHWATER FLUX
if any(strcmp(fieldnames(boxi.conseq),'freshwstd'))
if ~exist('p_compw') | p_compw
freshwater=boxi.conseq.freshw+bhat(pm.ifw);
dfreshwater=sqrt(P(pm.ifw,pm.ifw));
disp(sprintf('Adjusted Freswater flux (positive=P): %6.2g +/- %6.2g Sv',...
freshwater,dfreshwater))
else
freshwater=NaN;
dfreshwater=NaN;
end
end
%%%%%%%% COMPUTE POSITIVE /NEGATIVE TRANSPORTS FOR EACH SECTION
% CRITERIA FOR WATER MASS TRANSPORT SELECTION
pltwmass=0
if pltwmass
tmin=-5 %POTENTIAL TEMPERATURE
tmax=1.15
smin=0
smax=50
%RESIDUAL WATER MASS
aresb=[];
Resrel=0;
Resb1=0;
gicolb=[];
for isa=1:Nsec
%GET SECTION PAIR DATA
secid=[gsecs.name{isa} gsecs.namesuf{isa}];
datadir=gsecs.datadir{isa};
getpdat
%COMPUTE MATRIX WITH AREAS
p_trig=1; %USE TRIANGLES AT BOTTOM
Dep=sw_dpth(Presctd,Slat(1)); %more accurate with depth
Botdep=sw_dpth(Botp,Slat); %up to (0.6Sv difference)
A0=mk_A0(p_trig,MPres,Dep,Botdep,Slat,Slon,ptemp,Maxd,...
Maxdp,Npair);
%ABSOLUTE VELOCITY IN M/S
gicol=pm.gifcol(isa):(pm.gifcol(isa)+gsecs.npair(isa)-1 ...
-any(strcmp(fieldnames(pm),'giEkcol')));
gicolb=[gicolb,gicol];
bhatsec=bhat(gicol);
Va=grelvel{isa}+ones(MPres(1),1)*bhatsec'/100;
%DENSITY AND POTENTIAL TEMPERATURE
prhoi=sw_dens(psali,ptemp,Presctd);
tpot=sw_ptmp(psali,ptemp,Presctd,0);
%FIND LINEAR COMBINATION COEFF FOR bhat: apos/aneg
apos=zeros(size(bhatsec'));
aneg=zeros(size(bhatsec'));
Trelpos=0;
Trelneg=0;
Tbpos1=0;
ares=zeros(size(bhatsec'));
for id=1:MPres(1)
gipos=find((Va(id,:)>0) & (tpot(id,:)>tmin) & (tpot(id,:)<tmax) &...
(psali(id,:)>smin) & (psali(id,:)<smax));
if ~isempty(gipos)
arhop=A0(id,gipos).*prhoi(id,gipos);
Trelpos=Trelpos+sum(arhop.*grelvel{isa}(id,gipos));
apos(gipos)=apos(gipos)+arhop;
Tbpos1=Tbpos1+arhop*bhatsec(gipos)/100; %FOR TEST
end
gineg=find((Va(id,:)<0) & (tpot(id,:)>tmin) & (tpot(id,:)<tmax) &...
(psali(id,:)>smin) & (psali(id,:)<smax));
if ~isempty(gineg)
arhon=A0(id,gineg).*prhoi(id,gineg);
Trelneg=Trelneg+sum(arhon.*grelvel{isa}(id,gineg));
aneg(gineg)=aneg(gineg)+arhon;
end
%FIND RESIDUAL OF A WATER MASS
gires=find((tpot(id,:)>tmin) & (tpot(id,:)<tmax) &...
(psali(id,:)>smin) & (psali(id,:)<smax));
if ~isempty(gires)
ares_=A0(id,gires).*prhoi(id,gires);
Resrel=Resrel-gsecs.inboxdir(isa)*sum(ares_.*grelvel{isa}(id,gires));
ares(gires)=ares(gires)-gsecs.inboxdir(isa)*ares_;
Resb1=Resb1-gsecs.inboxdir(isa)*ares_*bhatsec(gires)/100;
end
end %id
aresb=[aresb,ares];
%TOTAL TRANPORT:
Tpos{isa}=Trelpos+apos*bhatsec/100;
Dtpos{isa}=sqrt(full(apos*P(gicol,gicol)*apos'))/100;
Tneg{isa}=Trelneg+aneg*bhatsec/100;
Dtneg{isa}=sqrt(full(aneg*P(gicol,gicol)*aneg'))/100;
end %isa
%Check if agrees with total transport
for isa=1:Nsec
disp((Tpos{isa}+Tneg{isa})/1e9+gsecs.EkmanT(isa))
end
%GET TOTAL RESIDUAL WATER MASS
Restot=Resrel+aresb*bhat(gicolb)/100;
Drestot=sqrt(full(aresb*P(gicolb,gicolb)*aresb'))/100;
end %pltwmass=0
%%%%%%%% HERE ADDITIONAL FLUX EQUATION TRANSPORTS
%Calculate Ekman transport based on the flux equations
%Ekman transport recomputation (if used)
if any(strcmp(fieldnames(gsecs),'dEkstd'))
%bxi_showek
elseif p_calculate_Ekman
disp('***************************************************')
disp('Net flux estimation without Ekman nor Rhs:')
disp('Same sign as Ekman flux in principle')
%looks up mass flux equations
for iflx=1:length(pm.eqname)
if findstr(pm.eqname{iflx},'aflxmassnet')
flx=Amat(iflx,:)*bhat+G(iflx);
dflx=sqrt(Amat(iflx,:)*P*Amat(iflx,:)');
disp(sprintf('%s: %4.2g +/- %4.2g', pm.eqname{iflx},flx,dflx))
end
end
end %if p_calculate_Ekman
%%%%%%%% PLOTS
if any(strcmp(fieldnames(boxi),'glevels'))
layids=boxi.glevels;
elseif any(strcmp(fieldnames(boxi),'sigint'))
layids=boxi.sigint;
elseif any(strcmp(fieldnames(boxi),'plevels'))
layids=boxi.plevels;
end
if p_plots
hydrotrans_plot
end
Document Actions
