| 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 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph