| 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 = ne*ne*6/(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)*nps*4; |
| 77 |
|
| 78 |
% Convert the horrible JPL mdsio layout to six sane faces |
| 79 |
fseek(sfid, offset, 'bof'); |
| 80 |
scube = unmangleJPL1( reshape(fread(sfid, nps, 'real*4'),tx*nt,ty), ... |
| 81 |
ne, tx ); |
| 82 |
|
| 83 |
fseek(tfid, offset, 'bof'); |
| 84 |
tcube = unmangleJPL1( reshape(fread(tfid, nps, 'real*4'),tx*nt,ty), ... |
| 85 |
ne, tx ); |
| 86 |
|
| 87 |
fseek(ufid, offset, 'bof'); |
| 88 |
ucube = unmangleJPL1( reshape(fread(ufid, nps, 'real*4'),tx*nt,ty), ... |
| 89 |
ne, tx ); |
| 90 |
|
| 91 |
fseek(vfid, offset, 'bof'); |
| 92 |
vcube = unmangleJPL1( reshape(fread(vfid, nps, 'real*4'),tx*nt,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) = -rhoConst*rC(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.998*10^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.998*10^2 * 9.81 * R(il-1) * 0.0001; |
| 113 |
bcubem1(:,:,icf) = eh3densjmd95( scube(:,:,icf), ... |
| 114 |
tcube(:,:,icf), p ); |
| 115 |
end |
| 116 |
end |
| 117 |
tmask = double(tcube ~= 0.0); |
| 118 |
bcube = bcube .* tmask; |
| 119 |
bcubem1 = bcubem1 .* tmask; |
| 120 |
|
| 121 |
% Fill the (ne+1)*(ne+1) grid with B,T values |
| 122 |
ni = ne; nip1 = ni + 1; |
| 123 |
nj = ne; njp1 = nj + 1; |
| 124 |
bcubep1 = zeros( [ nip1 njp1 6 ] ); |
| 125 |
for icf = 1:6 |
| 126 |
bcubep1(2:nip1,2:njp1,icf) = bcube(:,:,icf); |
| 127 |
end |
| 128 |
|
| 129 |
% Do the upwind-edge B exchanges |
| 130 |
bcubep1(1,2:njp1,1) = bcube(ni:-1:1,nj,5); % - |
| 131 |
bcubep1(2:nip1,1,1) = bcube(1:ni,nj,6); % - |
| 132 |
bcubep1(1,2:njp1,2) = bcube(ni,1:nj,1); % - |
| 133 |
bcubep1(2:nip1,1,2) = bcube(ni,nj:-1:1,6); % - |
| 134 |
bcubep1(1,2:njp1,3) = bcube(ni:-1:1,nj,1); % - |
| 135 |
bcubep1(2:nip1,1,3) = bcube(1:ni,nj,2); % - |
| 136 |
bcubep1(1,2:njp1,4) = bcube(ni,1:nj,3); % - |
| 137 |
bcubep1(2:nip1,1,4) = bcube(ni,nj:-1:1,2); % - |
| 138 |
bcubep1(1,2:njp1,5) = bcube(ni:-1:1,nj,3); % - |
| 139 |
bcubep1(2:nip1,1,5) = bcube(1:ni,nj,4); % - |
| 140 |
bcubep1(1,2:njp1,6) = bcube(ni,1:nj,5); % - |
| 141 |
bcubep1(2:nip1,1,6) = bcube(ni,nj:-1:1,4); % - |
| 142 |
|
| 143 |
% Get B values on the U grid points |
| 144 |
masku = 1.0 - abs(diff(double(bcubep1(2:nip1,:,:) == 0.0),1,2)); |
| 145 |
diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2); |
| 146 |
bonu = bcube(:,:,:) + masku.*diffu; |
| 147 |
|
| 148 |
% Get B values on the V grid points |
| 149 |
maskv = 1.0 - abs(diff(double(bcubep1(:,2:nip1,:) == 0.0),1,1)); |
| 150 |
diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1); |
| 151 |
bonv = bcube(:,:,:) + maskv.*diffv; |
| 152 |
|
| 153 |
if ((newyear ~= oldyear) && (nm > 0)) || (it == length(ilist)) |
| 154 |
|
| 155 |
% Write if its a new year or the last of all entries |
| 156 |
bout = sprintf('b-ave-%d', oldyear); |
| 157 |
mout = sprintf('bm1-ave-%d', oldyear); |
| 158 |
uout = sprintf('ub-ave-%d', oldyear); |
| 159 |
vout = sprintf('vb-ave-%d', oldyear); |
| 160 |
disp([' ==> Writing files "' bout '" "' mout '" "' uout '" "' vout '"']); |
| 161 |
bido = fopen(bout, 'a', 'ieee-be'); |
| 162 |
mido = fopen(mout, 'a', 'ieee-be'); |
| 163 |
uido = fopen(uout, 'a', 'ieee-be'); |
| 164 |
vido = fopen(vout, 'a', 'ieee-be'); |
| 165 |
mb = mb / nm; |
| 166 |
mbm1 = mbm1 / nm; |
| 167 |
mbu = mbu / nm; |
| 168 |
mbv = mbv / nm; |
| 169 |
fwrite(bido, mb, 'real*4'); |
| 170 |
fwrite(mido, mbm1, 'real*4'); |
| 171 |
fwrite(uido, mbu, 'real*4'); |
| 172 |
fwrite(vido, mbv, 'real*4'); |
| 173 |
fclose(bido); |
| 174 |
fclose(mido); |
| 175 |
fclose(uido); |
| 176 |
fclose(vido); |
| 177 |
|
| 178 |
% Re-zero the accumulators |
| 179 |
mb = zeros(ne,ne,6); |
| 180 |
mbm1 = zeros(ne,ne,6); |
| 181 |
mbu = zeros(ne,ne,6); |
| 182 |
mbv = zeros(ne,ne,6); |
| 183 |
nm = 0; |
| 184 |
|
| 185 |
end |
| 186 |
|
| 187 |
% Accumulate partial sums |
| 188 |
mb = mb + bcube; |
| 189 |
mbm1 = mbm1 + bcubem1; |
| 190 |
mbu = mbu + bonu.*ucube; |
| 191 |
mbv = mbv + bonv.*vcube; |
| 192 |
nm = nm + 1; |
| 193 |
|
| 194 |
oldyear = newyear; |
| 195 |
|
| 196 |
end |
| 197 |
end |
Parent Directory
| 
Revision Log
| 
Revision Graph