/[MITgcm]/MITgcm/verification/MLAdjust/input.0.leithD/genit.m
ViewVC logotype

Contents of /MITgcm/verification/MLAdjust/input.0.leithD/genit.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (show annotations) (download)
Thu Feb 13 00:06:19 2014 UTC (10 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +0 -0 lines
FILE REMOVED
dir renamed to "input.AhFlxF"

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=1000
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 prec='real*8';
24 ieee='b';
25
26 H=200; %Max depth
27
28 %-- Grid: x
29 dx=ones(1,nx); % uniform resolution
30 dx=dx*Lx/sum(dx);
31 xf=cumsum([0 dx]); % Face x points
32 xc=(xf(1:end-1)+xf(2:end))/2; % Centered x points
33
34 %-- Grid: y
35 dy=ones(1,ny); % uniform resolution
36 dy=dy*Ly/sum(dy);
37 yf=cumsum([0 dy]); % Face y-points
38 yc=(yf(1:end-1)+yf(2:end))/2; %Centered y-points
39 L=yc(end)-yc(1); % this takes into account the wall of topography!!!
40
41 %-- Grid: z
42 dh=H/nz*ones(1,nz);
43 zf=-round(cumsum([0 dh])/1)*1; % Face z points
44 dh=-diff(zf);
45 zc=(zf(1:end-1)+zf(2:end))/2; % centered z points
46 nz=length(dh);
47 H=sum(dh)
48
49 [XT,YT,ZT]=ndgrid(xc,yc,zc); % This is the centered, temperature grid.
50 [XU,YU,ZU]=ndgrid(xf,yc,zc); % This is the U grid.
51 [XV,YV,ZV]=ndgrid(xc,yf,zc); % This is the V grid.
52 [XW,YW,ZW]=ndgrid(xc,yc,zf); % This is the W grid.
53 [XB,YB]=ndgrid(xc,yc); % This is the Bathymetry grid.
54
55 % Stratification Params
56 Dml=200; % Mixed-layer Depth
57 D=Dml;
58 Lf=2*1e3; % Half-Width of Front
59 y0=yc(round(length(yc)/2)); % centered location of Front
60
61 Ms=(2e-4)^2 %db/dy 2e-4 is about 1 K/50 km
62 Ns=9*Ms^2/f0^2 %db/dz
63
64 deltatheta=Lf*Ms/g/tAlpha
65 T0=17
66 dthetadz=-Ns/g/tAlpha
67
68 %-- Bathymetry
69 % Bathymetry is on Xc, Yc
70 hh=ones(size(XB));
71 hh(:,end)=0*hh(:,end);
72 hh=-H*hh;
73
74 figure(1)
75 subplot(221)
76 pcolor(XB/1e3,YB/1e3,hh)
77 axis equal
78
79 subplot(223)
80 plot(xc/1e3,dx/1e3,'.'); xlabel('X (km)');ylabel('\Delta x (km)')
81 title('\Delta x')
82
83 subplot(222)
84 plot(dy/1e3,yc/1e3,'.'); ylabel('Y (km)');xlabel('\Delta y (km)')
85 title('\Delta y')
86
87 subplot(224)
88 plot(dh,zc,'.'); ylabel('Z (m)');xlabel('\Delta z (m)')
89 title('\Delta z')
90
91 fid=fopen('topo_sl.bin','w',ieee); fwrite(fid,hh,prec); fclose(fid);
92
93 % Initial Temp Profile
94 figure(2)
95 subplot(4,2,1)
96 theta=T0+dthetadz*(ZT-ZT(1,1,end))+deltatheta*tanh((YT-y0)/Lf)/2;
97
98 % Impose a strong initial and restoring mixed layer to Dml
99 i=1;
100 iml=1;
101
102 while ((ZT(1,1,iml)>-Dml)&(iml<length(ZT(1,1,:))))
103 iml=iml+1;
104 end
105
106 while (i<iml)
107 theta(:,:,i)=theta(:,:,iml);
108 i=i+1;
109 end
110
111 [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)]
112
113
114 subplot(4,2,1)
115 [h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)));
116 axis([(yc(round(length(yc)/2))-1.5*Lf)/1e3 (yc(round(length(yc)/2))+1.5*Lf)/1e3 -Dml 0])
117 colorbar
118 title('Potl Temp')
119 xlabel('y (km)');ylabel('z (m)')
120
121 subplot(4,2,2)
122 [h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)));
123 colorbar
124 title('Potl Temp')
125 xlabel('y (km)');ylabel('z (m)')
126
127
128 subplot(4,2,3)
129 dens=theta*rho*tAlpha+rho;
130 [h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)));
131 title('Density')
132 xlabel('y (km)');ylabel('z (m)')
133 axis([(yc(round(length(yc)/2))-1.5*Lf)/1e3 (yc(round(length(yc)/2))+1.5*Lf)/1e3 -Dml 0])
134
135 %Spice
136
137 x1=xc(round(length(xc)/4)); % centered location of Front #1
138 x2=xc(round(3*length(xc)/4)); % centered location of Front #2
139 spice=T0+dthetadz*(ZT-ZT(1,1,end))+deltatheta*(tanh((XT-x1)/Lf)-tanh((XT-x2)/Lf)-1)/2;
140
141 i=1
142 while (i<iml)
143 spice(:,:,i)=spice(:,:,iml);
144 i=i+1;
145 end
146
147
148 subplot(4,2,8)
149 [h,c]=contourf(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(spice(:,1,:)));
150 title('Spice')
151 xlabel('x (km)');ylabel('z (m)')
152
153 subplot(4,2,4)
154 [h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)));
155 title('Density')
156 xlabel('y (km)');ylabel('z (m)')
157
158 subplot(4,2,7)
159 decay=dyspacing;
160
161 figure(2)
162 plot(YB(1,:)/1e3,hh(1,:))
163 title('topography')
164 xlabel('x (km)');ylabel('value')
165
166 %-- Perturb Initial temperature
167 pert=rand(size(theta(:,:,1)));
168 pert=1e-5*(pert-0.5);
169 for i=1:length(theta(1,1,:))
170 theta(:,:,i)=theta(:,:,i)+pert;
171 end
172
173 fid=fopen('thetaInitial.bin','w',ieee);
174 fwrite(fid,theta,prec); fclose(fid);
175
176 fid=fopen('spiceInitial.bin','w',ieee);
177 fwrite(fid,spice,prec); fclose(fid);
178
179

  ViewVC Help
Powered by ViewVC 1.1.22