matlab_class/gcmfaces_calc/calc_UV_conv.m

function [fldDIV]=calc_UV_conv(fldU,fldV,varargin);
%object:    compute flow field convergent part (i.e. minus the divergence)
%inputs:    fldU and fldV are transport or velocity fields
%optional:  list_factors is the list of factors that need to 
%               be applied to fldU,fldV. By default it is empty (i.e. {}).
%               The most complete list would be {'dh','dz','hfac'}.
%output:    fldDIV is the convergence (integrated, not averaged, over grid cell area)
%
%notes:     fldU,fldV  that may be
%  either       [A] a 3D vector field
%  or           [B] a 2D vector field
%  
%  in case [A], layer thicknesses = mygrid.DRF; in case [B] layer thickness = 1
%  in any case, the global variable mygrid is supposed to be available

global mygrid;

%initialize output:
n3=max(size(fldU.f1,3),1); n4=max(size(fldV.f1,4),1);

%prepare fldU/fldV:
fldU(isnan(fldU))=0; fldV(isnan(fldV))=0;

%if nargin==3; list_factors=varargin; else; list_factors={'dh','dz','hfac'}; end;
if nargin==3; list_factors=varargin{1}; else; list_factors={}; end;

dxg=mk3D(mygrid.DXG,fldU); dyg=mk3D(mygrid.DYG,fldU); 
if size(fldU.f1,3)==length(mygrid.DRF); drf=mk3D(mygrid.DRF,fldU);
elseif size(fldU.f1,3)==1; drf=fldU; drf(:)=1;
else; error('error in calc_UV_conv: non supported field size\n');
end;
facW=drf; facW(:)=1; facS=facW;
for ii=1:length(list_factors); 
  tmp1=list_factors{ii}; 
  if strcmp(tmp1,'dh'); facW=facW.*dyg; facS=facS.*dxg; 
  elseif strcmp(tmp1,'dz'); facW=facW.*drf; facS=facS.*drf;
  elseif strcmp(tmp1,'hfac'); facW=facW.*mygrid.hFacW; facS=facS.*mygrid.hFacS;
  elseif isempty(tmp1); 1;
  else; fprintf('error in calc_UV_conv: non supported factor\n'); return;
  end;
end; 

for k4=1:n4;
fldU(:,:,:,k4)=fldU(:,:,:,k4).*facW;
fldV(:,:,:,k4)=fldV(:,:,:,k4).*facS;
end;

[FLDU,FLDV]=exch_UV(fldU,fldV);
FLDU(isnan(FLDU))=0; FLDV(isnan(FLDV))=0;

fldDIV=fldU;
for iFace=1:fldDIV.nFaces; 
fldDIV{iFace}=FLDU{iFace}(1:end-1,:,:,:)-FLDU{iFace}(2:end,:,:,:)+...
              FLDV{iFace}(:,1:end-1,:,:)-FLDV{iFace}(:,2:end,:,:);
end;


1	gforget	1.1	function [fldDIV]=calc_UV_conv(fldU,fldV,varargin);
2			%object: compute flow field convergent part (i.e. minus the divergence)
3			%inputs: fldU and fldV are transport or velocity fields
4			%optional: list_factors is the list of factors that need to
5			% be applied to fldU,fldV. By default it is empty (i.e. {}).
6			% The most complete list would be {'dh','dz','hfac'}.
7			%output: fldDIV is the convergence (integrated, not averaged, over grid cell area)
8			%
9			%notes: fldU,fldV that may be
10			% either [A] a 3D vector field
11			% or [B] a 2D vector field
12			%
13			% in case [A], layer thicknesses = mygrid.DRF; in case [B] layer thickness = 1
14			% in any case, the global variable mygrid is supposed to be available
15
16			global mygrid;
17
18			%initialize output:
19			n3=max(size(fldU.f1,3),1); n4=max(size(fldV.f1,4),1);
20
21			%prepare fldU/fldV:
22			fldU(isnan(fldU))=0; fldV(isnan(fldV))=0;
23
24			%if nargin==3; list_factors=varargin; else; list_factors={'dh','dz','hfac'}; end;
25			if nargin==3; list_factors=varargin{1}; else; list_factors={}; end;
26
27			dxg=mk3D(mygrid.DXG,fldU); dyg=mk3D(mygrid.DYG,fldU);
28			if size(fldU.f1,3)==length(mygrid.DRF); drf=mk3D(mygrid.DRF,fldU);
29			elseif size(fldU.f1,3)==1; drf=fldU; drf(:)=1;
30			else; error('error in calc_UV_conv: non supported field size\n');
31			end;
32			facW=drf; facW(:)=1; facS=facW;
33			for ii=1:length(list_factors);
34			tmp1=list_factors{ii};
35			if strcmp(tmp1,'dh'); facW=facW.dyg; facS=facS.dxg;
36			elseif strcmp(tmp1,'dz'); facW=facW.drf; facS=facS.drf;
37			elseif strcmp(tmp1,'hfac'); facW=facW.mygrid.hFacW; facS=facS.mygrid.hFacS;
38			elseif isempty(tmp1); 1;
39			else; fprintf('error in calc_UV_conv: non supported factor\n'); return;
40			end;
41			end;
42
43			for k4=1:n4;
44			fldU(:,:,:,k4)=fldU(:,:,:,k4).*facW;
45			fldV(:,:,:,k4)=fldV(:,:,:,k4).*facS;
46			end;
47
48			[FLDU,FLDV]=exch_UV(fldU,fldV);
49			FLDU(isnan(FLDU))=0; FLDV(isnan(FLDV))=0;
50
51			fldDIV=fldU;
52			for iFace=1:fldDIV.nFaces;
53			fldDIV{iFace}=FLDU{iFace}(1:end-1,:,:,:)-FLDU{iFace}(2:end,:,:,:)+...
54			FLDV{iFace}(:,1:end-1,:,:)-FLDV{iFace}(:,2:end,:,:);
55			end;
56
57