MLAdjust/input.1.leith/genit.m

%function geninit()
%  This function generates initial condition files for the mitgcm.
% if init_vel=1, then the front problem is used
% if init_vel=0, then resting (other than small symmetry-breaking
% random noise) initial conditions are used.

nx=50
ny=nx/2+1
nz=40
init_vel=1;
dxspacing=1e3
dyspacing=dxspacing

Lx=dxspacing*nx
Ly=dyspacing*ny 

%-- Params
g=9.81;
tAlpha=-2e-4;
f0=7.29e-5;
rho=1035
day=24*60^2;
time_scale=Inf %day/24/4
surf_time_scale=Inf
tau0=0.0           % Wind forcing magnitude...
daylength=day/3;    % Length of daytime heating
qconst=0          % Surface steady heat flux W/m^2 (positive for ocean cool)

H=2000;            %Max depth

%-- Grid: x
dx_ratio=30; dx_trans=.05; dx_min=100;
xn=(0.5:nx)/nx;
dx=ones(1,nx);                                  % uniform resolution
dx=dx*Lx/sum(dx); 
xf=cumsum([0 dx]); % Face x points
xc=(xf(1:end-1)+xf(2:end))/2; % Centered x points

%-- Grid: y
dy_ratio=30; dy_trans=.02; dy_min=100;
yn=(0.5:ny)/ny;
dy=ones(1,ny);                                  % uniform resolution
dy=dy*Ly/sum(dy); 
yf=cumsum([0 dy]);  % Face y-points
yc=(yf(1:end-1)+yf(2:end))/2;  %Centered y-points
L=yc(end)-yc(1);        % this takes into account the wall of topography!!!

%-- Grid: z

dh=H/nz*ones(1,nz);
%for i=21:38
%  dh(i)=dh(i-1)*1.1917;
%end

zf=-round(cumsum([0 dh])/1)*1;   % Face z points
dh=-diff(zf);
zc=(zf(1:end-1)+zf(2:end))/2;  % centered z points
nz=length(dh);
H=sum(dh)

[XT,YT,ZT]=ndgrid(xc,yc,zc); % This is the centered, temperature grid.
[XU,YU,ZU]=ndgrid(xf,yc,zc); % This is the U grid.
[XV,YV,ZV]=ndgrid(xc,yf,zc); % This is the V grid.
[XW,YW,ZW]=ndgrid(xc,yc,zf); % This is the W grid.
[XB,YB]=ndgrid(xc,yc); % This is the Bathymetry grid.

% Stratification Params
Dml=200;             % Mixed-layer Depth
D=Dml;
Lf=64*1e3;          % Half-Width of Front
%Lf=50*dyspacing;          % Half-Width of Front
y0=yc(round(length(yc)/2));     % centered location of Front


Ms=(2*f0)^2     %db/dy  2e-4 is about 1 K/50 km
Nsml=0*(f0)^2     %db/dz
Nsint=(32*f0)^2      %db/dz in interior

Ri=Nsml*f0^2/Ms^2;
Ld=max(Nsml^.5/f0*Dml, Ms/f0^2*Dml)
Ls=(f0^2/Ms/Dml*sqrt(5/2/(1+max(Ri,1))))^-1*2*pi

Riint=Nsint*f0^2/Ms^2;
Ldint=max(Nsint^.5/f0*H, Ms/f0^2*H)
Lsint=(f0^2/Ms/H*sqrt(5/2/(1+max(Riint,1))))^-1*2*pi

Resfacs=[Lf/Ld Lf/Ls Ls/dxspacing Ld/dxspacing]

dxspacing
optdxspacing=round(Ls*0.05)

deltatheta=2*Lf*Ms/g/tAlpha
T0=17
dthetadzml=-Nsml/g/tAlpha
dthetadzint=-Nsint/g/tAlpha

%-- Bathymetry
% Bathymetry is on Xc, Yc
hh=ones(size(XB));

% Variable Bathymetry
% hh(:,:)=hh(:,:)-2*max(0,abs(YB/L-1/2)-3/8);

hh(:,end)=0*hh(:,end);


hh=-H*hh;

Ltop=L/2;

figure(1)
subplot(221)
pcolor(XB/1e3,YB/1e3,hh)


subplot(223)
plot(xc/1e3,dx/1e3,'.'); xlabel('X (km)');ylabel('\Delta x (km)')
title('\Delta x')

subplot(222)
plot(dy/1e3,yc/1e3,'.'); ylabel('Y (km)');xlabel('\Delta y (km)')
title('\Delta y')

subplot(224)
plot(dh,zc,'.'); ylabel('Z (m)');xlabel('\Delta z (m)')
title('\Delta z')

% Surface maximum diurnal heat flux (gets multiplied by cos(2pi*t/day)

figure(3)
deltaT=3600;
t=deltaT:deltaT:86400*2;

q1=sum(max(cos(2*pi*t/day),cos(pi*daylength/day))-cos(pi*daylength/day))/length(t);

qdiur=-qconst/q1;

q1=qconst+qdiur*(max(cos(2*pi*t/day),cos(pi*daylength/day))-cos(pi*daylength/day));

plot(t,q1)

sprintf('Set data param fuconst to %0.5g',tau0)
sprintf('Set data param daylength to %0.5g',daylength)
sprintf('Set data param qdiur to %0.5g',qdiur)
sprintf('Set data param qconst to %0.5g',qconst)

xlabel('t (days)'); ylabel('Q Flux (pos=ocean cool)');title('Heat flux')

fid=fopen('dx_variable.bin','w','b'); fwrite(fid,dx,'real*8'); fclose(fid);
fid=fopen('dy_variable.bin','w','b'); fwrite(fid,dy,'real*8'); fclose(fid);
fid=fopen('topo_sl.bin','w','b'); fwrite(fid,hh,'real*8'); fclose(fid);

% Initial Temp Profile

figure(2)
subplot(4,2,1)

theta=T0+dthetadzint*(ZT+Dml)+deltatheta*tanh((YT-y0)/Lf)/2;

% Impose a strong Dml
i=1;
iml=1;

while ((ZT(1,1,iml)>-Dml)&(iml<length(ZT(1,1,:))))
  iml=iml+1;
end

while (i<iml)
  theta(:,:,i)=T0+dthetadzml*(ZT(:,:,i)+Dml)+deltatheta*tanh((YT(:,:,i)-y0)/Lf)/2;
  i=i+1;
end

[ZT(1,1,1) theta(1,1,1) theta(1,end,1) ;ZT(1,1,i) theta(1,1,i) theta(1,end,i);ZT(1,1,end) theta(1,1,end) theta(1,end,end)]

Tvec=15:.1:18;

subplot(4,2,1)
[h,c]=contourf(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(theta(:,1,:)),Tvec);
title('Potl Temp')
xlabel('x (km)');ylabel('z (m)')
fid=fopen('thetaRestoreFile.bin','w','b'); fwrite(fid,theta,'real*8'); fclose(fid);

subplot(4,2,2)
[h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)),Tvec);
colorbar
title('Potl Temp')
xlabel('y (km)');ylabel('z (m)')
fid=fopen('thetaRestoreFile.bin','w','b'); fwrite(fid,theta,'real*8'); fclose(fid);

Rvec=1031:.05:1032;

subplot(4,2,3)
dens=theta*rho*tAlpha+rho;
[h,c]=contour(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(dens(:,1,:)),Rvec);
title('Density')
xlabel('x (km)');ylabel('z (m)')

subplot(4,2,4)
[h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)),Rvec);
title('Density')
colorbar
xlabel('y (km)');ylabel('z (m)')

%Spice

x1=xc(round(length(xc)/4));     % centered location of Front #1
x2=xc(round(3*length(xc)/4));     % centered location of Front #2

%spice=((SL+SR)+(SR-SL)*tanh((ZT-Dml/2)/Dml))/2

spice=T0+dthetadzint*(ZT+Dml)+deltatheta*(tanh((XT-x1)/Lf)-tanh((XT-x2)/Lf)-1)/2;

i=1

while (i<iml)
  spice(:,:,i)=T0+dthetadzml*(ZT(:,:,i)+Dml)+deltatheta*(tanh((XT(:,:,i)-x1)/Lf)-tanh((XT(:,:,i)-x2)/Lf)-1)/2;
  i=i+1;
end

fid=fopen('spiceRestoreFile.bin','w','b'); fwrite(fid,spice,'real*8'); fclose(fid);

subplot(4,2,8)
[h,c]=contourf(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(spice(:,1,:)),Tvec);
%clabel(h,c)
title('Spice')
xlabel('x (km)');ylabel('z (m)')

%clabel(h,c);

subplot(4,2,7)
invtime=0*YT;
decay=dyspacing;
invtime=1/time_scale*(exp(-YT/decay+YT(1,1,1)/decay)+exp(YT/decay-YT(1,end-1,1)/decay));
% Surface Restoring
invtime(:,:,1)=invtime(:,:,1)+1/surf_time_scale;

figure(4)
subplot(2,1,1)
plot(squeeze(YT(1,1:end-1,:)/1e3),squeeze(invtime(1,1:end-1,:)*86400))
ylabel('inverse restoring time (1/day)') 
xlabel('Y (km)')

subplot(2,1,2)
plot(squeeze(invtime(1,1:end-1,:)*86400)',squeeze(ZT(1,1:end-1,:))')
xlabel('inverse restoring time (1/day)') 
ylabel('depth')

figure(2)
plot(YB(1,:)/1e3,hh(1,:))
title('topography')
xlabel('x (km)');ylabel('value')

fid=fopen('invtimeRestoreFile.bin','w','b'); fwrite(fid,invtime,'real*8'); fclose(fid);

%-- Perturb Initial temperature

pert=rand(size(theta(:,:,1)));

pert=1e-4*(pert-0.5);

for i=1:length(theta(1,1,:))
  theta(:,:,i)=theta(:,:,i)+pert;
end

fid=fopen('thetaInitial.bin','w','b');
fwrite(fid,theta,'real*8'); fclose(fid);

fid=fopen('spiceInitial.bin','w','b');
fwrite(fid,spice,'real*8'); fclose(fid);

%-- Thermal wind
U = 0*XU;

for k=(nz-1):-1:1
  U(1:end-1,1:end-1,k) = U(1:end-1,1:end-1,k+1)+dh(k)*g/f0/rho*diff(dens(:,:,k),1,2)./diff(YU(1:end-1,:,k+1),1,2);
  U(:,end,k) = 0;      % Get endpoint
  U(end,:,k) = U(1,:,k);             % In the wall
end

% Make initial state Momentumless (at least boussinesq, good approx)

[imax,jmax,kmax]=size(U);

for i=1:imax
 for j=1:jmax
  U(i,j,:)=U(i,j,:)-mean(U(i,j,:));
 end
end

timestep=dxspacing/max(abs(U(:)))/3.

subplot(4,2,5)
[h,c]=contour(squeeze(XU(:,1,:))/1e3,squeeze(ZU(:,1,:)),squeeze(U(:,1,:)),-1:.02:1);
title('U');xlabel('x (km)');ylabel('z (m)')

subplot(4,2,6)
[h,c]=contour(squeeze(YU(1,:,:))/1e3,squeeze(ZU(1,:,:)),squeeze(U(1,:,:)),-1:.02:1);

title('U');xlabel('x (km)');ylabel('z (m)')

U2=(U([1:end-1],:,:)+U([2:end],:,:))/2;

fid=fopen('uInitial.bin','w','b');
fwrite(fid,U2,'real*8'); fclose(fid);


%-- Wind stress

%tau=tau0*tanh((((squeeze(YU(:,:,1))-y0-Lf))/Lf));

% Constant Wind Stress
%tau=tau0*ones(size(squeeze(YU(:,:,1))));

%subplot(4,2,8)
%quiver(squeeze(XU(:,:,1))/1e3,squeeze(YU(:,:,1))/1e3,tau,0*tau)
%title('Taux');xlabel('x (km)');ylabel('y (km)')
%fid=fopen('taux_sl.bin','w','b');
%fwrite(fid,tau,'real*8'); fclose(fid);

%save 'zf' zf

1	baylor	1.1	%function geninit()
2			% This function generates initial condition files for the mitgcm.
3			% if init_vel=1, then the front problem is used
4			% if init_vel=0, then resting (other than small symmetry-breaking
5			% random noise) initial conditions are used.
6
7			nx=50
8			ny=nx/2+1
9			nz=40
10			init_vel=1;
11	baylor	1.2	dxspacing=1e3
12	baylor	1.1	dyspacing=dxspacing
13
14			Lx=dxspacing*nx
15	baylor	1.2	Ly=dyspacing*ny
16	baylor	1.1
17			%-- Params
18			g=9.81;
19			tAlpha=-2e-4;
20			f0=7.29e-5;
21			rho=1035
22			day=24*60^2;
23	baylor	1.2	time_scale=Inf %day/24/4
24			surf_time_scale=Inf
25			tau0=0.0 % Wind forcing magnitude...
26			daylength=day/3; % Length of daytime heating
27			qconst=0 % Surface steady heat flux W/m^2 (positive for ocean cool)
28	baylor	1.1
29	baylor	1.2	H=2000; %Max depth
30	baylor	1.1
31			%-- Grid: x
32	baylor	1.2	dx_ratio=30; dx_trans=.05; dx_min=100;
33			xn=(0.5:nx)/nx;
34	baylor	1.1	dx=ones(1,nx); % uniform resolution
35			dx=dx*Lx/sum(dx);
36			xf=cumsum([0 dx]); % Face x points
37			xc=(xf(1:end-1)+xf(2:end))/2; % Centered x points
38
39			%-- Grid: y
40	baylor	1.2	dy_ratio=30; dy_trans=.02; dy_min=100;
41			yn=(0.5:ny)/ny;
42	baylor	1.1	dy=ones(1,ny); % uniform resolution
43			dy=dy*Ly/sum(dy);
44			yf=cumsum([0 dy]); % Face y-points
45			yc=(yf(1:end-1)+yf(2:end))/2; %Centered y-points
46			L=yc(end)-yc(1); % this takes into account the wall of topography!!!
47
48			%-- Grid: z
49	baylor	1.2
50	baylor	1.1	dh=H/nz*ones(1,nz);
51	baylor	1.2	%for i=21:38
52			% dh(i)=dh(i-1)*1.1917;
53			%end
54
55	baylor	1.1	zf=-round(cumsum([0 dh])/1)*1; % Face z points
56			dh=-diff(zf);
57			zc=(zf(1:end-1)+zf(2:end))/2; % centered z points
58			nz=length(dh);
59			H=sum(dh)
60
61			[XT,YT,ZT]=ndgrid(xc,yc,zc); % This is the centered, temperature grid.
62			[XU,YU,ZU]=ndgrid(xf,yc,zc); % This is the U grid.
63			[XV,YV,ZV]=ndgrid(xc,yf,zc); % This is the V grid.
64			[XW,YW,ZW]=ndgrid(xc,yc,zf); % This is the W grid.
65			[XB,YB]=ndgrid(xc,yc); % This is the Bathymetry grid.
66
67			% Stratification Params
68			Dml=200; % Mixed-layer Depth
69			D=Dml;
70	baylor	1.2	Lf=64*1e3; % Half-Width of Front
71			%Lf=50*dyspacing; % Half-Width of Front
72	baylor	1.1	y0=yc(round(length(yc)/2)); % centered location of Front
73
74
75	baylor	1.2	Ms=(2*f0)^2 %db/dy 2e-4 is about 1 K/50 km
76			Nsml=0*(f0)^2 %db/dz
77			Nsint=(32*f0)^2 %db/dz in interior
78
79			Ri=Nsml*f0^2/Ms^2;
80			Ld=max(Nsml^.5/f0Dml, Ms/f0^2Dml)
81			Ls=(f0^2/Ms/Dmlsqrt(5/2/(1+max(Ri,1))))^-12*pi
82
83			Riint=Nsint*f0^2/Ms^2;
84			Ldint=max(Nsint^.5/f0H, Ms/f0^2H)
85			Lsint=(f0^2/Ms/Hsqrt(5/2/(1+max(Riint,1))))^-12*pi
86
87			Resfacs=[Lf/Ld Lf/Ls Ls/dxspacing Ld/dxspacing]
88
89			dxspacing
90			optdxspacing=round(Ls*0.05)
91
92			deltatheta=2LfMs/g/tAlpha
93	baylor	1.1	T0=17
94	baylor	1.2	dthetadzml=-Nsml/g/tAlpha
95			dthetadzint=-Nsint/g/tAlpha
96	baylor	1.1
97			%-- Bathymetry
98			% Bathymetry is on Xc, Yc
99			hh=ones(size(XB));
100	baylor	1.2
101			% Variable Bathymetry
102			% hh(:,:)=hh(:,:)-2*max(0,abs(YB/L-1/2)-3/8);
103
104	baylor	1.1	hh(:,end)=0*hh(:,end);
105	baylor	1.2
106
107	baylor	1.1	hh=-H*hh;
108
109	baylor	1.2	Ltop=L/2;
110
111	baylor	1.1	figure(1)
112			subplot(221)
113			pcolor(XB/1e3,YB/1e3,hh)
114	baylor	1.2
115	baylor	1.1
116			subplot(223)
117			plot(xc/1e3,dx/1e3,'.'); xlabel('X (km)');ylabel('\Delta x (km)')
118			title('\Delta x')
119
120			subplot(222)
121			plot(dy/1e3,yc/1e3,'.'); ylabel('Y (km)');xlabel('\Delta y (km)')
122			title('\Delta y')
123
124			subplot(224)
125			plot(dh,zc,'.'); ylabel('Z (m)');xlabel('\Delta z (m)')
126			title('\Delta z')
127
128	baylor	1.2	% Surface maximum diurnal heat flux (gets multiplied by cos(2pi*t/day)
129
130			figure(3)
131			deltaT=3600;
132			t=deltaT:deltaT:86400*2;
133
134			q1=sum(max(cos(2pit/day),cos(pidaylength/day))-cos(pidaylength/day))/length(t);
135
136			qdiur=-qconst/q1;
137
138			q1=qconst+qdiur(max(cos(2pit/day),cos(pidaylength/day))-cos(pi*daylength/day));
139
140			plot(t,q1)
141
142			sprintf('Set data param fuconst to %0.5g',tau0)
143			sprintf('Set data param daylength to %0.5g',daylength)
144			sprintf('Set data param qdiur to %0.5g',qdiur)
145			sprintf('Set data param qconst to %0.5g',qconst)
146
147			xlabel('t (days)'); ylabel('Q Flux (pos=ocean cool)');title('Heat flux')
148
149			fid=fopen('dx_variable.bin','w','b'); fwrite(fid,dx,'real*8'); fclose(fid);
150			fid=fopen('dy_variable.bin','w','b'); fwrite(fid,dy,'real*8'); fclose(fid);
151			fid=fopen('topo_sl.bin','w','b'); fwrite(fid,hh,'real*8'); fclose(fid);
152	baylor	1.1
153			% Initial Temp Profile
154	baylor	1.2
155	baylor	1.1	figure(2)
156			subplot(4,2,1)
157
158	baylor	1.2	theta=T0+dthetadzint(ZT+Dml)+deltathetatanh((YT-y0)/Lf)/2;
159
160			% Impose a strong Dml
161	baylor	1.1	i=1;
162			iml=1;
163
164			while ((ZT(1,1,iml)>-Dml)&(iml<length(ZT(1,1,:))))
165			iml=iml+1;
166			end
167
168			while (i<iml)
169	baylor	1.2	theta(:,:,i)=T0+dthetadzml(ZT(:,:,i)+Dml)+deltathetatanh((YT(:,:,i)-y0)/Lf)/2;
170	baylor	1.1	i=i+1;
171			end
172
173			[ZT(1,1,1) theta(1,1,1) theta(1,end,1) ;ZT(1,1,i) theta(1,1,i) theta(1,end,i);ZT(1,1,end) theta(1,1,end) theta(1,end,end)]
174
175	baylor	1.2	Tvec=15:.1:18;
176	baylor	1.1
177			subplot(4,2,1)
178	baylor	1.2	[h,c]=contourf(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(theta(:,1,:)),Tvec);
179	baylor	1.1	title('Potl Temp')
180	baylor	1.2	xlabel('x (km)');ylabel('z (m)')
181			fid=fopen('thetaRestoreFile.bin','w','b'); fwrite(fid,theta,'real*8'); fclose(fid);
182	baylor	1.1
183			subplot(4,2,2)
184	baylor	1.2	[h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)),Tvec);
185	baylor	1.1	colorbar
186			title('Potl Temp')
187			xlabel('y (km)');ylabel('z (m)')
188	baylor	1.2	fid=fopen('thetaRestoreFile.bin','w','b'); fwrite(fid,theta,'real*8'); fclose(fid);
189	baylor	1.1
190	baylor	1.2	Rvec=1031:.05:1032;
191	baylor	1.1
192			subplot(4,2,3)
193			dens=thetarhotAlpha+rho;
194	baylor	1.2	[h,c]=contour(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(dens(:,1,:)),Rvec);
195			title('Density')
196			xlabel('x (km)');ylabel('z (m)')
197
198			subplot(4,2,4)
199			[h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)),Rvec);
200	baylor	1.1	title('Density')
201	baylor	1.2	colorbar
202	baylor	1.1	xlabel('y (km)');ylabel('z (m)')
203
204			%Spice
205
206			x1=xc(round(length(xc)/4)); % centered location of Front #1
207			x2=xc(round(3*length(xc)/4)); % centered location of Front #2
208	baylor	1.2
209			%spice=((SL+SR)+(SR-SL)*tanh((ZT-Dml/2)/Dml))/2
210
211			spice=T0+dthetadzint(ZT+Dml)+deltatheta(tanh((XT-x1)/Lf)-tanh((XT-x2)/Lf)-1)/2;
212	baylor	1.1
213			i=1
214	baylor	1.2
215	baylor	1.1	while (i<iml)
216	baylor	1.2	spice(:,:,i)=T0+dthetadzml(ZT(:,:,i)+Dml)+deltatheta(tanh((XT(:,:,i)-x1)/Lf)-tanh((XT(:,:,i)-x2)/Lf)-1)/2;
217	baylor	1.1	i=i+1;
218			end
219
220	baylor	1.2	fid=fopen('spiceRestoreFile.bin','w','b'); fwrite(fid,spice,'real*8'); fclose(fid);
221	baylor	1.1
222			subplot(4,2,8)
223	baylor	1.2	[h,c]=contourf(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(spice(:,1,:)),Tvec);
224			%clabel(h,c)
225	baylor	1.1	title('Spice')
226			xlabel('x (km)');ylabel('z (m)')
227
228	baylor	1.2	%clabel(h,c);
229	baylor	1.1
230			subplot(4,2,7)
231	baylor	1.2	invtime=0*YT;
232	baylor	1.1	decay=dyspacing;
233	baylor	1.2	invtime=1/time_scale*(exp(-YT/decay+YT(1,1,1)/decay)+exp(YT/decay-YT(1,end-1,1)/decay));
234			% Surface Restoring
235			invtime(:,:,1)=invtime(:,:,1)+1/surf_time_scale;
236
237			figure(4)
238			subplot(2,1,1)
239			plot(squeeze(YT(1,1:end-1,:)/1e3),squeeze(invtime(1,1:end-1,:)*86400))
240			ylabel('inverse restoring time (1/day)')
241			xlabel('Y (km)')
242
243			subplot(2,1,2)
244			plot(squeeze(invtime(1,1:end-1,:)*86400)',squeeze(ZT(1,1:end-1,:))')
245			xlabel('inverse restoring time (1/day)')
246			ylabel('depth')
247	baylor	1.1
248			figure(2)
249			plot(YB(1,:)/1e3,hh(1,:))
250			title('topography')
251			xlabel('x (km)');ylabel('value')
252
253	baylor	1.2	fid=fopen('invtimeRestoreFile.bin','w','b'); fwrite(fid,invtime,'real*8'); fclose(fid);
254
255	baylor	1.1	%-- Perturb Initial temperature
256	baylor	1.2
257	baylor	1.1	pert=rand(size(theta(:,:,1)));
258	baylor	1.2
259			pert=1e-4*(pert-0.5);
260
261	baylor	1.1	for i=1:length(theta(1,1,:))
262			theta(:,:,i)=theta(:,:,i)+pert;
263			end
264
265	baylor	1.2	fid=fopen('thetaInitial.bin','w','b');
266			fwrite(fid,theta,'real*8'); fclose(fid);
267
268			fid=fopen('spiceInitial.bin','w','b');
269			fwrite(fid,spice,'real*8'); fclose(fid);
270
271			%-- Thermal wind
272			U = 0*XU;
273
274			for k=(nz-1):-1:1
275			U(1:end-1,1:end-1,k) = U(1:end-1,1:end-1,k+1)+dh(k)g/f0/rhodiff(dens(:,:,k),1,2)./diff(YU(1:end-1,:,k+1),1,2);
276			U(:,end,k) = 0; % Get endpoint
277			U(end,:,k) = U(1,:,k); % In the wall
278			end
279
280			% Make initial state Momentumless (at least boussinesq, good approx)
281
282			[imax,jmax,kmax]=size(U);
283
284			for i=1:imax
285			for j=1:jmax
286			U(i,j,:)=U(i,j,:)-mean(U(i,j,:));
287			end
288			end
289
290			timestep=dxspacing/max(abs(U(:)))/3.
291
292			subplot(4,2,5)
293			[h,c]=contour(squeeze(XU(:,1,:))/1e3,squeeze(ZU(:,1,:)),squeeze(U(:,1,:)),-1:.02:1);
294			title('U');xlabel('x (km)');ylabel('z (m)')
295
296			subplot(4,2,6)
297			[h,c]=contour(squeeze(YU(1,:,:))/1e3,squeeze(ZU(1,:,:)),squeeze(U(1,:,:)),-1:.02:1);
298
299			title('U');xlabel('x (km)');ylabel('z (m)')
300
301			U2=(U([1:end-1],:,:)+U([2:end],:,:))/2;
302
303			fid=fopen('uInitial.bin','w','b');
304			fwrite(fid,U2,'real*8'); fclose(fid);
305
306
307			%-- Wind stress
308
309			%tau=tau0*tanh((((squeeze(YU(:,:,1))-y0-Lf))/Lf));
310
311			% Constant Wind Stress
312			%tau=tau0*ones(size(squeeze(YU(:,:,1))));
313	baylor	1.1
314	baylor	1.2	%subplot(4,2,8)
315			%quiver(squeeze(XU(:,:,1))/1e3,squeeze(YU(:,:,1))/1e3,tau,0*tau)
316			%title('Taux');xlabel('x (km)');ylabel('y (km)')
317			%fid=fopen('taux_sl.bin','w','b');
318			%fwrite(fid,tau,'real*8'); fclose(fid);
319	baylor	1.1
320	baylor	1.2	%save 'zf' zf
321	baylor	1.1