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