/[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.3 - (show annotations) (download)
Wed Jan 9 22:10:20 2013 UTC (11 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64o, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64n, checkpoint64e, checkpoint64d, checkpoint64c, checkpoint64g, checkpoint64f, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l
Changes since 1.2: +54 -28 lines
- change to only one time-record in files (previously, some had 2);
- rename (edge100.bin,edge+20.bin ->  ice0_area.bin ice0_heff.bin) and
  narrow band of low ice-concentration near N & S edge to just 2 rows (instead of 3)
- re-adjust relaxation SST.

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

  ViewVC Help
Powered by ViewVC 1.1.22