You are here: Home / Alexandre Ganachaud / hydrosys / mkequats.m

mkequats.m

by Webmaster Legos last modified Jul 05, 2012 12:37 PM

Objective-C source code icon mkequats.m — Objective-C source code, 17 kB (17688 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

logo cnes logo IRD Logo université de Toulouse Logo université Paul Sabatier Logo CNRS
Logo bibliothèque OBS Logo Observatoire Midi Pyrénées