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),hFacLim*dz(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); |