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