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