mkequats.m
by
Webmaster Legos
—
last modified
Jul 05, 2012 12:37 PM
mkequats.m
—
Objective-C source code,
17 kB (17,688 bytes)
File contents
%function mkequats
% KEY:
% USAGE :
%
%
%
%
% DESCRIPTION :
%
%
% INPUT:
%
% OUTPUT:
%
% AUTHOR : A.Ganachaud (ganacho@gulf.mit.edu) , Jul 97
%
% SIDE EFFECTS :
%
% TESTS : see cmp_Amat.m (comparison with previous Amat from Alison)
%
% SEE ALSO :
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALLER:
% CALLEE: mk_set_layprops
Nsec=length(gsecs.name); %number of sections
Nlay=boxi.nlay;%number of layers
Nconseq=length(boxi.conseq.propid); %num. of cons. equations
if ~iscell(propnm) %conversion for compatibility with previous versions
for iprop=1:size(propnm,1)
propnm1{iprop}=propnm(iprop,:);
end
propnm=propnm1;
end
% OPENING A DIARY
drnm=[OPdir boxi.name '_' boxi.modelid '.dry'];
disp(['OPENING DIARY ' drnm])
if exist(drnm)==2
disp('removing old diary')
unix(['\rm ' drnm]);
end
diary(drnm)
clk=fix(clock);
disp([date sprintf(' %i:%i',clk(4:5))])
%Check that there is a top and bottom layer
if any(strcmp(fieldnames(boxi),'glevels'))
if boxi.glevels(1)>22 | boxi.glevels(boxi.nlay)<48
disp(['boxi.glevels must be finished with small and large ' ...
'values to get surface and bottom !'])
ppause
end
elseif any(strcmp(fieldnames(boxi),'sigint'))
if boxi.sigint(1)>22 | boxi.sigint(boxi.nlay)<48
disp(['boxi.sigint must be finished with small and large ' ...
'values to get surface and bottom !'])
ppause
end
elseif ~isstr(boxi.plevels) & (...
boxi.plevels(1)>1 | boxi.plevels(boxi.nlay)<6500)
disp(['boxi.plevels must be finished with small and large ' ...
'values to get surface and bottom !'])
ppause
end
if any(strcmp(fieldnames(gsecs),'dEkstd'))
disp('Ekman transport Adjustment')
p_dEk=1;
else
disp('No Ekman transport Adjustement')
p_dEk=0;
ppause
end
if ~iscell(boxi.conseq.norm)
tmp=boxi.conseq.norm;boxi.conseq=rmfield(boxi.conseq,'norm');
for iprop=1:Nconseq
boxi.conseq.norm{iprop}=tmp;
end
end
if any(strcmp(fieldnames(boxi.conseq),'freshwstd'))
disp('Get Freshwater flux')
p_getFW=1;
else
disp('Does not get Freshwater flux ...')
p_getFW=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PART I:
% FOR EACH SECTION
% COMPUTE INTERFACE POSITIONS
% REFERENCE VELOCITY
% FOR EACH PROPERTY
% GET INTEGRALS IN EACH LAYER (DX*DY*C AND DX*DY*C*VREL)
% GET AVERAGE LAYER PROPERTIES
% GET AVG PROPERTIES FOR THE WHOLE BOX
if ~exist('lprop')
mk_set_layprops
end
% RESULTING INTEGRALS FOR TRANSPORTS AND A MATRIX:
% lprop{isb,iprop}(ipair,ilay) = rho * C * DX * DZ
% lgvprop{isb,iprop}(ipair,ilay) = relvel * rho * C * DX * DZ
% lunits{iprop} = units for lgvprop and lgvprop*b (if b in cm/s)
% lscale(iprop) = value by which was divided lgvprop to get right units
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PART II:
% WRITE CONSERVATION EQUATIONS
% GROUPED BY PROPERTY, NOT BY LAYER AS IN ALISON'S
pm.gifcol=1;
irow=1; %current row
Amat=zeros(Nlay*Nconseq+Nsec*Nconseq,...
sum(gsecs.npair)+... %ref. velocities
any(strcmp(fieldnames(gsecs),'dEkstd'))*length(gsecs.npair)+... %dEkman
(Nlay-2)*(1+any(strcmp(fieldnames(boxi.conseq),'Kzstd')))+...%W/Kzz
p_getFW); %Last is freshwater flux
Gunsgn=zeros(size(Amat));
G=zeros(Nlay*Nconseq+Nsec*Nconseq,1);
Ekman=zeros(Nlay*Nconseq+Nsec*Nconseq,1);
Norm=ones(Nlay*Nconseq+Nsec*Nconseq,1);
Rwght=zeros(Nlay*Nconseq+Nsec*Nconseq,1);
Rhs=zeros(Nlay*Nconseq+Nsec*Nconseq,1);
disp('WRITING CONSERVATION EQUATIONS ...')
for iprop=1:Nconseq
%FILL CONSERVATION EQUATIONS FOR CURRENT PROPERTY
%SECTIONS
for isb=1:Nsec
icol=pm.gifcol(isb); %current column
gicol=icol+gsecs.gip2select{isb}-1; %pair selected
gil=1:Nlay;
girow=irow:irow+Nlay-1;
if any(strcmp(fieldnames(boxi),'gil2mask'))& ...
length(boxi.gil2mask)>=isb&~isempty(boxi.gil2mask{isb})
disp(sprintf('Section %s: masking layers %i to %i in cons. eq.',...
gsecs.name{isb},min(boxi.gil2mask{isb}),max(boxi.gil2mask{isb})))
gil(boxi.gil2mask{isb})=[];
girow=girow(gil);
end
Amat(girow,gicol)=...
-gsecs.inboxdir(isb)*(lprop{isb,iprop}(:,gil)');
Gunsgn(girow,gicol)=...
abs(gsecs.inboxdir(isb))*lgvprop{isb,iprop}(:,gil)';
%abs(gsecs.inboxdir(isb)) USED ONLY TO CANCEL A SECTION IF 0
G(girow)=G(girow)...
-gsecs.inboxdir(isb)*(sum(lgvprop{isb,iprop}(:,gil))');
switch boxi.conseq.norm{iprop}
case 'std' %for use with property anomalies
giboggy=find(imag(boxi.lstd(:,iprop)));
if ~isempty(giboggy)
disp(sprintf('Complex standard deviation found layer %i',giboggy))
disp('Take the one from layer above');ppause
boxi.lstd(giboggy,iprop)=boxi.lstd(giboggy(1)-1,iprop);
end
Norm(irow:irow+Nlay-1)=1e9/1000*boxi.lstd(:,iprop)/lscale(iprop);
%disp(Norm(irow:irow+Nlay-1));ppause
if ~any(strcmp(fieldnames(boxi.conseq),'anom'))...
| length(boxi.conseq.anom)<iprop ...
| boxi.conseq.anom(iprop)~=1
disp(['Warning: using std norm for non-anomaly ' propnm{iprop} ...
' equation']);
end
case 'rms'
Norm(irow:irow+Nlay-1)=1e9/1000*boxi.lrms(:,iprop)/lscale(iprop);
case 'one'
Norm(irow:irow+Nlay-1)=ones(Nlay,1);
otherwise
error(['norm ' boxi.conseq.norm{iprop} ' not programmed'])
end
%EKMAN TRANSPORT IS perEk*TotEk*avglayer
%Remove Zero width layers
glay2sum=find(Lays{isb,1}.lverarea(1:Nlay-1));
glayzero=find(~Lays{isb,1}.lverarea(1:Nlay-1));
ekperc=1e6*gsecs.EkmanT(isb)*[gsecs.perEk(:,isb);...
zeros(Nlay-1-length(gsecs.perEk(:,isb)),1)];
gsecs.Ekt(glay2sum,isb,iprop)=Lays{isb,iprop}.lavgprop(glay2sum)'.*...
ekperc(glay2sum)/lscale(iprop);
gsecs.Ekt(glayzero,isb,iprop)=0;
gsecs.Ekt(Nlay,isb,iprop)=sum(gsecs.Ekt(glay2sum,isb,iprop));
Ekman(irow:irow+Nlay-1)=Ekman(irow:irow+Nlay-1) ...
-gsecs.inboxdir(isb)*gsecs.Ekt(:,isb,iprop);
%WEIGHT
Rwght(irow:irow+Nlay-1)=boxi.conseq.weight(:,iprop);
%FILL PARAMETERS
pm.gilcol(isb)=icol+gsecs.npair(isb)-1;
if p_dEk
pm.giEkcol(isb)=pm.gilcol(isb)+1;
if gsecs.EkmanT(isb)==0
gsecs.EkmanT(isb)=1e-6;
end
Amat(irow:irow+Nlay-1,pm.giEkcol(isb))=...
-gsecs.inboxdir(isb)*gsecs.Ekt(:,isb,iprop)./gsecs.EkmanT(isb);
end
if isb<Nsec
pm.gifcol(isb+1)=pm.gilcol(isb)+1+(p_dEk==1);
end
end %isb
%RHS (OPTIONAL)
if any(strcmp(fieldnames(boxi.conseq),'rhs'))&...
size(boxi.conseq.rhs,2)>=iprop
Rhs([irow:irow+Nlay-1])=boxi.conseq.rhs(:,iprop);
end
%FRESHWATER (OPTIONAL)
if iprop==MASS
if ~isempty(boxi.conseq.freshw)&boxi.conseq.freshw
disp('Adding the a priori freshwater flux to RHS of ')
disp('total and first layer mass equations')
Rhs([irow,irow+Nlay-1])=Rhs([irow,irow+Nlay-1])+boxi.conseq.freshw;
end
end
%VERTICAL ADVECTION
icol=icol+gsecs.npair(isb)+(p_dEk==1);
pm.ifwcol=icol;
for il=2:Nlay-1
%TOP TO BOTTOM HAS NOT W TERMS
%BOTTOM FLUX IN UPPER LAYER
if isinf(boxi.lverarea(il))|boxi.lverarea(il)==0
Amat(irow+il-2,icol)=0;
boxi.conseq.wstd(il-1)=0;
if any(strcmp(fieldnames(boxi.conseq),'Kzstd'))
boxi.conseq.Kzstd(il-1)=0;
end
Rwght(irow+il-1)=0;
if iprop==MASS
disp(sprintf(['NULL BOTTOM LAYER %i: Rwght(%i)=0,wstd(%i)=0 '...
'Kzstd(%i)=0'], il,irow+il-1,il-1))
end
elseif isinf(boxi.lverarea(il-1))|boxi.lverarea(il-1)==0
Amat(irow+il-2,icol)=0;
boxi.conseq.wstd(il-1)=0;
if any(strcmp(fieldnames(boxi.conseq),'Kzstd'))
boxi.conseq.Kzstd(il-1)=0;
end
Rwght(irow+il-2)=0;
if iprop==MASS
disp(sprintf(['NULL SURFACE LAYER %i: Rwght(%i)=0,wstd(%i)=0 '...
'Kzstd(%i)=0'],il,irow+il-2,il-1,il-1))
end
else
Amat(irow+il-2,icol)=-boxi.harea(il)*boxi.liavg(il,iprop)/...
lscale(iprop)/100;
end
Amat(irow+il-1,icol)=-Amat(irow+il-2,icol);
icol=icol+1;
end %il
pm.ilwcol=icol-1;
%FRESHWATER FLUX DETERMINATION
if p_getFW
if iprop==MASS
%Take the Ekman perc outcropping for repartion in the surface layers
percFW=[mean(gsecs.perEk,2);zeros(Nlay-1-size(gsecs.perEk,1),1);1];
Amat(irow:irow+Nlay-1,icol)=-percFW;
pm.ifw=icol;
end
icol=icol+1;
end
%VERTICAL DIFFUSION
if any(strcmp(fieldnames(boxi.conseq),'Kzstd'))&...
~isempty(boxi.conseq.Kzstd) & iprop~=MASS %No mass diffusion
pm.ifKzcol=icol;
for il=2:Nlay-1
%TOP TO BOTTOM HAS NOT W TERMS
%BOTTOM FLUX IN UPPER LAYER
Amat(irow+il-2,icol)=+boxi.harea(il)*boxi.liavgdCrdz(il-1,iprop)/...
lscale(iprop)/1e4; %Kz in cm^2/s
Amat(irow+il-1,icol)=-Amat(irow+il-2,icol);
icol=icol+1;
end %il
pm.ilKzcol=icol-1;
end %if any(strcmp(fieldnames(boxi.conseq),'Kzstd'))
%CREATE ANOMALY EQUATIONS
if iprop~=MASS & any(strcmp(fieldnames(boxi.conseq),'anom'))...
&length(boxi.conseq.anom)>=iprop ...
&boxi.conseq.anom(iprop)==1
if ~strcmp(boxi.conseq.norm{iprop},'std')
disp('warning, boxi.conseq.norm is not std for property anomalies')
switch killblank(boxi.name)
case {'georgia','agulhas','kerguelen','A12I5I9',...
'adelaide','maquarie','sanmartin','P6I5'}
disp(['OK for ' boxi.name ])
otherwise
ppause
end
end
disp(['Anomaly Equation for ' propnm{iprop}])
nAmat=size(Amat,2);
C0=boxi.lavg(:,iprop)*1e9/1000/lscale(iprop);
girowprop=irow:irow+Nlay-1;
girowmass=1:Nlay;
Amat(girowprop,:)=Amat(girowprop,:)-C0*ones(1,nAmat).*Amat(girowmass,:);
Gunsgn(girowprop,:)=Gunsgn(girowprop,:)...
-C0*ones(1,nAmat).*Gunsgn(girowmass,:);
G(girowprop)=G(girowprop)-C0.*G(girowmass);
Ekman(girowprop)=Ekman(girowprop)-C0.*Ekman(girowmass);
Rhs(girowprop)=Rhs(girowprop)-C0.*Rhs(girowmass);
end
%FILL PARAMETERS
for il=1:Nlay-1
pm.eqname{irow+il-1}=sprintf('cons%s-%03i',propnm{iprop},il);
end
if iprop==MASS
pm.eqname{irow+Nlay-1}=sprintf('%s-%s',propnm{iprop},boxi.name);
else
pm.eqname{irow+Nlay-1}=sprintf('cons%s-tot',propnm{iprop});
end
irow=irow+Nlay;
end %on iprop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PART III:
% WRITE FLUX EQUATIONS
% Flux equations are not used as constraints. They contain top-to bottom
% flux for each section and each property. They are not really used
% afterwards and are more here for historical reasons. They also are
% useful to check the Amat coefficients.
disp('WRITTING FLUX EQUATIONS ...')
for isb=1:Nsec
icol=pm.gifcol(isb);
pm.giffrow(isb)=irow;
gicol=icol+gsecs.gip2select{isb}-1; %pair selected
for iprop=1:Nconseq
Amat(irow,gicol)=...
lprop{isb,iprop}(:,Nlay)';
Gunsgn(irow,gicol)=...
lgvprop{isb,iprop}(:,Nlay)';
G(irow)=sum(lgvprop{isb,iprop}(:,Nlay))';
switch boxi.conseq.norm{iprop}
case 'rms'
Norm(irow)=1e9/1000*boxi.lrms(Nlay,iprop)/lscale(iprop);
case 'one'
Norm(irow)=1;
case 'std'
Norm(irow)=1;
if iprop==MASS
disp(['Norm 1 used for flux equation, STD used only in Flux'...
' equations'])
end
end
Ekman(irow)=gsecs.Ekt(Nlay,isb,iprop);
if p_dEk
Amat(irow,pm.giEkcol(isb))=gsecs.Ekt(Nlay,isb,iprop)/gsecs.EkmanT(isb);
end %if p_dEk
pm.eqname{irow}=sprintf('flux%s%s',propnm{iprop},gsecs.name{isb});
irow=irow+1;
end %iprop
pm.gilfrow(isb)=irow-1;
end %isb
% WRITE ADDITIONAL FLUX EQUATIONS
if exist('flxeq')& ~isempty(flxeq)
isbused=[];
for iflx=1:length(flxeq)
if (any(strcmp(fieldnames(flxeq(iflx)),'putEk')) & flxeq(iflx).putEk)|...
(any(strcmp(fieldnames(flxeq(iflx)),'putek')) & flxeq(iflx).putek)|...
(any(strcmp(fieldnames(flxeq(iflx)),'puteq')) & flxeq(iflx).puteq)
disp(sprintf(...
'Include Ekman transport in add. equation %s',flxeq(iflx).name))
putek=1;
else
disp(sprintf(...
'No Ekman transport in add. equation %s',flxeq(iflx).name))
putek=0;
end
if ~isempty(flxeq(iflx).secs)
pm.gifafrow(iflx)=irow;
for iprop_f=1:length(flxeq(iflx).prop)
if iprop==1
if ~isempty(isbused)&any(isb==isbused)
disp(sprintf(...
'WARNING: section %i already used in flux equation %s',...
isb,flxeq(iflx).name))
end
isbused=[isbused,isb];
end
iprop=flxeq(iflx).prop(iprop_f);
G=[G;0];
Ekman=[Ekman;0];
for isf=1:length(flxeq(iflx).secs)
isb=flxeq(iflx).secs(isf);
icol=pm.gifcol(isb);
gicol=flxeq(iflx).fpair(isf):flxeq(iflx).lpair(isf);
if flxeq(iflx).fpair(isf)<min(gsecs.gip2select{isb}) |...
flxeq(iflx).lpair(isf)>max(gsecs.gip2select{isb})
error(['flxeq(iflx).fpair, .lpair must be given in absolute '...
'pair #'])
end
gipairs=find(gsecs.gip2select{isb}==flxeq(iflx).fpair(isf)):...
find(gsecs.gip2select{isb}==flxeq(iflx).lpair(isf));
gilay2sum=flxeq(iflx).toplay:flxeq(iflx).botlay;
Amat(irow,icol+gicol-1)=abs(gsecs.inboxdir(isb))*...
flxeq(iflx).sign(isf)*sum(...
lprop{isb,iprop}(gipairs,gilay2sum)',1);
%MULTIPLY BY abs(gsecs.inboxdir(isb)) TO ALLOW SECTION REMOVAL
Gunsgn(irow,icol+gicol-1)=abs(gsecs.inboxdir(isb))*sum(...
lgvprop{isb,iprop}(gipairs,gilay2sum)',1);
G(irow)=G(irow)+abs(gsecs.inboxdir(isb))*...
flxeq(iflx).sign(isf)*sum(Gunsgn(irow,icol+gicol-1));
if putek
Ekman(irow)=abs(gsecs.inboxdir(isb))*...
Ekman(irow)+flxeq(iflx).sign(isf)*sum(...
gsecs.Ekt(gilay2sum,isb,iprop));
if p_dEk
Amat(irow,pm.giEkcol(isb))=abs(gsecs.inboxdir(isb))*...
flxeq(iflx).sign(isf)*sum(...
gsecs.Ekt(gilay2sum,isb,iprop))/gsecs.EkmanT(isb);
end
else
Ekman(irow)=0;
if p_dEk
Amat(irow,pm.giEkcol(isb))=0;
end
end
end %isf
Rhs(irow)=flxeq(iflx).frhs(iprop_f);
switch boxi.conseq.norm{iprop}
case 'rms'
Norm(irow)=rms(1e9/1000/lscale(iprop)*...
boxi.lrms(flxeq(iflx).toplay:flxeq(iflx).botlay,iprop));
case 'one'
Norm(irow)=1;
case 'std'
Norm(irow)=1;
end
Rwght(irow)=flxeq(iflx).weight(iprop_f)*any(Amat(irow,:));
pm.eqname{irow}=...
sprintf('aflx%s%s',propnm{iprop},flxeq(iflx).name);
irow=irow+1;
end %iprop_f
pm.gilafrow(iflx)=irow-1;
end %~isempty
end %iflx
else
pm.gifafrow=[];
pm.gilafrow=[];
end %exist flxeq
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PART III: FILL COLUMN WEIGHTS
if any(strcmp(fieldnames(gsecs),'colnorm')) %NORMALIZES
if strcmp(gsecs.colnorm,'length') %LENGTH AS RINTOUL'S
girows=find(Rwght);
for icol=1:size(Amat,2)
Cwght(icol)=sqrt(sqrt(1./...
((Amat(girows,icol)./Norm(girows).*Rwght(girows))'*...
(Amat(girows,icol)./Norm(girows).*Rwght(girows)))));
end
else
error(['Norm ' gsecs.colnorm ' not defined'])
end
else %Norm is only the a priori size
Cwght=ones(1,size(Amat,2));
end
%MULTIPLIES THE NORM BY THE A PRIORI SIZE
for isb=1:Nsec
gicol=pm.gifcol(isb):pm.gilcol(isb);
if length(gsecs.bstd{isb})~=gsecs.npair(isb)
error('gsecs.bstd and gsecs.binit must be given for whole section')
end
Cwght(gicol)=Cwght(gicol).*gsecs.bstd{isb}';
if p_dEk
Cwght(pm.giEkcol(isb))=Cwght(pm.giEkcol(isb))*gsecs.dEkstd(isb);
end
end
gicol=pm.ifwcol:pm.ilwcol;
Cwght(gicol)=Cwght(gicol).*boxi.conseq.wstd';
if p_getFW
Cwght(pm.ifw)=Cwght(pm.ifw).*boxi.conseq.freshwstd;
end
if any(strcmp(fieldnames(pm),'ifKzcol'))
gicol=pm.ifKzcol:pm.ilKzcol;
Cwght(gicol)=Cwght(gicol).*boxi.conseq.Kzstd';
end
%FILL IN INITIAL GUESS
for isb=1:Nsec
Binit(pm.gifcol(isb):pm.gilcol(isb))=gsecs.binit{isb};
end
Binit(pm.ifwcol:pm.ilwcol)=zeros(pm.ilwcol-pm.ifwcol+1,1);
if any(strcmp(fieldnames(pm),'ifKzcol'))
Binit(pm.ifKzcol:pm.ilKzcol)=0;
end
if any(strcmp(fieldnames(pm),'ifw'))
Binit(pm.ifw)=0;
end
if p_dEk
disp('change the parameters as if we had an extra pair for each section')
disp('when there is Ekman unknown')
pm.gilcol=pm.gilcol+1;
gsecs.npair=gsecs.npair+1;
end
%CHECK THAT THE DELETED PAIRS WERE DELETED (just an anti boggus)
for isb=1:Nsec
if any(strcmp(fieldnames(gsecs),'pair2mask'))&...
length(gsecs.pair2mask)>=isb&...
~isempty(gsecs.pair2mask{isb})
gcoldeleted=pm.gifcol(isb)-1+gsecs.pair2mask{isb};
if any(Amat(:,gcoldeleted))
error('columns not deleted in Amat: bug in the program')
elseif any(Gunsgn(:,gcoldeleted))
error('columns not deleted in Gunsgn: bug in the program')
elseif any(Cwght(gcoldeleted))|any(Binit(gcoldeleted))
error('columns not deleted in Cwght/Binit: bug in the program')
end
end
end
Amat=sparse(Amat);
Gunsgn=sparse(Gunsgn);
Binit=Binit';
Cwght=Cwght';
disp(['SAVING EQUATIONS IN:'])
opfname=[OPdir boxi.name '_' boxi.modelid '_equn.mat '];
disp(opfname)
eval(['save ' opfname ' gsecs boxi flxeq propnm '...
'Amat Gunsgn G Ekman Norm '...
'Rhs Rwght Cwght Binit pm '...
'lipres rlpres grelvel lunits lscale '...
'Laybs Lays Layprops Laypropbs '])
%spy(Amat,'k');grid on;set(gcf,'papero','land');setlargefig
%spy(Gunsgn,'k');grid on;set(gcf,'papero','land');setlargefig
diary off
disp(['CLOSE DIARY ' drnm ])
Document Actions
