BoxAdj/input/gendata.m

clear all
close all
ieee='b';accuracy='single';

dz = [500 500 500 500 500 500 500 500 ]
zf=-cumsum([0 dz]);
zc=(zf(1:end-1)+zf(2:end))/2;

Ho=-zf(length(zf));
nx=64;ny=64;nz = 8;

if 1
    % set up input for model
    % Flat bottom at z=-Ho
    h=-Ho*ones(nx,ny);
    %h(end,:)=0;%h(:,end)=0;%WALLS
    fid=fopen('topog.box','w',ieee); fwrite(fid,h,accuracy); fclose(fid);
    
    T = [20.,16.,12.,10., 9., 8., 7., 6.]    
    T2D = single(ones(64,1)*T);       
%    plot((T2D(1,:)),zc,'-ro','linewidth',2)
    T3D = zeros(nx,ny,nz,'single');
    for i = 1:64;   T3D(i,:,:)= T2D;  end
  
  % T OBSERVED
    fid = fopen('FinalThetaObs.bin','w',ieee);
       fwrite(fid,T3D,accuracy);
    fclose(fid)

  %VPROFILE
    x = linspace(-1,1,ny)';%x = [-1 linspace(-1,1,ny-2) 1]
    V = single(1-(x.^2));
    V2D =single(V*ones(1,nz))./10;
  %V2D = zeros(64,nz,'single');
  for i = 1:64
    V3D(i,:,:) = V2D;
  end

% Lets make a zonal jet solution 
  %INITIAL CONDITION 
  fid = fopen('Vini.bin','w',ieee);
  fwrite(fid,V3D.*0,accuracy);
  fclose(fid)
  fid = fopen('Uini.bin','w',ieee);
  fwrite(fid,V3D,accuracy);
  fclose(fid)
  
  %NORTHERN BOUNDARY
  fid = fopen('Unbc.bin','w',ieee);
  fwrite(fid,zeros(nx,nz,accuracy),accuracy);
  fclose(fid)
  fid = fopen('Vnbc.bin','w',ieee);
  fwrite(fid,zeros(nx,nz,accuracy),accuracy);
  fclose(fid)
  fid = fopen('Snbc.bin','w',ieee);
  fwrite(fid,35.*ones(nx,nz,accuracy),accuracy);
  fclose(fid)
  fid = fopen('Tnbc.bin','w',ieee);
  fwrite(fid,T2D,accuracy);
  fclose(fid) 

  %SOUTHERN BOUNDARY
  fid = fopen('Usbc.bin','w',ieee);
  fwrite(fid,zeros(nx,nz,accuracy),accuracy);
  fclose(fid)
  fid = fopen('Vsbc.bin','w',ieee);
  fwrite(fid,zeros(nx,nz,accuracy),accuracy);
  fclose(fid)
  fid = fopen('Ssbc.bin','w',ieee);
  fwrite(fid,35.*ones(nx,nz,accuracy),accuracy);
  fclose(fid)
  fid = fopen('Tsbc.bin','w',ieee);
  fwrite(fid,T2D,accuracy);
  fclose(fid)

  %WESTERN BOUNDARY
  fid = fopen('Uwbc.bin','w',ieee);
  fwrite(fid,V2D,accuracy);
  fclose(fid)
  fid = fopen('Vwbc.bin','w',ieee);
  fwrite(fid,zeros(ny,nz,accuracy),accuracy);
  fclose(fid)
  fid = fopen('Swbc.bin','w',ieee);
  fwrite(fid,35.*ones(ny,nz,accuracy),accuracy);
  fclose(fid)
  fid = fopen('Twbc.bin','w',ieee);
  fwrite(fid,T2D,accuracy);
  fclose(fid)


  %EASTERN BOUNDARY
  fid = fopen('Uebc.bin','w',ieee);
  fwrite(fid,V2D,accuracy);
  fclose(fid)
  fid = fopen('Vebc.bin','w',ieee);
  fwrite(fid,zeros(ny,nz,accuracy),accuracy);
  fclose(fid)
  fid = fopen('Sebc.bin','w',ieee);
  fwrite(fid,35.*ones(ny,nz,accuracy),accuracy);
  fclose(fid)
  fid = fopen('Tebc.bin','w',ieee);
  fwrite(fid,T2D,accuracy);
  fclose(fid)
  % end of model setup if-then
end


if 0
  % get vertical modes for obcs decomp
  rhonil = 1035;g = 9.81;
  RC = squeeze(rdmds('../GRID/RC'));
  RF = squeeze(rdmds('../GRID/RF'));  
  DRF = squeeze(rdmds('../GRID/DRF'));  
  mask = squeeze(rdmds('../GRID/maskInC'));
  RLOW=RF(end);
  zmid = [0; (RC(1:end-1)+RC(2:end))/2;RLOW];
  NZ = length(RC);
  
  % all the dRhdz* come down to Nz
  dRhodz_bar = rdmds('../run_fwd/dRhodz_5');
  dRhodz_bar(mask==0)=0;dRhodz_bar(dRhodz_bar==0) = nan;
  dRhodz_bar = squeeze(nanmean(nanmean(dRhodz_bar(32:33,32:33,:),1),2));
  dRhodz_bar(1)=dRhodz_bar(2);dRhodz_bar(NZ+1)=dRhodz_bar(NZ);
  Nz = (-g./rhonil.*dRhodz_bar).^0.5;
  
  modesv = zeros(NZ,NZ,NZ);
  % when only one depth, just have 1 for the mode
  modesv(1,1,1) = 1.0;
  % for more than 1 depth, solve eigenvalue problem
  for k = 2:NZ;
    % iz = vector of depth levels
    iz = 1:k;
    % print out a progress number
    NZ-k,
    % regularly spaced depths to the bottom of layer k;
    % leave out the surface: VERT_FSFB2.m will supply it.
    % 5 m spacing was selected as being small enough
    zreg=-5:-5:RF(k+1);
    Nzreg = interp1(zmid,Nz.^2,zreg,'linear',0);
    [c2, F, G, N2, Pmid] = VERT_FSFB2(Nzreg,zreg);
    %VERT_FSFB2.m adds a surface point
    zreg = [0,zreg];
    % sort from largest to smallest phase speed
    [c2s indx] = sort(abs(c2),'descend');
    % now interp back to our grid
    for mds = 1:k
      YI = interp1(zreg,F(:,indx(mds)),RC(iz),'linear',0);
      modesv(iz,mds,k) = YI;
    end
    % have to weight by dz!  Make a vector of fractional dz's
    zwt = DRF(iz)/sum(DRF(iz),1);
    %ensure first mode is barotropic (constant in depth)
    avm1=sum(modesv(iz,1,k).*zwt,1);
    modesv(iz,1,k)=avm1;
    %make all modes orthogonal weighted by delta z
    % use gram-schmidt leaving first one the same
    for mds = 1:k-1
        R = sum((modesv(iz,mds,k).^2).*zwt,1);
        R2 = (modesv(iz,mds,k).*zwt)'*modesv(iz,mds+1:k,k)/R;
        modesv(iz,mds+1:k,k) = modesv(iz,mds+1:k,k) - ...
            modesv(iz,mds,k)*R2;
    end
    %All now orthognal, now normalize
    for mds = 1:k
       R = sqrt(sum((modesv(iz,mds,k).^2).*zwt,1));
       if R < 1e-8
         fprintf('Small R!! for mds = %2i and k = %2i\n',mds,k);
         R = inf;
       end
       modesv(iz,mds,k) = modesv(iz,mds,k)./R;
     end
  end;% end of loop over depth level k
  fid = fopen('baro_invmodes.bin','w','b');
  fwrite(fid,modesv,'double');
  fclose(fid)
  if 1 %plot first 5 modes for deepest case
    figure
    clf;plot(modesv(:,[1:5],NZ),RC(:));
    title('output modes at deepest point')
  end
  if 1  % test orthogonality
    % do whole matrix (need to include dz!)
    k = NZ;
    cmm = (squeeze(modesv(iz,iz,k)).*repmat(zwt,[1 k]))'* ...
          squeeze(modesv(iz,iz,k));
    figure;imagesc(cmm);colorbar
    title([num2str(k) ' mode orthogonality min, max diag: ' ...
           num2str(min(diag(cmm))) ', ' ...
           num2str(max(diag(cmm)))])
  end
  save baro_invmodes modesv
end

1	mmazloff	1.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