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