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 ])