/[MITgcm]/MITgcm_contrib/verification_other/offline_cheapaml/input/gendata.m
ViewVC logotype

Contents of /MITgcm_contrib/verification_other/offline_cheapaml/input/gendata.m

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


Revision 1.2 - (show annotations) (download)
Tue Jun 11 01:23:52 2013 UTC (12 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64k, checkpoint64x, checkpoint64z, checkpoint64o, checkpoint64p, checkpoint64r, checkpoint64w, checkpoint64v, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint64m, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65u, checkpoint65j, checkpoint67a, checkpoint67b, checkpoint65i, checkpoint67d, checkpoint65m, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65e, checkpoint64i, checkpoint65h, checkpoint65, checkpoint64j, checkpoint65n, HEAD
Changes since 1.1: +14 -4 lines
- new wind V-component wind (to produce small uniform convergence)
- rename windx.bin --> windx_10ms.bin ; add more diagnostics.

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 yc=[1:ny]-.5; ymid=mean(yc);
8 yv=yc-0.5;
9
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 namf=['windx_',int2str(windx),'ms.bin'];
23 uwind=windx*ones(nx,ny,nt);
24 if kwr > 0,
25 fprintf('write to file: %s\n',namf);
26 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 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 + dtx*sin(pi*(1+2*xc'/nx)) + dtx;
138 to_x=tfreeze + dtx*sin(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 H0*1.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 %--

  ViewVC Help
Powered by ViewVC 1.1.22