| 1 |
%======================================================= |
| 2 |
% |
| 3 |
% $Header: $ |
| 4 |
% |
| 5 |
% Ed Hill |
| 6 |
% |
| 7 |
|
| 8 |
% Check the (u'b'),(v'b') values using the linear EoS: |
| 9 |
% |
| 10 |
% v'b' =approx= g * Alpha * v't' - g * Beta * v's' |
| 11 |
% |
| 12 |
% where: |
| 13 |
% g = 9.81 |
| 14 |
% Alpha = 2.10^{-4} |
| 15 |
% Beta = 7.10^{-4} |
| 16 |
|
| 17 |
clear all ; close all |
| 18 |
|
| 19 |
fnam = [ 'U_'; 'V_'; 'S_'; 'T_'; 'B_'; 'US'; 'UT'; 'UB'; 'VS'; 'VT'; 'VB' ]; |
| 20 |
vnam = lower(fnam); |
| 21 |
|
| 22 |
for i = 1:size(fnam,1) |
| 23 |
c = sprintf(['%sfid = ' ... |
| 24 |
'fopen(''%s_tave_1996'',''r'',''ieee-be'');'], ... |
| 25 |
vnam(i,:),fnam(i,:) ); |
| 26 |
disp(c); |
| 27 |
eval(c); |
| 28 |
end |
| 29 |
|
| 30 |
delR = [ |
| 31 |
10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ... |
| 32 |
10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ... |
| 33 |
31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ... |
| 34 |
93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ... |
| 35 |
139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ... |
| 36 |
341.50,364.50,387.50,410.50,433.50,456.50 ]; |
| 37 |
R = cumsum(delR) - 0.5*delR; |
| 38 |
|
| 39 |
nztot = 50; |
| 40 |
ne = 510; |
| 41 |
nps = 6 * ne * ne; |
| 42 |
iz = 1; |
| 43 |
for iz = 1:nztot |
| 44 |
|
| 45 |
offset = (iz - 1)*nps*4; |
| 46 |
|
| 47 |
for i = 1:size(fnam,1) |
| 48 |
c = sprintf('fseek(%sfid, offset, ''bof'');',vnam(i,:)); |
| 49 |
disp(c); |
| 50 |
eval(c); |
| 51 |
c = sprintf(['%s = reshape(fread(%sfid, nps,' ... |
| 52 |
'''real*4''),ne,ne,6);'],vnam(i,:),vnam(i,:)); |
| 53 |
disp(c); |
| 54 |
eval(c); |
| 55 |
end |
| 56 |
|
| 57 |
% Fill the (ne+1)*(ne+1) grid with T,S,B values |
| 58 |
ni = ne; nip1 = ni + 1; |
| 59 |
nj = ne; njp1 = nj + 1; |
| 60 |
tcubep1 = zeros( [ nip1 njp1 6 ] ); |
| 61 |
scubep1 = zeros( [ nip1 njp1 6 ] ); |
| 62 |
bcubep1 = zeros( [ nip1 njp1 6 ] ); |
| 63 |
for icf = 1:6 |
| 64 |
tcubep1(2:nip1,2:njp1,icf) = t_(:,:,icf); |
| 65 |
scubep1(2:nip1,2:njp1,icf) = s_(:,:,icf); |
| 66 |
bcubep1(2:nip1,2:njp1,icf) = b_(:,:,icf); |
| 67 |
end |
| 68 |
|
| 69 |
% Do the upwind-edge T exchanges |
| 70 |
tcubep1(1,2:njp1,1) = t_(ni:-1:1,nj,5); % - |
| 71 |
tcubep1(2:nip1,1,1) = t_(1:ni,nj,6); % - |
| 72 |
tcubep1(1,2:njp1,2) = t_(ni,1:nj,1); % - |
| 73 |
tcubep1(2:nip1,1,2) = t_(ni,nj:-1:1,6); % - |
| 74 |
tcubep1(1,2:njp1,3) = t_(ni:-1:1,nj,1); % - |
| 75 |
tcubep1(2:nip1,1,3) = t_(1:ni,nj,2); % - |
| 76 |
tcubep1(1,2:njp1,4) = t_(ni,1:nj,3); % - |
| 77 |
tcubep1(2:nip1,1,4) = t_(ni,nj:-1:1,2); % - |
| 78 |
tcubep1(1,2:njp1,5) = t_(ni:-1:1,nj,3); % - |
| 79 |
tcubep1(2:nip1,1,5) = t_(1:ni,nj,4); % - |
| 80 |
tcubep1(1,2:njp1,6) = t_(ni,1:nj,5); % - |
| 81 |
tcubep1(2:nip1,1,6) = t_(ni,nj:-1:1,4); % - |
| 82 |
|
| 83 |
% Do the upwind-edge S exchanges |
| 84 |
scubep1(1,2:njp1,1) = s_(ni:-1:1,nj,5); % - |
| 85 |
scubep1(2:nip1,1,1) = s_(1:ni,nj,6); % - |
| 86 |
scubep1(1,2:njp1,2) = s_(ni,1:nj,1); % - |
| 87 |
scubep1(2:nip1,1,2) = s_(ni,nj:-1:1,6); % - |
| 88 |
scubep1(1,2:njp1,3) = s_(ni:-1:1,nj,1); % - |
| 89 |
scubep1(2:nip1,1,3) = s_(1:ni,nj,2); % - |
| 90 |
scubep1(1,2:njp1,4) = s_(ni,1:nj,3); % - |
| 91 |
scubep1(2:nip1,1,4) = s_(ni,nj:-1:1,2); % - |
| 92 |
scubep1(1,2:njp1,5) = s_(ni:-1:1,nj,3); % - |
| 93 |
scubep1(2:nip1,1,5) = s_(1:ni,nj,4); % - |
| 94 |
scubep1(1,2:njp1,6) = s_(ni,1:nj,5); % - |
| 95 |
scubep1(2:nip1,1,6) = s_(ni,nj:-1:1,4); % - |
| 96 |
|
| 97 |
% Do the upwind-edge B exchanges |
| 98 |
bcubep1(1,2:njp1,1) = b_(ni:-1:1,nj,5); % - |
| 99 |
bcubep1(2:nip1,1,1) = b_(1:ni,nj,6); % - |
| 100 |
bcubep1(1,2:njp1,2) = b_(ni,1:nj,1); % - |
| 101 |
bcubep1(2:nip1,1,2) = b_(ni,nj:-1:1,6); % - |
| 102 |
bcubep1(1,2:njp1,3) = b_(ni:-1:1,nj,1); % - |
| 103 |
bcubep1(2:nip1,1,3) = b_(1:ni,nj,2); % - |
| 104 |
bcubep1(1,2:njp1,4) = b_(ni,1:nj,3); % - |
| 105 |
bcubep1(2:nip1,1,4) = b_(ni,nj:-1:1,2); % - |
| 106 |
bcubep1(1,2:njp1,5) = b_(ni:-1:1,nj,3); % - |
| 107 |
bcubep1(2:nip1,1,5) = b_(1:ni,nj,4); % - |
| 108 |
bcubep1(1,2:njp1,6) = b_(ni,1:nj,5); % - |
| 109 |
bcubep1(2:nip1,1,6) = b_(ni,nj:-1:1,4); % - |
| 110 |
|
| 111 |
% Get T values on the U,V grid points |
| 112 |
masku = 1.0 - abs(diff(double(tcubep1(2:nip1,:,:) == 0.0),1,2)); |
| 113 |
maskv = 1.0 - abs(diff(double(tcubep1(:,2:nip1,:) == 0.0),1,1)); |
| 114 |
diffu = 0.5*diff(tcubep1(2:nip1,:,:),1,2); |
| 115 |
diffv = 0.5*diff(tcubep1(:,2:nip1,:),1,1); |
| 116 |
tonu = t_(:,:,:) + masku.*diffu; |
| 117 |
tonv = t_(:,:,:) + maskv.*diffv; |
| 118 |
|
| 119 |
% Get S values on the U,V grid points |
| 120 |
diffu = 0.5*diff(scubep1(2:nip1,:,:),1,2); |
| 121 |
diffv = 0.5*diff(scubep1(:,2:nip1,:),1,1); |
| 122 |
sonu = s_(:,:,:) + masku.*diffu; |
| 123 |
sonv = s_(:,:,:) + maskv.*diffv; |
| 124 |
|
| 125 |
% Get B values on the U grid points |
| 126 |
diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2); |
| 127 |
diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1); |
| 128 |
bonu = b_(:,:,:) + masku.*diffu; |
| 129 |
bonv = b_(:,:,:) + maskv.*diffv; |
| 130 |
|
| 131 |
% Calculate the primes |
| 132 |
uptp = ut - u_.*tonu; |
| 133 |
vptp = vt - v_.*tonv; |
| 134 |
upsp = us - u_.*sonu; |
| 135 |
vpsp = vs - v_.*sonv; |
| 136 |
upbp = ub - u_.*bonu; |
| 137 |
vpbp = vb - v_.*bonv; |
| 138 |
|
| 139 |
% Calculate the linear-EoS approximation |
| 140 |
Alpha = 2.0*10^-4; |
| 141 |
Beta = 7.0*10^-4; |
| 142 |
lin_upbp = 9.81.*( Alpha.*uptp - Beta.*upsp ); |
| 143 |
lin_vpbp = 9.81.*( Alpha.*vptp - Beta.*vpsp ); |
| 144 |
|
| 145 |
lns_upbp = 9.81.*( Alpha.*uptp ); |
| 146 |
ljs_upbp = 9.81.*( -Beta.*upsp ); |
| 147 |
|
| 148 |
% Quick plots |
| 149 |
figure(1), surf( vptp(:,:,6)), view(2), shading interp, colorbar, ... |
| 150 |
title('v''T''') |
| 151 |
|
| 152 |
figure(1), surf( upbp(:,:,1)), view(2), shading interp, colorbar, ... |
| 153 |
title('u''B'' from JMD95 EoS') |
| 154 |
figure(2), surf(lin_upbp(:,:,1)), view(2), shading interp, colorbar, ... |
| 155 |
title('u''B'' from Linear EoS') |
| 156 |
figure(3), surf(lns_upbp(:,:,1)), view(2), shading interp, colorbar, ... |
| 157 |
title('u''B'' by forgetting salt') |
| 158 |
figure(4), surf(ljs_upbp(:,:,1)), view(2), shading interp, colorbar, ... |
| 159 |
title('u''B'' with just salt') |
| 160 |
|
| 161 |
|
| 162 |
figure(3), surf( uptp(:,:,1)), view(2), shading interp, colorbar, ... |
| 163 |
title('u''T''') |
| 164 |
figure(4), surf( upsp(:,:,1)), view(2), shading interp, colorbar, ... |
| 165 |
title('u''S''') |
| 166 |
|
| 167 |
figure(3), surf( s_(:,:,1)), view(2), shading interp, colorbar, ... |
| 168 |
title('S') |
| 169 |
figure(4), surf( t_(:,:,1)), view(2), shading interp, colorbar, ... |
| 170 |
title('T') |
| 171 |
|
| 172 |
figure(3), surf( sonu(:,:,1)), view(2), shading interp, colorbar, ... |
| 173 |
title('S (on U points)') |
| 174 |
figure(4), surf( tonu(:,:,1)), view(2), shading interp, colorbar, ... |
| 175 |
title('T (on U points)') |
| 176 |
|
| 177 |
figure(3), surf( ut(:,:,1)), view(2), shading interp, colorbar, ... |
| 178 |
title('UT') |
| 179 |
figure(4), surf( us(:,:,1)), view(2), shading interp, colorbar, ... |
| 180 |
title('US') |
| 181 |
|
| 182 |
% Check bar(B) with bar(T),bar(S) |
| 183 |
p = 0.09998 * 9.81 * R(iz); |
| 184 |
for icf = 1:6 |
| 185 |
bcheck(:,:,icf) = ... |
| 186 |
eh3densjmd95( squeeze(s_(:,:,icf)), ... |
| 187 |
squeeze(t_(:,:,icf)), p ) * -9.81 / 1000; |
| 188 |
end |
| 189 |
figure(3), surf( b_(:,:,1)), view(2), shading interp, ... |
| 190 |
caxis([-10.1 -9.8]), colorbar, title('bar(B) from model') |
| 191 |
figure(4), surf( bcheck(:,:,1)), view(2), shading interp, ... |
| 192 |
caxis([-10.1 -9.8]), colorbar, title('bar(B) from bar(T),bar(S)') |
| 193 |
|
| 194 |
|
| 195 |
|
| 196 |
end |
Parent Directory
| 
Revision Log
| 
Revision Graph