mlosch/interp_llc/hfac.m

function [hfacC,ddz] = hfac(dz,H,hfacmin,dzmin)
%
% Create hfacC(i,j,k) and dz(i,j,k) from dz(k), H(i,j)
% using MITgcmUV model parameters hfacmin and dzmin
%
% e.g.
%   [hfacC,ddz] = hfac(dz,H,hfacmin,dzmin)
% e.g. To create a "p-mask" from above
%   pmask=zeros(size(hfacC)); pmask( find(hfacC>0) )=1;

[nx,ny,nt]=size(H);
nz=prod(size(dz));
N=[nx ny nt nz];

zf=[0 -cumsum(dz)];

hfacC=zeros(N);
ddz=zeros(N);
for k=1:nz,
 hFacLim=max([ hfacmin min([1 dzmin/dz(k)]) ]);
 if hFacLim == 1
  report('Level k=%3i  dz(k)=%8.2f  Full cell\n',k,dz(k))
 else
  report('Level k=%3i  dz(k)=%8.2f  Lopping with min[hfacC(:,k)*dz(k)]=%8.2f\n',k,dz(k),hFacLim*dz(k))
 end
 ddd=(zf(k)-H)/dz(k);
 ddd(find(ddd > 1)) = 1;
%ddd(find(ddd < 0)) = 0;
%ddd(find(ddd < hFacLim/2 & ddd ~= 0)) = 0;
 ddd(find(ddd < hFacLim/2)) = 0; % This should do the job of the above 2 lines
 ddd(find(ddd >= hFacLim/2 & ddd < hFacLim)) = hFacLim;
 hfacC(:,:,:,k)=ddd;
 ddz(:,:,:,k)=ddd*dz(k);
end

hfacC=squeeze(hfacC);
ddz=squeeze(ddz);
1	mlosch	1.1	function [hfacC,ddz] = hfac(dz,H,hfacmin,dzmin)
2			%
3			% Create hfacC(i,j,k) and dz(i,j,k) from dz(k), H(i,j)
4			% using MITgcmUV model parameters hfacmin and dzmin
5			%
6			% e.g.
7			% [hfacC,ddz] = hfac(dz,H,hfacmin,dzmin)
8			% e.g. To create a "p-mask" from above
9			% pmask=zeros(size(hfacC)); pmask( find(hfacC>0) )=1;
10
11			[nx,ny,nt]=size(H);
12			nz=prod(size(dz));
13			N=[nx ny nt nz];
14
15			zf=[0 -cumsum(dz)];
16
17			hfacC=zeros(N);
18			ddz=zeros(N);
19			for k=1:nz,
20			hFacLim=max([ hfacmin min([1 dzmin/dz(k)]) ]);
21			if hFacLim == 1
22			report('Level k=%3i dz(k)=%8.2f Full cell\n',k,dz(k))
23			else
24			report('Level k=%3i dz(k)=%8.2f Lopping with min[hfacC(:,k)dz(k)]=%8.2f\n',k,dz(k),hFacLimdz(k))
25			end
26			ddd=(zf(k)-H)/dz(k);
27			ddd(find(ddd > 1)) = 1;
28			%ddd(find(ddd < 0)) = 0;
29			%ddd(find(ddd < hFacLim/2 & ddd ~= 0)) = 0;
30			ddd(find(ddd < hFacLim/2)) = 0; % This should do the job of the above 2 lines
31			ddd(find(ddd >= hFacLim/2 & ddd < hFacLim)) = hFacLim;
32			hfacC(:,:,:,k)=ddd;
33			ddz(:,:,:,k)=ddd*dz(k);
34			end
35
36			hfacC=squeeze(hfacC);
37			ddz=squeeze(ddz);