MITgcm_contrib/mitgcm_tools/hfac.m

function [hfacC,ddz] = hfac(dz,H,hfacmin,dzmin)
%[hfacc]=hfac(dz,h,hfacmin,dzmin);
%[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
%
% hfacc corresponds to the non-dimensional thickness hFacC in MITgcm
% ddz   is the dimensional thickness, ddz(i,j,k)=hfacc(i,j,k)*dz(k)
%
%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;
%
%Written by adcroft@mit.edu, 2001
%$Header:

[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
  disp(sprintf('Level k=%3i  dz(k)=%8.2f  Full cell',k,dz(k)))
 else
  disp(sprintf('Level k=%3i  dz(k)=%8.2f  Lopping with min[hfacC(:,k)*dz(k)]=%8.2f',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	adcroft	1.1	function [hfacC,ddz] = hfac(dz,H,hfacmin,dzmin)
2			%[hfacc]=hfac(dz,h,hfacmin,dzmin);
3			%[hfacc,ddz]=hfac(dz,h,hfacmin,dzmin);
4			%
5			%Create hfacC(i,j,k) and dz(i,j,k) from dz(k), H(i,j)
6			%using MITgcmUV model parameters hfacmin and dzmin
7			%
8			% hfacc corresponds to the non-dimensional thickness hFacC in MITgcm
9			% ddz is the dimensional thickness, ddz(i,j,k)=hfacc(i,j,k)*dz(k)
10			%
11			%e.g.
12			% [hfacC,ddz] = hfac(dz,H,hfacmin,dzmin)
13			%e.g. To create a "p-mask" from above
14			% pmask=zeros(size(hfacC)); pmask( find(hfacC>0) )=1;
15			%
16			%Written by adcroft@mit.edu, 2001
17			%$Header:
18
19			[nx,ny,nt]=size(H);
20			nz=prod(size(dz));
21			N=[nx ny nt nz];
22
23			zf=[0 -cumsum(dz)];
24
25			hfacC=zeros(N);
26			ddz=zeros(N);
27			for k=1:nz,
28			hFacLim=max([ hfacmin min([1 dzmin/dz(k)]) ]);
29			if hFacLim == 1
30			disp(sprintf('Level k=%3i dz(k)=%8.2f Full cell',k,dz(k)))
31			else
32			disp(sprintf('Level k=%3i dz(k)=%8.2f Lopping with min[hfacC(:,k)dz(k)]=%8.2f',k,dz(k),hFacLimdz(k)))
33			end
34			ddd=(zf(k)-H)/dz(k);
35			ddd(find(ddd > 1)) = 1;
36			%ddd(find(ddd < 0)) = 0;
37			%ddd(find(ddd < hFacLim/2 & ddd ~= 0)) = 0;
38			ddd(find(ddd < hFacLim/2)) = 0; % This should do the job of the above 2 lines
39			ddd(find(ddd >= hFacLim/2 & ddd < hFacLim)) = hFacLim;
40			hfacC(:,:,:,k)=ddd;
41			ddz(:,:,:,k)=ddd*dz(k);
42			end
43
44			hfacC=squeeze(hfacC);
45			ddz=squeeze(ddz);