offline_cheapaml/input/gendata.m


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

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

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

windx=10.;
H0=-100.;

namf='channel.bin';
depth=H0*ones(nx,ny); depth(:,1)=0.; depth(:,ny)=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=['windx_',int2str(windx),'ms.bin'];
uwind=windx*ones(nx,ny,nt);
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,uwind,'real*8'); fclose(fid);
end

namf='windy_conv.bin';
dvdy=-1.e-6*5.e+3; %- uniform convergence: wWind = 10^-6 m/s
yy=yv-ymid; vwind=dvdy*yy; vwind(1)=0;
fld=ones(nx,1)*vwind;
if kwr > 0,
 fprintf('write to file: %s\n',namf);
 fid=fopen(namf,'w','b'); fwrite(fid,fld,'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='ice0_area.bin'; iceC0=1.;
iceConc=zeros(nx,ny);
iceConc(21:60,:)=iceC0;
iceConc(:,1)=0; iceConc(:,ny)=0;
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=70; dswA=30;
ymid=(yc(2)+yc(ny))/2;
ysn=yc-ymid ; ysn=ysn/(yc(ny)-yc(2)); ysn=sin(pi*ysn);

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

dlw0=270; dlwA=-20;
namf=['dlw_',int2str(dlw0),'y.bin'];
dlw_y=dlw0+dlwA*ysn;
fld=ones(nx,1)*dlw_y;
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; taC=-10;
%ta0=cel2K+taC; 
%- cheapAML works in deg.C
ta0=taC;
ta=ta0*ones(nx,ny);
namf=['tair_',int2str(taC),'.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
taK=cel2K+taC;
tmpbulk = cvapor_fac*exp(-cvapor_exp./taK);
qa0= (rh/100.)*tmpbulk/atmrho ;
qa=qa0*ones(nx,ny);
namf=['qa',int2str(rh),'_',int2str(taC),'.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);
dtx=1; %- dtx = amplitude of relaxing ocean temp variations in X-dir
%to_x=tfreeze + dtx*sin(pi*(1+2*xc'/nx)) + dtx;
to_x=tfreeze + dtx*sin(pi*(2*xc'/nx)) + dtx;
  mnV=min(to_x); MxV=max(to_x); Avr=mean(to_x);
  fprintf(' SST* av,mn,Mx: %9.6f , %9.6f , %9.6f , %9.6f\n',Avr,mnV,MxV,MxV-mnV);
to=repmat(to_x,[1 ny nt]);
namf=['tocn_',int2str(dtx),'x.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)]);

%--

figure(2);clf;
subplot(211)
plot(xc,to_x,'r-'); hold on;
%plot(xc,dewPt-cel2K,'b-');
plot(xc,tfreeze*ones(nx,1),'k-');
hold off;
AA=axis; axis([-nx/2 nx/2 AA(3:4)]);
legend('To','Tfrz');
grid
title('Ocean Temp ^oC');
%title(['del-Temp-X= ',int2str(dtx),' ; RH= ',int2str(rh),' ; Air Temp (^oC)']);
subplot(212)
plot(yc,dsw_y,'r-'); hold on;
plot(yc,dlw_y-200,'b-');
hold off;
%AA=axis; axis([0 ny AA(3:4)]);
axis([0 ny 30 110]);
grid
legend('sw','lw','Location','South');
title('Downward SW and LW (-200) [W/m^2]');

%--

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');

iceVol=iceH0*iceConc;
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');
%--
1	jmc	1.1
2			kwr=0;
3			%kwr=-1;
4			nx=80; ny=42; nr=3; nt=1;
5
6			xc=[1:nx]; xc=xc-mean(xc);
7	jmc	1.2	yc=[1:ny]-.5; ymid=mean(yc);
8			yv=yc-0.5;
9	jmc	1.1
10			%------------------------------------------------------
11
12			windx=10.;
13			H0=-100.;
14
15			namf='channel.bin';
16			depth=H0*ones(nx,ny); depth(:,1)=0.; depth(:,ny)=0.;
17			if kwr > 0,
18			fprintf('write to file: %s\n',namf);
19			fid=fopen(namf,'w','b'); fwrite(fid,depth,'real*8'); fclose(fid);
20			end
21
22	jmc	1.2	namf=['windx_',int2str(windx),'ms.bin'];
23			uwind=windx*ones(nx,ny,nt);
24	jmc	1.1	if kwr > 0,
25			fprintf('write to file: %s\n',namf);
26	jmc	1.2	fid=fopen(namf,'w','b'); fwrite(fid,uwind,'real*8'); fclose(fid);
27			end
28
29			namf='windy_conv.bin';
30			dvdy=-1.e-6*5.e+3; %- uniform convergence: wWind = 10^-6 m/s
31			yy=yv-ymid; vwind=dvdy*yy; vwind(1)=0;
32			fld=ones(nx,1)*vwind;
33			if kwr > 0,
34			fprintf('write to file: %s\n',namf);
35			fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
36	jmc	1.1	end
37
38			%- file name convention: "const_{xx}.bin" <-> uniform value = xx (in percent)
39			namf='const_00.bin';
40			fld=0*ones(nx,ny,nt);
41			if kwr > 0,
42			fprintf('write to file: %s\n',namf);
43			fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
44			end
45
46			namf='const100.bin'; w0=1.;
47			%var=w0*ones(nx,ny);
48			%if kwr > 0,
49			% fprintf('write to file: %s\n',namf);
50			% fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
51			%end
52
53			namf='const+20.bin'; w0=0.2;
54			var=w0*ones(nx,ny);
55			if kwr > 0,
56			fprintf('write to file: %s\n',namf);
57			fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
58			end
59
60			namf='ice0_area.bin'; iceC0=1.;
61			iceConc=zeros(nx,ny);
62			iceConc(21:60,:)=iceC0;
63			iceConc(:,1)=0; iceConc(:,ny)=0;
64			if kwr > 0,
65			fprintf('write to file: %s\n',namf);
66			fid=fopen(namf,'w','b'); fwrite(fid,iceConc,'real*8'); fclose(fid);
67			end
68
69			namf='ice0_heff.bin'; iceH0=0.2;
70			%iceVol=iceConc*iceH0;
71			%if kwr > 0,
72			% fprintf('write to file: %s\n',namf);
73			% fid=fopen(namf,'w','b'); fwrite(fid,iceVol,'real*8'); fclose(fid);
74			%end
75
76			%------------------------------------------------------
77
78			dsw0=70; dswA=30;
79			ymid=(yc(2)+yc(ny))/2;
80			ysn=yc-ymid ; ysn=ysn/(yc(ny)-yc(2)); ysn=sin(pi*ysn);
81
82			namf=['dsw_',int2str(dsw0),'y.bin'];
83			dsw_y=dsw0+dswA*ysn;
84			fld=ones(nx,1)*dsw_y;
85			if kwr > 0,
86			fprintf('write to file: %s\n',namf);
87			fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
88			end
89
90			dlw0=270; dlwA=-20;
91			namf=['dlw_',int2str(dlw0),'y.bin'];
92			dlw_y=dlw0+dlwA*ysn;
93			fld=ones(nx,1)*dlw_y;
94			if kwr > 0,
95			fprintf('write to file: %s\n',namf);
96			fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
97			end
98
99			cel2K=273.15; taC=-10;
100			%ta0=cel2K+taC;
101			%- cheapAML works in deg.C
102			ta0=taC;
103			ta=ta0*ones(nx,ny);
104			namf=['tair_',int2str(taC),'.bin'];
105			if kwr > 0,
106			fprintf('write to file: %s\n',namf);
107			fid=fopen(namf,'w','b'); fwrite(fid,ta,'real*8'); fclose(fid);
108			end;
109
110			cvapor_fac = 640380.000 ;
111			cvapor_exp = 5107.400 ;
112			atmrho = 1.200 ;
113			rh=70; %- specific humid <--> 70.% relative humid
114			taK=cel2K+taC;
115			tmpbulk = cvapor_fac*exp(-cvapor_exp./taK);
116			qa0= (rh/100.)*tmpbulk/atmrho ;
117			qa=qa0*ones(nx,ny);
118			namf=['qa',int2str(rh),'_',int2str(taC),'.bin'];
119			if kwr > 0,
120			fprintf('write to file: %s\n',namf);
121			fid=fopen(namf,'w','b'); fwrite(fid,qa,'real*8'); fclose(fid);
122			end;
123
124			%- salinity
125			sCst=30;
126			so=sCst*ones(nx,ny,nt);
127			namf='socn.bin';
128			%if kwr > 0,
129			% fprintf('write to file: %s\n',namf);
130			% fid=fopen(namf,'w','b'); fwrite(fid,so,'real*8'); fclose(fid);
131			%end;
132
133			muTf = 5.4e-2;
134			tfreeze=-muTf*sCst;
135			fprintf('T-freeze = %10.6f\n',tfreeze);
136			dtx=1; %- dtx = amplitude of relaxing ocean temp variations in X-dir
137			%to_x=tfreeze + dtxsin(pi(1+2*xc'/nx)) + dtx;
138			to_x=tfreeze + dtxsin(pi(2*xc'/nx)) + dtx;
139			mnV=min(to_x); MxV=max(to_x); Avr=mean(to_x);
140			fprintf(' SST* av,mn,Mx: %9.6f , %9.6f , %9.6f , %9.6f\n',Avr,mnV,MxV,MxV-mnV);
141			to=repmat(to_x,[1 ny nt]);
142			namf=['tocn_',int2str(dtx),'x.bin'];
143			if kwr > 0,
144			fprintf('write to file: %s\n',namf);
145			fid=fopen(namf,'w','b'); fwrite(fid,to,'real*8'); fclose(fid);
146			end;
147
148			%-- make some plots to check: ----------------
149
150			hScal=[-1.1 0.1]*abs(H0);
151			figure(1); clf;
152			subplot(211);
153			var=depth;
154			imagesc(xc,yc,var'); set(gca,'YDir','normal');
155			%caxis(hScal);
156			%change_colmap(-1);
157			colorbar;
158			title('Depth [m]');
159
160			subplot(413);
161			var=depth;
162			j1=2;
163			j2=ny/2;
164			j3=j2+1;
165			plot(xc,var(:,j1),'k-')
166			hold on; j=j+1;
167			plot(xc,var(:,j2),'ro-')
168			plot(xc,var(:,j3),'b-')
169			hold off;
170			axis([-nx/2 nx/2 hScal]);
171			grid
172			legend(int2str(j1),int2str(j2),int2str(j3));
173			title('Depth @ j= cst');
174
175			subplot(414);
176			i=nx/2;
177			plot(yc,var(i,:),'k-')
178			axis([0 ny H01.1 -H0.1]);
179			grid
180			title(['Depth @ i=',int2str(i)]);
181
182			%--
183
184			figure(2);clf;
185			subplot(211)
186			plot(xc,to_x,'r-'); hold on;
187			%plot(xc,dewPt-cel2K,'b-');
188			plot(xc,tfreeze*ones(nx,1),'k-');
189			hold off;
190			AA=axis; axis([-nx/2 nx/2 AA(3:4)]);
191			legend('To','Tfrz');
192			grid
193			title('Ocean Temp ^oC');
194			%title(['del-Temp-X= ',int2str(dtx),' ; RH= ',int2str(rh),' ; Air Temp (^oC)']);
195			subplot(212)
196			plot(yc,dsw_y,'r-'); hold on;
197			plot(yc,dlw_y-200,'b-');
198			hold off;
199			%AA=axis; axis([0 ny AA(3:4)]);
200			axis([0 ny 30 110]);
201			grid
202			legend('sw','lw','Location','South');
203			title('Downward SW and LW (-200) [W/m^2]');
204
205			%--
206
207			figure(3);clf;
208			subplot(311)
209			var=iceConc; ccB=[-1 12]/10;
210			imagesc(xc,yc,var'); set(gca,'YDir','normal');
211			caxis(ccB);
212			%change_colmap(-1);
213			colorbar;
214			title('Ice Concentration in Channel');
215
216			iceVol=iceH0*iceConc;
217			subplot(312)
218			var=iceVol; ccB=[-1 12]/50;
219			imagesc(xc,yc,var'); set(gca,'YDir','normal');
220			caxis(ccB);
221			%change_colmap(-1);
222			colorbar;
223			title('Effective ice thickness in Channel');
224
225			subplot(313)
226			var=iceConc(1,:);
227			%plot(yc,var,'b-x'); hold on;
228			semilogy(yc,var,'b-x'); hold on;
229			var=iceVol(1,:);
230			%plot(yc,var,'r-x'); hold off;
231			semilogy(yc,var,'r-x'); hold on;
232			AA=axis; axis([0 ny [0 2]*iceC0]);
233			grid
234			legend('iceC','hEff','Location','South');
235			title('Initial ice in Channel : y-section');
236			%--