| 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 |
vnam = [ 'U_'; 'V_'; 'T_'; 'S_'; 'B_'; 'B1'; ... |
| 47 |
'UT'; 'VT'; 'US'; 'VS'; 'UB'; 'VB' ]; |
| 48 |
|
| 49 |
il = 1; |
| 50 |
for il = 1:nlay |
| 51 |
|
| 52 |
% a = sprintf('Slice: %g',il); |
| 53 |
% disp(a); |
| 54 |
|
| 55 |
% Zero the accumulators |
| 56 |
for io = 1:size(vnam,1) |
| 57 |
comm = sprintf('%s_acc = zeros(ne,ne,6);', vnam(io,:)); |
| 58 |
eval(comm) |
| 59 |
end |
| 60 |
|
| 61 |
it = 1; |
| 62 |
for it = 1:length(ilist) |
| 63 |
|
| 64 |
iter = ilist(it); |
| 65 |
str = datestr(startDate + (iter*deltatT-tavefreq/2)/60/60/24); |
| 66 |
dv = datevec(startDate + (iter*deltatT-tavefreq/2)/60/60/24); |
| 67 |
newyear = dv(1); |
| 68 |
|
| 69 |
sfname = sprintf(Spat, ilist(it)); |
| 70 |
tfname = sprintf(Tpat, ilist(it)); |
| 71 |
ufname = sprintf(Upat, ilist(it)); |
| 72 |
vfname = sprintf(Vpat, ilist(it)); |
| 73 |
disp([ a ' ' str ' ' tfname ]); |
| 74 |
sfid = fopen(sfname, 'r', 'ieee-be'); |
| 75 |
tfid = fopen(tfname, 'r', 'ieee-be'); |
| 76 |
ufid = fopen(ufname, 'r', 'ieee-be'); |
| 77 |
vfid = fopen(vfname, 'r', 'ieee-be'); |
| 78 |
offset = (il - 1)*nps*4; |
| 79 |
|
| 80 |
% Convert the horrible JPL mdsio layout to six sane faces |
| 81 |
fseek(sfid, offset, 'bof'); |
| 82 |
scube = unmangleJPL1( reshape(fread(sfid, nps, 'real*4'),tx*nt,ty), ... |
| 83 |
ne, tx ); |
| 84 |
|
| 85 |
fseek(tfid, offset, 'bof'); |
| 86 |
tcube = unmangleJPL1( reshape(fread(tfid, nps, 'real*4'),tx*nt,ty), ... |
| 87 |
ne, tx ); |
| 88 |
|
| 89 |
fseek(ufid, offset, 'bof'); |
| 90 |
ucube = unmangleJPL1( reshape(fread(ufid, nps, 'real*4'),tx*nt,ty), ... |
| 91 |
ne, tx ); |
| 92 |
|
| 93 |
fseek(vfid, offset, 'bof'); |
| 94 |
vcube = unmangleJPL1( reshape(fread(vfid, nps, 'real*4'),tx*nt,ty), ... |
| 95 |
ne, tx ); |
| 96 |
|
| 97 |
fclose(sfid); |
| 98 |
fclose(tfid); |
| 99 |
fclose(ufid); |
| 100 |
fclose(vfid); |
| 101 |
|
| 102 |
% Compute the density |
| 103 |
bcube = zeros(size(scube)); |
| 104 |
bcubem1 = zeros(size(scube)); |
| 105 |
for icf = 1:6 |
| 106 |
% pressure_for_eos.F : |
| 107 |
% locPres(i,j) = -rhoConst*rC(k)*gravity*maskC(i,j,k,bi,bj) |
| 108 |
% rhoConst = 9.997999999999999E+02 /* Ref density ( kg/m^3 ) */ |
| 109 |
% gravity = 9.810000000000000E+00 /* Grav accel ( m/s^2 ) */ |
| 110 |
p = 0.09998 * 9.81 * R(il); |
| 111 |
bcube(:,:,icf) = ... |
| 112 |
eh3densjmd95( squeeze(scube(:,:,icf)), ... |
| 113 |
squeeze(tcube(:,:,icf)), p ) * -9.81 / 1000; |
| 114 |
if il > 1 |
| 115 |
p = 0.09998 * 9.81 * R(il-1); |
| 116 |
bcubem1(:,:,icf) = ... |
| 117 |
eh3densjmd95( scube(:,:,icf), ... |
| 118 |
tcube(:,:,icf), p ) * -9.81 / 1000; |
| 119 |
end |
| 120 |
end |
| 121 |
tmask = double(tcube ~= 0.0); |
| 122 |
scube = scube .* tmask; |
| 123 |
bcube = bcube .* tmask; |
| 124 |
bcubem1 = bcubem1 .* tmask; |
| 125 |
|
| 126 |
% Fill the (ne+1)*(ne+1) grid with T,S,B values |
| 127 |
ni = ne; nip1 = ni + 1; |
| 128 |
nj = ne; njp1 = nj + 1; |
| 129 |
tcubep1 = zeros( [ nip1 njp1 6 ] ); |
| 130 |
scubep1 = zeros( [ nip1 njp1 6 ] ); |
| 131 |
bcubep1 = zeros( [ nip1 njp1 6 ] ); |
| 132 |
for icf = 1:6 |
| 133 |
tcubep1(2:nip1,2:njp1,icf) = tcube(:,:,icf); |
| 134 |
scubep1(2:nip1,2:njp1,icf) = scube(:,:,icf); |
| 135 |
bcubep1(2:nip1,2:njp1,icf) = bcube(:,:,icf); |
| 136 |
end |
| 137 |
|
| 138 |
% Do the upwind-edge T exchanges |
| 139 |
tcubep1(1,2:njp1,1) = tcube(ni:-1:1,nj,5); % - |
| 140 |
tcubep1(2:nip1,1,1) = tcube(1:ni,nj,6); % - |
| 141 |
tcubep1(1,2:njp1,2) = tcube(ni,1:nj,1); % - |
| 142 |
tcubep1(2:nip1,1,2) = tcube(ni,nj:-1:1,6); % - |
| 143 |
tcubep1(1,2:njp1,3) = tcube(ni:-1:1,nj,1); % - |
| 144 |
tcubep1(2:nip1,1,3) = tcube(1:ni,nj,2); % - |
| 145 |
tcubep1(1,2:njp1,4) = tcube(ni,1:nj,3); % - |
| 146 |
tcubep1(2:nip1,1,4) = tcube(ni,nj:-1:1,2); % - |
| 147 |
tcubep1(1,2:njp1,5) = tcube(ni:-1:1,nj,3); % - |
| 148 |
tcubep1(2:nip1,1,5) = tcube(1:ni,nj,4); % - |
| 149 |
tcubep1(1,2:njp1,6) = tcube(ni,1:nj,5); % - |
| 150 |
tcubep1(2:nip1,1,6) = tcube(ni,nj:-1:1,4); % - |
| 151 |
|
| 152 |
% Do the upwind-edge S exchanges |
| 153 |
scubep1(1,2:njp1,1) = scube(ni:-1:1,nj,5); % - |
| 154 |
scubep1(2:nip1,1,1) = scube(1:ni,nj,6); % - |
| 155 |
scubep1(1,2:njp1,2) = scube(ni,1:nj,1); % - |
| 156 |
scubep1(2:nip1,1,2) = scube(ni,nj:-1:1,6); % - |
| 157 |
scubep1(1,2:njp1,3) = scube(ni:-1:1,nj,1); % - |
| 158 |
scubep1(2:nip1,1,3) = scube(1:ni,nj,2); % - |
| 159 |
scubep1(1,2:njp1,4) = scube(ni,1:nj,3); % - |
| 160 |
scubep1(2:nip1,1,4) = scube(ni,nj:-1:1,2); % - |
| 161 |
scubep1(1,2:njp1,5) = scube(ni:-1:1,nj,3); % - |
| 162 |
scubep1(2:nip1,1,5) = scube(1:ni,nj,4); % - |
| 163 |
scubep1(1,2:njp1,6) = scube(ni,1:nj,5); % - |
| 164 |
scubep1(2:nip1,1,6) = scube(ni,nj:-1:1,4); % - |
| 165 |
|
| 166 |
% Do the upwind-edge B exchanges |
| 167 |
bcubep1(1,2:njp1,1) = bcube(ni:-1:1,nj,5); % - |
| 168 |
bcubep1(2:nip1,1,1) = bcube(1:ni,nj,6); % - |
| 169 |
bcubep1(1,2:njp1,2) = bcube(ni,1:nj,1); % - |
| 170 |
bcubep1(2:nip1,1,2) = bcube(ni,nj:-1:1,6); % - |
| 171 |
bcubep1(1,2:njp1,3) = bcube(ni:-1:1,nj,1); % - |
| 172 |
bcubep1(2:nip1,1,3) = bcube(1:ni,nj,2); % - |
| 173 |
bcubep1(1,2:njp1,4) = bcube(ni,1:nj,3); % - |
| 174 |
bcubep1(2:nip1,1,4) = bcube(ni,nj:-1:1,2); % - |
| 175 |
bcubep1(1,2:njp1,5) = bcube(ni:-1:1,nj,3); % - |
| 176 |
bcubep1(2:nip1,1,5) = bcube(1:ni,nj,4); % - |
| 177 |
bcubep1(1,2:njp1,6) = bcube(ni,1:nj,5); % - |
| 178 |
bcubep1(2:nip1,1,6) = bcube(ni,nj:-1:1,4); % - |
| 179 |
|
| 180 |
% Get T values on the U,V grid points |
| 181 |
masku = 1.0 - abs(diff(double(tcubep1(2:nip1,:,:) == 0.0),1,2)); |
| 182 |
maskv = 1.0 - abs(diff(double(tcubep1(:,2:nip1,:) == 0.0),1,1)); |
| 183 |
diffu = 0.5*diff(tcubep1(2:nip1,:,:),1,2); |
| 184 |
diffv = 0.5*diff(tcubep1(:,2:nip1,:),1,1); |
| 185 |
tonu = tcube(:,:,:) + masku.*diffu; |
| 186 |
tonv = tcube(:,:,:) + maskv.*diffv; |
| 187 |
|
| 188 |
% Get S values on the U,V grid points |
| 189 |
diffu = 0.5*diff(scubep1(2:nip1,:,:),1,2); |
| 190 |
diffv = 0.5*diff(scubep1(:,2:nip1,:),1,1); |
| 191 |
sonu = scube(:,:,:) + masku.*diffu; |
| 192 |
sonv = scube(:,:,:) + maskv.*diffv; |
| 193 |
|
| 194 |
% Get B values on the U grid points |
| 195 |
diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2); |
| 196 |
diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1); |
| 197 |
bonu = bcube(:,:,:) + masku.*diffu; |
| 198 |
bonv = bcube(:,:,:) + maskv.*diffv; |
| 199 |
|
| 200 |
if ((newyear ~= oldyear) && (nm > 0)) || (it == length(ilist)) |
| 201 |
|
| 202 |
% Write if its a new year or the last of all entries |
| 203 |
disp([' ==> Writing files for ==> ' sprintf('%d',oldyear)]); |
| 204 |
for io = 1:size(vnam,1) |
| 205 |
comm = sprintf( ... |
| 206 |
'%s_ido = fopen(''%s_tave_%d'', ''a'', ''ieee-be'');', ... |
| 207 |
vnam(io,:), vnam(io,:), oldyear ); |
| 208 |
eval(comm) |
| 209 |
comm = sprintf( '%s_acc = %s_acc ./ nm;',vnam(io,:), ... |
| 210 |
vnam(io,:)); |
| 211 |
eval(comm) |
| 212 |
comm = sprintf('fwrite(%s_ido, %s_acc, ''real*4'');', ... |
| 213 |
vnam(io,:), vnam(io,:) ); |
| 214 |
eval(comm) |
| 215 |
comm = sprintf( 'fclose(%s_ido);', vnam(io,:) ); |
| 216 |
eval(comm) |
| 217 |
|
| 218 |
% Re-zero the accumulators |
| 219 |
comm = sprintf('%s_acc = zeros(ne,ne,6);', vnam(io,:)); |
| 220 |
eval(comm) |
| 221 |
end |
| 222 |
|
| 223 |
nm = 0; |
| 224 |
|
| 225 |
end |
| 226 |
|
| 227 |
% Accumulate partial sums |
| 228 |
U__acc = U__acc + ucube; |
| 229 |
V__acc = V__acc + vcube; |
| 230 |
T__acc = T__acc + tcube; |
| 231 |
S__acc = S__acc + scube; |
| 232 |
B__acc = B__acc + bcube; |
| 233 |
B1_acc = B1_acc + bcubem1; |
| 234 |
|
| 235 |
UT_acc = UT_acc + tonu.*ucube; |
| 236 |
VT_acc = VT_acc + tonv.*vcube; |
| 237 |
US_acc = US_acc + sonu.*ucube; |
| 238 |
VS_acc = VS_acc + sonv.*vcube; |
| 239 |
UB_acc = UB_acc + bonu.*ucube; |
| 240 |
VB_acc = VB_acc + bonv.*vcube; |
| 241 |
|
| 242 |
nm = nm + 1; |
| 243 |
|
| 244 |
oldyear = newyear; |
| 245 |
|
| 246 |
end |
| 247 |
end |
| 248 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph