offline_exf_seaice/input/gendata.m


kwr=1; kprt=0;
kwr=0; kprt=0;
nx=80; ny=42; nr=3; nt=1;

xc=[1:nx]; xc=xc-mean(xc);
yc=[1:ny]-.5;

%------------------------------------------------------

windx=10.;
H0=-100.;

namf='channel.bin';
depth=H0*ones(nx,ny); depth(:,1)=0.;
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,depth,'real*8'); fclose(fid);
end

namf='bathy_3c.bin';
msk=abs(xc)'*ones(1,ny)+ones(nx,1)*yc;
depth=H0*ones(nx,ny); depth(:,1)=0.;
depth(find(msk < 24))=0.;
y2d=ones(nx,1)*yc;
depth(find(y2d > ny/2))=H0;
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,depth,'real*8'); fclose(fid);
end

namf='windx.bin';
wnd=windx*ones(nx,ny,nt);
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,wnd,'real*8'); fclose(fid);
end

%- file name convention: "const_{xx}.bin" <-> uniform value = xx (in percent)
namf='const_00.bin';
fld=0*ones(nx,ny,nt);
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
end

namf='const100.bin'; w0=1.;
var=w0*ones(nx,ny);
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
end

namf='const+20.bin'; w0=0.2;
var=w0*ones(nx,ny);
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
end

%namf='const+40.bin';
%u0=0.4;
%var=u0*ones(nx,ny);
%if kwr > 0,
% fprintf('write to file: %s\n',namf);
% fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
%end

%namf='const-10.bin';
%v0=-0.1;
%var=v0*ones(nx,ny);
%if kwr > 0,
% fprintf('write to file: %s\n',namf);
% fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
%end
%------------------------------------------------------

namf='ice0_area.bin'; iceC0=1.;
iceConc=iceC0*ones(nx,ny); iceConc(:,1)=0;
iceConc(:,2)=0.00*iceC0;
iceConc(:,3)=0.10*iceC0;
iceConc(:,end)  =0.00*iceC0;
iceConc(:,end-1)=0.01*iceC0;
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,iceConc,'real*8'); fclose(fid);
end

namf='ice0_heff.bin'; iceH0=0.2;
iceVol=iceConc*iceH0;
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,iceVol,'real*8'); fclose(fid);
end

%------------------------------------------------------

dsw0=100;
namf=['dsw_',int2str(dsw0),'.bin'];
fld=dsw0*ones(nx,ny,nt);
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
end

dlw0=250;
namf=['dlw_',int2str(dlw0),'.bin'];
fld=dlw0*ones(nx,ny,nt);
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
end

cel2K=273.15; dtx=4; %- dtx = amplitude of air temp variations in X-dir
ta_x=cel2K + dtx*sin(pi*(1+2*xc'/nx));
ta=repmat(ta_x,[1 ny nt]);
namf=['tair_',int2str(dtx),'x.bin'];
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,ta,'real*8'); fclose(fid);
end;

cvapor_fac     =   640380.000  ;
cvapor_exp     =     5107.400  ;
atmrho         =        1.200  ;
rh=70; %- specific humid <--> 70.% relative humid
tmpbulk = cvapor_fac*exp(-cvapor_exp./ta_x);
qa_x = (rh/100.)*tmpbulk/atmrho ;
qa=repmat(qa_x,[1 ny nt]);
namf=['qa',int2str(rh),'_',int2str(dtx),'x.bin'];
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,qa,'real*8'); fclose(fid);
end;

%- salinity
sCst=30;
so=sCst*ones(nx,ny,nt);
namf='socn.bin';
%if kwr > 0,
% fprintf('write to file: %s\n',namf);
% fid=fopen(namf,'w','b'); fwrite(fid,so,'real*8'); fclose(fid);
%end;

muTf = 5.4e-2;
tfreeze=-muTf*sCst;
fprintf('T-freeze = %10.6f\n',tfreeze);
%- parabolic profile in Y, max @ j=4, min @ j=ny, amplitude=1.K
to_y=(yc-3.5)/(ny-4);
to_y=tfreeze+0.5-to_y.*to_y;
  mnV=min(to_y); MxV=max(to_y); Avr=mean(to_y(2:end));
  fprintf(' SST* av,mn,Mx: %9.6f , %9.6f , %9.6f , %9.6f\n',Avr,mnV,MxV,MxV-mnV);
to=repmat(to_y,[nx 1 nt]);
namf='tocn.bin';
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,to,'real*8'); fclose(fid);
end;

%-- make some plots to check: ----------------

hScal=[-1.1 0.1]*abs(H0);
figure(1); clf;
subplot(211);
var=depth;
imagesc(xc,yc,var'); set(gca,'YDir','normal');
%caxis(hScal);
%change_colmap(-1);
colorbar;
title('Depth [m]');

subplot(413);
var=depth;
j1=2;
j2=ny/2;
j3=j2+1;
plot(xc,var(:,j1),'k-')
hold on; j=j+1;
plot(xc,var(:,j2),'ro-')
plot(xc,var(:,j3),'b-')
hold off;
axis([-nx/2 nx/2 hScal]);
grid
legend(int2str(j1),int2str(j2),int2str(j3));
title('Depth @ j= cst');

subplot(414);
i=nx/2;
plot(yc,var(i,:),'k-')
axis([0 ny H0*1.1 -H0*.1]);
grid
title(['Depth @ i=',int2str(i)]);

%--
dewPt=(qa_x*atmrho)/cvapor_fac;
dewPt=-cvapor_exp./log(dewPt);

figure(2);clf;
subplot(311)
P(1)=plot(xc,ta_x-cel2K,'r-'); hold on;
P(2)=plot(xc,dewPt-cel2K,'b-');
P(3)=plot(xc,tfreeze*ones(nx,1),'k-');
set(P,'LineWidth',1);
hold off; AA=axis;
axis([-nx/2 nx/2 AA(3:4)]);
legend('ta','dew');
grid
xlabel('X')
title(['Air Temp (^oC): del-Temp-X = ',int2str(dtx),' , RH= ',int2str(rh)]);
subplot(312)
P(1)=plot(yc,to_y,'b-'); hold on;
P(2)=plot(yc,tfreeze*ones(ny,1),'k-');
set(P,'LineWidth',1);
hold off; AA=axis;
L=line([1 1],AA(3:4)); set(L,'LineWidth',2.,'Color',[0 0 0]);
axis([0 ny AA(3:4)]);
grid
xlabel('Y')
title('Ocean Temp ^oC');

subplot(313)
var=iceConc(1,:);
P(1)=semilogy(yc,var,'b-x'); hold on;
%plot(yc,var,'b-x'); hold on;
var=iceVol(1,:);
P(2)=semilogy(yc,var,'r-x');
%plot(yc,var,'r-x');
set(P,'LineWidth',1);
hold off; AA=axis;
L=line([1 1],AA(3:4)); set(L,'LineWidth',2.,'Color',[0 0 0]);
axis([0 ny [0 2]*iceC0]);
grid
xlabel('Y')
legend('iceC','hEff','Location','South');
title('Initial ice in Channel : y-section');
%-----
if kprt == 1, f=2;
 namfig=sprintf('forcing_%2.2i',f);
 fprintf([' print fig= %2i to file: ',namfig,' '],f);
 set(f,'PaperOrientation','portrait')
%set(f,'PaperPosition',[0.25 1.5 6. 8.]);
 set(f,'PaperPosition',[0.25 1.5 5.25 7.]);
 print(f,'-depsc2',namfig); fprintf('\n');
end
%-----

figure(3);clf;
subplot(311)
var=iceConc; ccB=[-1 12]/10;
imagesc(xc,yc,var'); set(gca,'YDir','normal');
caxis(ccB);
%change_colmap(-1);
colorbar;
title('Ice Concentration in Channel');

subplot(312)
var=iceVol; ccB=[-1 12]/50;
imagesc(xc,yc,var'); set(gca,'YDir','normal');
caxis(ccB);
%change_colmap(-1);
colorbar;
title('Effective ice thickness in Channel');

subplot(313)
var=iceConc(1,:);
%plot(yc,var,'b-x'); hold on;
semilogy(yc,var,'b-x'); hold on;
var=iceVol(1,:);
%plot(yc,var,'r-x'); hold off;
semilogy(yc,var,'r-x'); hold on;
AA=axis; axis([0 ny [0 2]*iceC0]);
grid
legend('iceC','hEff','Location','South');
title('Initial ice in Channel : y-section');
%-----

return
1	jmc	1.1
2	jmc	1.4	kwr=1; kprt=0;
3			kwr=0; kprt=0;
4	jmc	1.3	nx=80; ny=42; nr=3; nt=1;
5	jmc	1.1
6			xc=[1:nx]; xc=xc-mean(xc);
7			yc=[1:ny]-.5;
8
9			%------------------------------------------------------
10
11			windx=10.;
12			H0=-100.;
13
14			namf='channel.bin';
15			depth=H0*ones(nx,ny); depth(:,1)=0.;
16	gforget	1.2	if kwr > 0,
17			fprintf('write to file: %s\n',namf);
18			fid=fopen(namf,'w','b'); fwrite(fid,depth,'real*8'); fclose(fid);
19			end
20	jmc	1.1
21			namf='bathy_3c.bin';
22			msk=abs(xc)'ones(1,ny)+ones(nx,1)yc;
23			depth=H0*ones(nx,ny); depth(:,1)=0.;
24			depth(find(msk < 24))=0.;
25			y2d=ones(nx,1)*yc;
26			depth(find(y2d > ny/2))=H0;
27			if kwr > 0,
28			fprintf('write to file: %s\n',namf);
29			fid=fopen(namf,'w','b'); fwrite(fid,depth,'real*8'); fclose(fid);
30			end
31
32			namf='windx.bin';
33			wnd=windx*ones(nx,ny,nt);
34			if kwr > 0,
35			fprintf('write to file: %s\n',namf);
36			fid=fopen(namf,'w','b'); fwrite(fid,wnd,'real*8'); fclose(fid);
37			end
38
39			%- file name convention: "const_{xx}.bin" <-> uniform value = xx (in percent)
40			namf='const_00.bin';
41			fld=0*ones(nx,ny,nt);
42			if kwr > 0,
43			fprintf('write to file: %s\n',namf);
44			fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
45			end
46
47			namf='const100.bin'; w0=1.;
48			var=w0*ones(nx,ny);
49			if kwr > 0,
50			fprintf('write to file: %s\n',namf);
51			fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
52			end
53
54			namf='const+20.bin'; w0=0.2;
55			var=w0*ones(nx,ny);
56			if kwr > 0,
57			fprintf('write to file: %s\n',namf);
58			fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
59			end
60
61			%namf='const+40.bin';
62			%u0=0.4;
63			%var=u0*ones(nx,ny);
64			%if kwr > 0,
65			% fprintf('write to file: %s\n',namf);
66			% fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
67			%end
68
69			%namf='const-10.bin';
70			%v0=-0.1;
71			%var=v0*ones(nx,ny);
72			%if kwr > 0,
73			% fprintf('write to file: %s\n',namf);
74			% fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
75			%end
76			%------------------------------------------------------
77
78	jmc	1.3	namf='ice0_area.bin'; iceC0=1.;
79			iceConc=iceC0*ones(nx,ny); iceConc(:,1)=0;
80			iceConc(:,2)=0.00*iceC0;
81			iceConc(:,3)=0.10*iceC0;
82			iceConc(:,end) =0.00*iceC0;
83			iceConc(:,end-1)=0.01*iceC0;
84			if kwr > 0,
85	gforget	1.2	fprintf('write to file: %s\n',namf);
86	jmc	1.3	fid=fopen(namf,'w','b'); fwrite(fid,iceConc,'real*8'); fclose(fid);
87	gforget	1.2	end
88
89	jmc	1.3	namf='ice0_heff.bin'; iceH0=0.2;
90			iceVol=iceConc*iceH0;
91			if kwr > 0,
92	gforget	1.2	fprintf('write to file: %s\n',namf);
93	jmc	1.3	fid=fopen(namf,'w','b'); fwrite(fid,iceVol,'real*8'); fclose(fid);
94	gforget	1.2	end
95
96			%------------------------------------------------------
97
98	jmc	1.1	dsw0=100;
99			namf=['dsw_',int2str(dsw0),'.bin'];
100			fld=dsw0*ones(nx,ny,nt);
101			if kwr > 0,
102			fprintf('write to file: %s\n',namf);
103			fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
104			end
105
106	jmc	1.3	dlw0=250;
107	jmc	1.1	namf=['dlw_',int2str(dlw0),'.bin'];
108			fld=dlw0*ones(nx,ny,nt);
109			if kwr > 0,
110			fprintf('write to file: %s\n',namf);
111			fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
112			end
113
114			cel2K=273.15; dtx=4; %- dtx = amplitude of air temp variations in X-dir
115			ta_x=cel2K + dtxsin(pi(1+2*xc'/nx));
116			ta=repmat(ta_x,[1 ny nt]);
117			namf=['tair_',int2str(dtx),'x.bin'];
118			if kwr > 0,
119			fprintf('write to file: %s\n',namf);
120			fid=fopen(namf,'w','b'); fwrite(fid,ta,'real*8'); fclose(fid);
121			end;
122
123			cvapor_fac = 640380.000 ;
124			cvapor_exp = 5107.400 ;
125			atmrho = 1.200 ;
126			rh=70; %- specific humid <--> 70.% relative humid
127			tmpbulk = cvapor_fac*exp(-cvapor_exp./ta_x);
128			qa_x = (rh/100.)*tmpbulk/atmrho ;
129			qa=repmat(qa_x,[1 ny nt]);
130			namf=['qa',int2str(rh),'_',int2str(dtx),'x.bin'];
131			if kwr > 0,
132			fprintf('write to file: %s\n',namf);
133			fid=fopen(namf,'w','b'); fwrite(fid,qa,'real*8'); fclose(fid);
134			end;
135
136			%- salinity
137			sCst=30;
138			so=sCst*ones(nx,ny,nt);
139			namf='socn.bin';
140			%if kwr > 0,
141			% fprintf('write to file: %s\n',namf);
142			% fid=fopen(namf,'w','b'); fwrite(fid,so,'real*8'); fclose(fid);
143			%end;
144
145			muTf = 5.4e-2;
146	jmc	1.3	tfreeze=-muTf*sCst;
147	jmc	1.1	fprintf('T-freeze = %10.6f\n',tfreeze);
148	jmc	1.3	%- parabolic profile in Y, max @ j=4, min @ j=ny, amplitude=1.K
149			to_y=(yc-3.5)/(ny-4);
150	jmc	1.1	to_y=tfreeze+0.5-to_y.*to_y;
151	jmc	1.3	mnV=min(to_y); MxV=max(to_y); Avr=mean(to_y(2:end));
152			fprintf(' SST* av,mn,Mx: %9.6f , %9.6f , %9.6f , %9.6f\n',Avr,mnV,MxV,MxV-mnV);
153	jmc	1.1	to=repmat(to_y,[nx 1 nt]);
154			namf='tocn.bin';
155			if kwr > 0,
156			fprintf('write to file: %s\n',namf);
157			fid=fopen(namf,'w','b'); fwrite(fid,to,'real*8'); fclose(fid);
158			end;
159
160			%-- make some plots to check: ----------------
161
162			hScal=[-1.1 0.1]*abs(H0);
163			figure(1); clf;
164			subplot(211);
165			var=depth;
166			imagesc(xc,yc,var'); set(gca,'YDir','normal');
167			%caxis(hScal);
168			%change_colmap(-1);
169			colorbar;
170			title('Depth [m]');
171
172			subplot(413);
173			var=depth;
174			j1=2;
175			j2=ny/2;
176			j3=j2+1;
177			plot(xc,var(:,j1),'k-')
178			hold on; j=j+1;
179			plot(xc,var(:,j2),'ro-')
180			plot(xc,var(:,j3),'b-')
181			hold off;
182			axis([-nx/2 nx/2 hScal]);
183			grid
184			legend(int2str(j1),int2str(j2),int2str(j3));
185			title('Depth @ j= cst');
186
187			subplot(414);
188			i=nx/2;
189			plot(yc,var(i,:),'k-')
190			axis([0 ny H01.1 -H0.1]);
191			grid
192			title(['Depth @ i=',int2str(i)]);
193
194			%--
195			dewPt=(qa_x*atmrho)/cvapor_fac;
196			dewPt=-cvapor_exp./log(dewPt);
197
198			figure(2);clf;
199	jmc	1.4	subplot(311)
200			P(1)=plot(xc,ta_x-cel2K,'r-'); hold on;
201			P(2)=plot(xc,dewPt-cel2K,'b-');
202			P(3)=plot(xc,tfreeze*ones(nx,1),'k-');
203			set(P,'LineWidth',1);
204			hold off; AA=axis;
205			axis([-nx/2 nx/2 AA(3:4)]);
206	jmc	1.1	legend('ta','dew');
207			grid
208	jmc	1.4	xlabel('X')
209			title(['Air Temp (^oC): del-Temp-X = ',int2str(dtx),' , RH= ',int2str(rh)]);
210			subplot(312)
211			P(1)=plot(yc,to_y,'b-'); hold on;
212			P(2)=plot(yc,tfreeze*ones(ny,1),'k-');
213			set(P,'LineWidth',1);
214			hold off; AA=axis;
215			L=line([1 1],AA(3:4)); set(L,'LineWidth',2.,'Color',[0 0 0]);
216			axis([0 ny AA(3:4)]);
217	jmc	1.1	grid
218	jmc	1.4	xlabel('Y')
219	jmc	1.1	title('Ocean Temp ^oC');
220	jmc	1.3
221	jmc	1.4	subplot(313)
222			var=iceConc(1,:);
223			P(1)=semilogy(yc,var,'b-x'); hold on;
224			%plot(yc,var,'b-x'); hold on;
225			var=iceVol(1,:);
226			P(2)=semilogy(yc,var,'r-x');
227			%plot(yc,var,'r-x');
228			set(P,'LineWidth',1);
229			hold off; AA=axis;
230			L=line([1 1],AA(3:4)); set(L,'LineWidth',2.,'Color',[0 0 0]);
231			axis([0 ny [0 2]*iceC0]);
232			grid
233			xlabel('Y')
234			legend('iceC','hEff','Location','South');
235			title('Initial ice in Channel : y-section');
236			%-----
237			if kprt == 1, f=2;
238			namfig=sprintf('forcing_%2.2i',f);
239			fprintf([' print fig= %2i to file: ',namfig,' '],f);
240			set(f,'PaperOrientation','portrait')
241			%set(f,'PaperPosition',[0.25 1.5 6. 8.]);
242			set(f,'PaperPosition',[0.25 1.5 5.25 7.]);
243			print(f,'-depsc2',namfig); fprintf('\n');
244			end
245			%-----
246	jmc	1.3
247			figure(3);clf;
248			subplot(311)
249			var=iceConc; ccB=[-1 12]/10;
250			imagesc(xc,yc,var'); set(gca,'YDir','normal');
251			caxis(ccB);
252			%change_colmap(-1);
253			colorbar;
254			title('Ice Concentration in Channel');
255
256			subplot(312)
257			var=iceVol; ccB=[-1 12]/50;
258			imagesc(xc,yc,var'); set(gca,'YDir','normal');
259			caxis(ccB);
260			%change_colmap(-1);
261			colorbar;
262			title('Effective ice thickness in Channel');
263
264			subplot(313)
265			var=iceConc(1,:);
266			%plot(yc,var,'b-x'); hold on;
267			semilogy(yc,var,'b-x'); hold on;
268			var=iceVol(1,:);
269			%plot(yc,var,'r-x'); hold off;
270			semilogy(yc,var,'r-x'); hold on;
271			AA=axis; axis([0 ny [0 2]*iceC0]);
272			grid
273			legend('iceC','hEff','Location','South');
274			title('Initial ice in Channel : y-section');
275	jmc	1.4	%-----
276
277			return