1 |
clear all |
2 |
close all |
3 |
ieee='b';accuracy='single'; |
4 |
|
5 |
dz = [500 500 500 500 500 500 500 500 ] |
6 |
zf=-cumsum([0 dz]); |
7 |
zc=(zf(1:end-1)+zf(2:end))/2; |
8 |
|
9 |
Ho=-zf(length(zf)); |
10 |
nx=64;ny=64;nz = 8; |
11 |
|
12 |
if 1 |
13 |
% set up input for model |
14 |
% Flat bottom at z=-Ho |
15 |
h=-Ho*ones(nx,ny); |
16 |
%h(end,:)=0;%h(:,end)=0;%WALLS |
17 |
fid=fopen('topog.box','w',ieee); fwrite(fid,h,accuracy); fclose(fid); |
18 |
|
19 |
T = [20.,16.,12.,10., 9., 8., 7., 6.] |
20 |
T2D = single(ones(64,1)*T); |
21 |
% plot((T2D(1,:)),zc,'-ro','linewidth',2) |
22 |
T3D = zeros(nx,ny,nz,'single'); |
23 |
for i = 1:64; T3D(i,:,:)= T2D; end |
24 |
|
25 |
% T OBSERVED |
26 |
fid = fopen('FinalThetaObs.bin','w',ieee); |
27 |
fwrite(fid,T3D,accuracy); |
28 |
fclose(fid) |
29 |
|
30 |
%VPROFILE |
31 |
x = linspace(-1,1,ny)';%x = [-1 linspace(-1,1,ny-2) 1] |
32 |
V = single(1-(x.^2)); |
33 |
V2D =single(V*ones(1,nz))./10; |
34 |
%V2D = zeros(64,nz,'single'); |
35 |
for i = 1:64 |
36 |
V3D(i,:,:) = V2D; |
37 |
end |
38 |
|
39 |
% Lets make a zonal jet solution |
40 |
%INITIAL CONDITION |
41 |
fid = fopen('Vini.bin','w',ieee); |
42 |
fwrite(fid,V3D.*0,accuracy); |
43 |
fclose(fid) |
44 |
fid = fopen('Uini.bin','w',ieee); |
45 |
fwrite(fid,V3D,accuracy); |
46 |
fclose(fid) |
47 |
|
48 |
%NORTHERN BOUNDARY |
49 |
fid = fopen('Unbc.bin','w',ieee); |
50 |
fwrite(fid,zeros(nx,nz,accuracy),accuracy); |
51 |
fclose(fid) |
52 |
fid = fopen('Vnbc.bin','w',ieee); |
53 |
fwrite(fid,zeros(nx,nz,accuracy),accuracy); |
54 |
fclose(fid) |
55 |
fid = fopen('Snbc.bin','w',ieee); |
56 |
fwrite(fid,35.*ones(nx,nz,accuracy),accuracy); |
57 |
fclose(fid) |
58 |
fid = fopen('Tnbc.bin','w',ieee); |
59 |
fwrite(fid,T2D,accuracy); |
60 |
fclose(fid) |
61 |
|
62 |
%SOUTHERN BOUNDARY |
63 |
fid = fopen('Usbc.bin','w',ieee); |
64 |
fwrite(fid,zeros(nx,nz,accuracy),accuracy); |
65 |
fclose(fid) |
66 |
fid = fopen('Vsbc.bin','w',ieee); |
67 |
fwrite(fid,zeros(nx,nz,accuracy),accuracy); |
68 |
fclose(fid) |
69 |
fid = fopen('Ssbc.bin','w',ieee); |
70 |
fwrite(fid,35.*ones(nx,nz,accuracy),accuracy); |
71 |
fclose(fid) |
72 |
fid = fopen('Tsbc.bin','w',ieee); |
73 |
fwrite(fid,T2D,accuracy); |
74 |
fclose(fid) |
75 |
|
76 |
%WESTERN BOUNDARY |
77 |
fid = fopen('Uwbc.bin','w',ieee); |
78 |
fwrite(fid,V2D,accuracy); |
79 |
fclose(fid) |
80 |
fid = fopen('Vwbc.bin','w',ieee); |
81 |
fwrite(fid,zeros(ny,nz,accuracy),accuracy); |
82 |
fclose(fid) |
83 |
fid = fopen('Swbc.bin','w',ieee); |
84 |
fwrite(fid,35.*ones(ny,nz,accuracy),accuracy); |
85 |
fclose(fid) |
86 |
fid = fopen('Twbc.bin','w',ieee); |
87 |
fwrite(fid,T2D,accuracy); |
88 |
fclose(fid) |
89 |
|
90 |
|
91 |
%EASTERN BOUNDARY |
92 |
fid = fopen('Uebc.bin','w',ieee); |
93 |
fwrite(fid,V2D,accuracy); |
94 |
fclose(fid) |
95 |
fid = fopen('Vebc.bin','w',ieee); |
96 |
fwrite(fid,zeros(ny,nz,accuracy),accuracy); |
97 |
fclose(fid) |
98 |
fid = fopen('Sebc.bin','w',ieee); |
99 |
fwrite(fid,35.*ones(ny,nz,accuracy),accuracy); |
100 |
fclose(fid) |
101 |
fid = fopen('Tebc.bin','w',ieee); |
102 |
fwrite(fid,T2D,accuracy); |
103 |
fclose(fid) |
104 |
% end of model setup if-then |
105 |
end |
106 |
|
107 |
|
108 |
if 0 |
109 |
% get vertical modes for obcs decomp |
110 |
rhonil = 1035;g = 9.81; |
111 |
RC = squeeze(rdmds('../GRID/RC')); |
112 |
RF = squeeze(rdmds('../GRID/RF')); |
113 |
DRF = squeeze(rdmds('../GRID/DRF')); |
114 |
mask = squeeze(rdmds('../GRID/maskInC')); |
115 |
RLOW=RF(end); |
116 |
zmid = [0; (RC(1:end-1)+RC(2:end))/2;RLOW]; |
117 |
NZ = length(RC); |
118 |
|
119 |
% all the dRhdz* come down to Nz |
120 |
dRhodz_bar = rdmds('../run_fwd/dRhodz_5'); |
121 |
dRhodz_bar(mask==0)=0;dRhodz_bar(dRhodz_bar==0) = nan; |
122 |
dRhodz_bar = squeeze(nanmean(nanmean(dRhodz_bar(32:33,32:33,:),1),2)); |
123 |
dRhodz_bar(1)=dRhodz_bar(2);dRhodz_bar(NZ+1)=dRhodz_bar(NZ); |
124 |
Nz = (-g./rhonil.*dRhodz_bar).^0.5; |
125 |
|
126 |
modesv = zeros(NZ,NZ,NZ); |
127 |
% when only one depth, just have 1 for the mode |
128 |
modesv(1,1,1) = 1.0; |
129 |
% for more than 1 depth, solve eigenvalue problem |
130 |
for k = 2:NZ; |
131 |
% iz = vector of depth levels |
132 |
iz = 1:k; |
133 |
% print out a progress number |
134 |
NZ-k, |
135 |
% regularly spaced depths to the bottom of layer k; |
136 |
% leave out the surface: VERT_FSFB2.m will supply it. |
137 |
% 5 m spacing was selected as being small enough |
138 |
zreg=-5:-5:RF(k+1); |
139 |
Nzreg = interp1(zmid,Nz.^2,zreg,'linear',0); |
140 |
[c2, F, G, N2, Pmid] = VERT_FSFB2(Nzreg,zreg); |
141 |
%VERT_FSFB2.m adds a surface point |
142 |
zreg = [0,zreg]; |
143 |
% sort from largest to smallest phase speed |
144 |
[c2s indx] = sort(abs(c2),'descend'); |
145 |
% now interp back to our grid |
146 |
for mds = 1:k |
147 |
YI = interp1(zreg,F(:,indx(mds)),RC(iz),'linear',0); |
148 |
modesv(iz,mds,k) = YI; |
149 |
end |
150 |
% have to weight by dz! Make a vector of fractional dz's |
151 |
zwt = DRF(iz)/sum(DRF(iz),1); |
152 |
%ensure first mode is barotropic (constant in depth) |
153 |
avm1=sum(modesv(iz,1,k).*zwt,1); |
154 |
modesv(iz,1,k)=avm1; |
155 |
%make all modes orthogonal weighted by delta z |
156 |
% use gram-schmidt leaving first one the same |
157 |
for mds = 1:k-1 |
158 |
R = sum((modesv(iz,mds,k).^2).*zwt,1); |
159 |
R2 = (modesv(iz,mds,k).*zwt)'*modesv(iz,mds+1:k,k)/R; |
160 |
modesv(iz,mds+1:k,k) = modesv(iz,mds+1:k,k) - ... |
161 |
modesv(iz,mds,k)*R2; |
162 |
end |
163 |
%All now orthognal, now normalize |
164 |
for mds = 1:k |
165 |
R = sqrt(sum((modesv(iz,mds,k).^2).*zwt,1)); |
166 |
if R < 1e-8 |
167 |
fprintf('Small R!! for mds = %2i and k = %2i\n',mds,k); |
168 |
R = inf; |
169 |
end |
170 |
modesv(iz,mds,k) = modesv(iz,mds,k)./R; |
171 |
end |
172 |
end;% end of loop over depth level k |
173 |
fid = fopen('baro_invmodes.bin','w','b'); |
174 |
fwrite(fid,modesv,'double'); |
175 |
fclose(fid) |
176 |
if 1 %plot first 5 modes for deepest case |
177 |
figure |
178 |
clf;plot(modesv(:,[1:5],NZ),RC(:)); |
179 |
title('output modes at deepest point') |
180 |
end |
181 |
if 1 % test orthogonality |
182 |
% do whole matrix (need to include dz!) |
183 |
k = NZ; |
184 |
cmm = (squeeze(modesv(iz,iz,k)).*repmat(zwt,[1 k]))'* ... |
185 |
squeeze(modesv(iz,iz,k)); |
186 |
figure;imagesc(cmm);colorbar |
187 |
title([num2str(k) ' mode orthogonality min, max diag: ' ... |
188 |
num2str(min(diag(cmm))) ', ' ... |
189 |
num2str(max(diag(cmm)))]) |
190 |
end |
191 |
save baro_invmodes modesv |
192 |
end |
193 |
|