high_res_cube/eddy_flux/calcUVBave.m

function m = calcUVBave(Spat, Tpat, Upat, Vpat, ilist, nlay, ne)

% Function m = calcTBave(Spat, Tpat, Upat, Vpat, ilist, nlay, ne)
% 
%
% INPUTS
% Tpat     string containing T-data file pattern
% Upat     string containing U-data file pattern
% Vpat     string containing V-data file pattern
% ilist    list of iteration values
% nlay     number of z layers
% ne       number of values alone each edge of each face
%
% OUTPUTS
% m        output array of dimension ns*nps

if nargin ~= 7
  disp('Error: wrong number of arguments')
  exit
end

deltatT = 1200;                % model time step size (s)
tavefreq = 259200;             % averaging period (s)
startDate = datenum(1992,1,1); % model integration starting date

oldyear = -1;
newyear = oldyear;

nps = 6 * ne*ne;
tx = 85;
ty = 85;
nt = ne*ne*6/(tx*ty);
cx = ne;
cy = ne;
nz = 1;

delR   = [ 
    10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
    10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
    31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
    93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
    139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
    341.50,364.50,387.50,410.50,433.50,456.50 ];
R = cumsum(delR) - 0.5*delR;

il = 1;
for il = 1:nlay
  
  a = sprintf('Slice: %g',il);
  %  disp(a);
  
  %  Zero the accumulators
  nm = 0;
  mb = zeros(ne,ne,6);
  mbm1 = zeros(ne,ne,6);
  mbu = zeros(ne,ne,6);
  mbv = zeros(ne,ne,6);
  
  it = 1;
  for it = 1:length(ilist)
    
    iter = ilist(it);
    str = datestr(startDate + (iter*deltatT-tavefreq/2)/60/60/24);
    dv = datevec(startDate + (iter*deltatT-tavefreq/2)/60/60/24);
    newyear = dv(1);
    
    sfname = sprintf(Spat, ilist(it));
    tfname = sprintf(Tpat, ilist(it));
    ufname = sprintf(Upat, ilist(it));
    vfname = sprintf(Vpat, ilist(it));
    disp([ a '  ' str '  ' tfname ]);
    sfid = fopen(sfname, 'r', 'ieee-be');
    tfid = fopen(tfname, 'r', 'ieee-be');
    ufid = fopen(ufname, 'r', 'ieee-be');
    vfid = fopen(vfname, 'r', 'ieee-be');
    offset = (il - 1)*nps*4;
    
    %  Convert the horrible JPL mdsio layout to six sane faces
    fseek(sfid, offset, 'bof');
    scube = unmangleJPL1( reshape(fread(sfid, nps, 'real*4'),tx*nt,ty), ...
                          ne, tx );
    
    fseek(tfid, offset, 'bof');
    tcube = unmangleJPL1( reshape(fread(tfid, nps, 'real*4'),tx*nt,ty), ...
                          ne, tx );

    fseek(ufid, offset, 'bof');
    ucube = unmangleJPL1( reshape(fread(ufid, nps, 'real*4'),tx*nt,ty), ...
                          ne, tx );
    
    fseek(vfid, offset, 'bof');
    vcube = unmangleJPL1( reshape(fread(vfid, nps, 'real*4'),tx*nt,ty), ...
                          ne, tx );

    fclose(sfid);
    fclose(tfid);
    fclose(ufid);
    fclose(vfid);
    
    %  Compute the density
    bcube   = zeros(size(scube));
    bcubem1 = zeros(size(scube));
    for icf = 1:6
      % pressure_for_eos.F :
      %   locPres(i,j) = -rhoConst*rC(k)*gravity*maskC(i,j,k,bi,bj)
      %   rhoConst = 9.997999999999999E+02  /* Ref density ( kg/m^3 ) */
      %   gravity  = 9.810000000000000E+00  /* Grav accel ( m/s^2 ) */
      p = 9.998*10^2 * 9.81 * R(il) * 0.0001;
      bcube(:,:,icf) = eh3densjmd95( squeeze(scube(:,:,icf)), ...
                                     squeeze(tcube(:,:,icf)), p );
      if il > 1
        p = -9.998*10^2 * 9.81 * R(il-1) * 0.0001;
        bcubem1(:,:,icf) = eh3densjmd95( scube(:,:,icf), ...
                                         tcube(:,:,icf), p );
      end
    end
    
    %  Fill the (ne+1)*(ne+1) grid with T values
    ni = ne;    nip1 = ni + 1;
    nj = ne;    njp1 = nj + 1;
    bcubep1 = zeros( [ nip1 njp1 6 ] );
    for icf = 1:6
      bcubep1(2:nip1,2:njp1,icf) = bcube(:,:,icf);
    end

    %  Do the upwind-edge B exchanges
    bcubep1(1,2:njp1,1) = bcube(ni:-1:1,nj,5);   % - 
    bcubep1(2:nip1,1,1) = bcube(1:ni,nj,6);      % -
    bcubep1(1,2:njp1,2) = bcube(ni,1:nj,1);      % -
    bcubep1(2:nip1,1,2) = bcube(ni,nj:-1:1,6);   % -
    bcubep1(1,2:njp1,3) = bcube(ni:-1:1,nj,1);   % -
    bcubep1(2:nip1,1,3) = bcube(1:ni,nj,2);      % -
    bcubep1(1,2:njp1,4) = bcube(ni,1:nj,3);      % -
    bcubep1(2:nip1,1,4) = bcube(ni,nj:-1:1,2);   % -
    bcubep1(1,2:njp1,5) = bcube(ni:-1:1,nj,3);   % -
    bcubep1(2:nip1,1,5) = bcube(1:ni,nj,4);      % -
    bcubep1(1,2:njp1,6) = bcube(ni,1:nj,5);      % -
    bcubep1(2:nip1,1,6) = bcube(ni,nj:-1:1,4);   % -

    %  Get B values on the U grid points
    masku = 1.0 - abs(diff(double(bcubep1(2:nip1,:,:) == 0.0),1,2));
    diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2);
    bonu = bcube(:,:,:) + masku.*diffu;
    
    %  Get B values on the V grid points
    maskv = 1.0 - abs(diff(double(bcubep1(:,2:nip1,:) == 0.0),1,1));
    diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1);
    bonv = bcube(:,:,:) + maskv.*diffv;
    
    if ((newyear ~= oldyear) && (nm > 0)) || (it == length(ilist))

      %  Write if its a new year or the last of all entries
      bout = sprintf('b-ave-%d', oldyear);
      mout = sprintf('bm1-ave-%d', oldyear);
      uout = sprintf('ub-ave-%d', oldyear);
      vout = sprintf('vb-ave-%d', oldyear);
      disp(['  ==>  Writing files "' bout '" "' mout '" "' uout '" "' vout '"']);
      bido = fopen(bout, 'a', 'ieee-be');
      mido = fopen(mout, 'a', 'ieee-be');
      uido = fopen(uout, 'a', 'ieee-be');
      vido = fopen(vout, 'a', 'ieee-be');
      mb = mb / nm;
      mbm1 = mbm1 / nm;
      mbu = mbu / nm;
      mbv = mbv / nm;
      fwrite(bido, mb,   'real*4');
      fwrite(mido, mbm1, 'real*4');
      fwrite(uido, mbu,  'real*4');
      fwrite(vido, mbv,  'real*4');
      fclose(bido);
      fclose(mido);
      fclose(uido);
      fclose(vido);
      
      %  Re-zero the accumulators
      mb = zeros(ne,ne,6);
      mbm1 = zeros(ne,ne,6);
      mbu = zeros(ne,ne,6);
      mbv = zeros(ne,ne,6);
      nm = 0;
      
    end
    
    %  Accumulate partial sums
    mb = mb + bcube;
    mbm1 = mbm1 + bcubem1;
    mbu = mbu + bonu.*ucube;
    mbv = mbv + bonv.*vcube;
    nm = nm + 1;
    
    oldyear = newyear;
    
  end
end
1	function m = calcUVBave(Spat, Tpat, Upat, Vpat, ilist, nlay, ne)
2
3	% Function m = calcTBave(Spat, Tpat, Upat, Vpat, ilist, nlay, ne)
4	%
5	%
6	% INPUTS
7	% Tpat string containing T-data file pattern
8	% Upat string containing U-data file pattern
9	% Vpat string containing V-data file pattern
10	% ilist list of iteration values
11	% nlay number of z layers
12	% ne number of values alone each edge of each face
13	%
14	% OUTPUTS
15	% m output array of dimension ns*nps
16
17	if nargin ~= 7
18	disp('Error: wrong number of arguments')
19	exit
20	end
21
22	deltatT = 1200; % model time step size (s)
23	tavefreq = 259200; % averaging period (s)
24	startDate = datenum(1992,1,1); % model integration starting date
25
26	oldyear = -1;
27	newyear = oldyear;
28
29	nps = 6 * ne*ne;
30	tx = 85;
31	ty = 85;
32	nt = nene6/(tx*ty);
33	cx = ne;
34	cy = ne;
35	nz = 1;
36
37	delR = [
38	10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
39	10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
40	31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
41	93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
42	139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
43	341.50,364.50,387.50,410.50,433.50,456.50 ];
44	R = cumsum(delR) - 0.5*delR;
45
46	il = 1;
47	for il = 1:nlay
48
49	a = sprintf('Slice: %g',il);
50	% disp(a);
51
52	% Zero the accumulators
53	nm = 0;
54	mb = zeros(ne,ne,6);
55	mbm1 = zeros(ne,ne,6);
56	mbu = zeros(ne,ne,6);
57	mbv = zeros(ne,ne,6);
58
59	it = 1;
60	for it = 1:length(ilist)
61
62	iter = ilist(it);
63	str = datestr(startDate + (iter*deltatT-tavefreq/2)/60/60/24);
64	dv = datevec(startDate + (iter*deltatT-tavefreq/2)/60/60/24);
65	newyear = dv(1);
66
67	sfname = sprintf(Spat, ilist(it));
68	tfname = sprintf(Tpat, ilist(it));
69	ufname = sprintf(Upat, ilist(it));
70	vfname = sprintf(Vpat, ilist(it));
71	disp([ a ' ' str ' ' tfname ]);
72	sfid = fopen(sfname, 'r', 'ieee-be');
73	tfid = fopen(tfname, 'r', 'ieee-be');
74	ufid = fopen(ufname, 'r', 'ieee-be');
75	vfid = fopen(vfname, 'r', 'ieee-be');
76	offset = (il - 1)nps4;
77
78	% Convert the horrible JPL mdsio layout to six sane faces
79	fseek(sfid, offset, 'bof');
80	scube = unmangleJPL1( reshape(fread(sfid, nps, 'real4'),txnt,ty), ...
81	ne, tx );
82
83	fseek(tfid, offset, 'bof');
84	tcube = unmangleJPL1( reshape(fread(tfid, nps, 'real4'),txnt,ty), ...
85	ne, tx );
86
87	fseek(ufid, offset, 'bof');
88	ucube = unmangleJPL1( reshape(fread(ufid, nps, 'real4'),txnt,ty), ...
89	ne, tx );
90
91	fseek(vfid, offset, 'bof');
92	vcube = unmangleJPL1( reshape(fread(vfid, nps, 'real4'),txnt,ty), ...
93	ne, tx );
94
95	fclose(sfid);
96	fclose(tfid);
97	fclose(ufid);
98	fclose(vfid);
99
100	% Compute the density
101	bcube = zeros(size(scube));
102	bcubem1 = zeros(size(scube));
103	for icf = 1:6
104	% pressure_for_eos.F :
105	% locPres(i,j) = -rhoConstrC(k)gravity*maskC(i,j,k,bi,bj)
106	% rhoConst = 9.997999999999999E+02 /* Ref density ( kg/m^3 ) */
107	% gravity = 9.810000000000000E+00 /* Grav accel ( m/s^2 ) */
108	p = 9.99810^2 9.81 * R(il) * 0.0001;
109	bcube(:,:,icf) = eh3densjmd95( squeeze(scube(:,:,icf)), ...
110	squeeze(tcube(:,:,icf)), p );
111	if il > 1
112	p = -9.99810^2 9.81 * R(il-1) * 0.0001;
113	bcubem1(:,:,icf) = eh3densjmd95( scube(:,:,icf), ...
114	tcube(:,:,icf), p );
115	end
116	end
117
118	% Fill the (ne+1)*(ne+1) grid with T values
119	ni = ne; nip1 = ni + 1;
120	nj = ne; njp1 = nj + 1;
121	bcubep1 = zeros( [ nip1 njp1 6 ] );
122	for icf = 1:6
123	bcubep1(2:nip1,2:njp1,icf) = bcube(:,:,icf);
124	end
125
126	% Do the upwind-edge B exchanges
127	bcubep1(1,2:njp1,1) = bcube(ni:-1:1,nj,5); % -
128	bcubep1(2:nip1,1,1) = bcube(1:ni,nj,6); % -
129	bcubep1(1,2:njp1,2) = bcube(ni,1:nj,1); % -
130	bcubep1(2:nip1,1,2) = bcube(ni,nj:-1:1,6); % -
131	bcubep1(1,2:njp1,3) = bcube(ni:-1:1,nj,1); % -
132	bcubep1(2:nip1,1,3) = bcube(1:ni,nj,2); % -
133	bcubep1(1,2:njp1,4) = bcube(ni,1:nj,3); % -
134	bcubep1(2:nip1,1,4) = bcube(ni,nj:-1:1,2); % -
135	bcubep1(1,2:njp1,5) = bcube(ni:-1:1,nj,3); % -
136	bcubep1(2:nip1,1,5) = bcube(1:ni,nj,4); % -
137	bcubep1(1,2:njp1,6) = bcube(ni,1:nj,5); % -
138	bcubep1(2:nip1,1,6) = bcube(ni,nj:-1:1,4); % -
139
140	% Get B values on the U grid points
141	masku = 1.0 - abs(diff(double(bcubep1(2:nip1,:,:) == 0.0),1,2));
142	diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2);
143	bonu = bcube(:,:,:) + masku.*diffu;
144
145	% Get B values on the V grid points
146	maskv = 1.0 - abs(diff(double(bcubep1(:,2:nip1,:) == 0.0),1,1));
147	diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1);
148	bonv = bcube(:,:,:) + maskv.*diffv;
149
150	if ((newyear ~= oldyear) && (nm > 0)) \|\| (it == length(ilist))
151
152	% Write if its a new year or the last of all entries
153	bout = sprintf('b-ave-%d', oldyear);
154	mout = sprintf('bm1-ave-%d', oldyear);
155	uout = sprintf('ub-ave-%d', oldyear);
156	vout = sprintf('vb-ave-%d', oldyear);
157	disp([' ==> Writing files "' bout '" "' mout '" "' uout '" "' vout '"']);
158	bido = fopen(bout, 'a', 'ieee-be');
159	mido = fopen(mout, 'a', 'ieee-be');
160	uido = fopen(uout, 'a', 'ieee-be');
161	vido = fopen(vout, 'a', 'ieee-be');
162	mb = mb / nm;
163	mbm1 = mbm1 / nm;
164	mbu = mbu / nm;
165	mbv = mbv / nm;
166	fwrite(bido, mb, 'real*4');
167	fwrite(mido, mbm1, 'real*4');
168	fwrite(uido, mbu, 'real*4');
169	fwrite(vido, mbv, 'real*4');
170	fclose(bido);
171	fclose(mido);
172	fclose(uido);
173	fclose(vido);
174
175	% Re-zero the accumulators
176	mb = zeros(ne,ne,6);
177	mbm1 = zeros(ne,ne,6);
178	mbu = zeros(ne,ne,6);
179	mbv = zeros(ne,ne,6);
180	nm = 0;
181
182	end
183
184	% Accumulate partial sums
185	mb = mb + bcube;
186	mbm1 = mbm1 + bcubem1;
187	mbu = mbu + bonu.*ucube;
188	mbv = mbv + bonv.*vcube;
189	nm = nm + 1;
190
191	oldyear = newyear;
192
193	end
194	end