/[MITgcm]/MITgcm/verification/offline_exf_seaice/input/gendata.m
ViewVC logotype

Contents of /MITgcm/verification/offline_exf_seaice/input/gendata.m

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


Revision 1.4 - (show annotations) (download)
Fri Mar 7 22:30:26 2014 UTC (10 years, 2 months ago) by jmc
Branch: MAIN
Changes since 1.3: +47 -15 lines
improve the plots

1
2 kwr=1; kprt=0;
3 kwr=0; kprt=0;
4 nx=80; ny=42; nr=3; nt=1;
5
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 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
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 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 fprintf('write to file: %s\n',namf);
86 fid=fopen(namf,'w','b'); fwrite(fid,iceConc,'real*8'); fclose(fid);
87 end
88
89 namf='ice0_heff.bin'; iceH0=0.2;
90 iceVol=iceConc*iceH0;
91 if kwr > 0,
92 fprintf('write to file: %s\n',namf);
93 fid=fopen(namf,'w','b'); fwrite(fid,iceVol,'real*8'); fclose(fid);
94 end
95
96 %------------------------------------------------------
97
98 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 dlw0=250;
107 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 + dtx*sin(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 tfreeze=-muTf*sCst;
147 fprintf('T-freeze = %10.6f\n',tfreeze);
148 %- parabolic profile in Y, max @ j=4, min @ j=ny, amplitude=1.K
149 to_y=(yc-3.5)/(ny-4);
150 to_y=tfreeze+0.5-to_y.*to_y;
151 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 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 H0*1.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 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 legend('ta','dew');
207 grid
208 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 grid
218 xlabel('Y')
219 title('Ocean Temp ^oC');
220
221 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
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 %-----
276
277 return

  ViewVC Help
Powered by ViewVC 1.1.22