| 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 |
dxspacing=1e3 |
| 12 |
dyspacing=dxspacing |
| 13 |
|
| 14 |
Lx=dxspacing*nx |
| 15 |
Ly=dyspacing*ny |
| 16 |
|
| 17 |
%-- Params |
| 18 |
g=9.81; |
| 19 |
tAlpha=-2e-4; |
| 20 |
f0=7.29e-5; |
| 21 |
rho=1035 |
| 22 |
day=24*60^2; |
| 23 |
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 |
|
| 29 |
H=2000; %Max depth |
| 30 |
|
| 31 |
%-- Grid: x |
| 32 |
dx_ratio=30; dx_trans=.05; dx_min=100; |
| 33 |
xn=(0.5:nx)/nx; |
| 34 |
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 |
dy_ratio=30; dy_trans=.02; dy_min=100; |
| 41 |
yn=(0.5:ny)/ny; |
| 42 |
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 |
|
| 50 |
dh=H/nz*ones(1,nz); |
| 51 |
%for i=21:38 |
| 52 |
% dh(i)=dh(i-1)*1.1917; |
| 53 |
%end |
| 54 |
|
| 55 |
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 |
Lf=64*1e3; % Half-Width of Front |
| 71 |
%Lf=50*dyspacing; % Half-Width of Front |
| 72 |
y0=yc(round(length(yc)/2)); % centered location of Front |
| 73 |
|
| 74 |
|
| 75 |
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/f0*Dml, Ms/f0^2*Dml) |
| 81 |
Ls=(f0^2/Ms/Dml*sqrt(5/2/(1+max(Ri,1))))^-1*2*pi |
| 82 |
|
| 83 |
Riint=Nsint*f0^2/Ms^2; |
| 84 |
Ldint=max(Nsint^.5/f0*H, Ms/f0^2*H) |
| 85 |
Lsint=(f0^2/Ms/H*sqrt(5/2/(1+max(Riint,1))))^-1*2*pi |
| 86 |
|
| 87 |
Resfacs=[Lf/Ld Lf/Ls Ls/dxspacing Ld/dxspacing] |
| 88 |
|
| 89 |
dxspacing |
| 90 |
optdxspacing=round(Ls*0.05) |
| 91 |
|
| 92 |
deltatheta=2*Lf*Ms/g/tAlpha |
| 93 |
T0=17 |
| 94 |
dthetadzml=-Nsml/g/tAlpha |
| 95 |
dthetadzint=-Nsint/g/tAlpha |
| 96 |
|
| 97 |
%-- Bathymetry |
| 98 |
% Bathymetry is on Xc, Yc |
| 99 |
hh=ones(size(XB)); |
| 100 |
|
| 101 |
% Variable Bathymetry |
| 102 |
% hh(:,:)=hh(:,:)-2*max(0,abs(YB/L-1/2)-3/8); |
| 103 |
|
| 104 |
hh(:,end)=0*hh(:,end); |
| 105 |
|
| 106 |
|
| 107 |
hh=-H*hh; |
| 108 |
|
| 109 |
Ltop=L/2; |
| 110 |
|
| 111 |
figure(1) |
| 112 |
subplot(221) |
| 113 |
pcolor(XB/1e3,YB/1e3,hh) |
| 114 |
|
| 115 |
|
| 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 |
% 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(2*pi*t/day),cos(pi*daylength/day))-cos(pi*daylength/day))/length(t); |
| 135 |
|
| 136 |
qdiur=-qconst/q1; |
| 137 |
|
| 138 |
q1=qconst+qdiur*(max(cos(2*pi*t/day),cos(pi*daylength/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 |
|
| 153 |
% Initial Temp Profile |
| 154 |
|
| 155 |
figure(2) |
| 156 |
subplot(4,2,1) |
| 157 |
|
| 158 |
theta=T0+dthetadzint*(ZT+Dml)+deltatheta*tanh((YT-y0)/Lf)/2; |
| 159 |
|
| 160 |
% Impose a strong Dml |
| 161 |
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 |
theta(:,:,i)=T0+dthetadzml*(ZT(:,:,i)+Dml)+deltatheta*tanh((YT(:,:,i)-y0)/Lf)/2; |
| 170 |
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 |
Tvec=15:.1:18; |
| 176 |
|
| 177 |
subplot(4,2,1) |
| 178 |
[h,c]=contourf(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(theta(:,1,:)),Tvec); |
| 179 |
title('Potl Temp') |
| 180 |
xlabel('x (km)');ylabel('z (m)') |
| 181 |
fid=fopen('thetaRestoreFile.bin','w','b'); fwrite(fid,theta,'real*8'); fclose(fid); |
| 182 |
|
| 183 |
subplot(4,2,2) |
| 184 |
[h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)),Tvec); |
| 185 |
colorbar |
| 186 |
title('Potl Temp') |
| 187 |
xlabel('y (km)');ylabel('z (m)') |
| 188 |
fid=fopen('thetaRestoreFile.bin','w','b'); fwrite(fid,theta,'real*8'); fclose(fid); |
| 189 |
|
| 190 |
Rvec=1031:.05:1032; |
| 191 |
|
| 192 |
subplot(4,2,3) |
| 193 |
dens=theta*rho*tAlpha+rho; |
| 194 |
[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 |
title('Density') |
| 201 |
colorbar |
| 202 |
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 |
|
| 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 |
|
| 213 |
i=1 |
| 214 |
|
| 215 |
while (i<iml) |
| 216 |
spice(:,:,i)=T0+dthetadzml*(ZT(:,:,i)+Dml)+deltatheta*(tanh((XT(:,:,i)-x1)/Lf)-tanh((XT(:,:,i)-x2)/Lf)-1)/2; |
| 217 |
i=i+1; |
| 218 |
end |
| 219 |
|
| 220 |
fid=fopen('spiceRestoreFile.bin','w','b'); fwrite(fid,spice,'real*8'); fclose(fid); |
| 221 |
|
| 222 |
subplot(4,2,8) |
| 223 |
[h,c]=contourf(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(spice(:,1,:)),Tvec); |
| 224 |
%clabel(h,c) |
| 225 |
title('Spice') |
| 226 |
xlabel('x (km)');ylabel('z (m)') |
| 227 |
|
| 228 |
%clabel(h,c); |
| 229 |
|
| 230 |
subplot(4,2,7) |
| 231 |
invtime=0*YT; |
| 232 |
decay=dyspacing; |
| 233 |
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 |
|
| 248 |
figure(2) |
| 249 |
plot(YB(1,:)/1e3,hh(1,:)) |
| 250 |
title('topography') |
| 251 |
xlabel('x (km)');ylabel('value') |
| 252 |
|
| 253 |
fid=fopen('invtimeRestoreFile.bin','w','b'); fwrite(fid,invtime,'real*8'); fclose(fid); |
| 254 |
|
| 255 |
%-- Perturb Initial temperature |
| 256 |
|
| 257 |
pert=rand(size(theta(:,:,1))); |
| 258 |
|
| 259 |
pert=1e-4*(pert-0.5); |
| 260 |
|
| 261 |
for i=1:length(theta(1,1,:)) |
| 262 |
theta(:,:,i)=theta(:,:,i)+pert; |
| 263 |
end |
| 264 |
|
| 265 |
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/rho*diff(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 |
|
| 314 |
%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 |
|
| 320 |
%save 'zf' zf |
| 321 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph