/[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.2 - (show annotations) (download)
Thu Dec 27 14:10:03 2012 UTC (11 years, 4 months ago) by gforget
Branch: MAIN
Changes since 1.1: +33 -4 lines
- add fields for thermodynamics only test.

1
2 kwr=1;
3 nx=80; ny=42; nr=3; nt=2;
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='edge100.bin'; w0=1.;
78 var=w0*ones(nx,ny);
79 var(:,2)=0*var(:,end);
80 var(:,3)=0.01*var(:,end-1);
81 var(:,4)=0.1*var(:,end-2);
82 var(:,end)=0*var(:,end);
83 var(:,end-1)=0.01*var(:,end-1);
84 var(:,end-2)=0.1*var(:,end-2);
85 if kwr >0,
86 fprintf('write to file: %s\n',namf);
87 fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
88 end
89
90 namf='edge+20.bin'; w0=0.2;
91 var=w0*ones(nx,ny);
92 var(:,2)=0*var(:,end);
93 var(:,3)=0.01*var(:,end-1);
94 var(:,4)=0.1*var(:,end-2);
95 var(:,end)=0*var(:,end);
96 var(:,end-1)=0.01*var(:,end-1);
97 var(:,end-2)=0.1*var(:,end-2);
98 if kwr >0,
99 fprintf('write to file: %s\n',namf);
100 fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
101 end
102
103 %------------------------------------------------------
104
105
106 dsw0=100;
107 namf=['dsw_',int2str(dsw0),'.bin'];
108 fld=dsw0*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 dlw0=250;
115 namf=['dlw_',int2str(dlw0),'.bin'];
116 fld=dlw0*ones(nx,ny,nt);
117 if kwr > 0,
118 fprintf('write to file: %s\n',namf);
119 fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
120 end
121
122 cel2K=273.15; dtx=4; %- dtx = amplitude of air temp variations in X-dir
123 ta_x=cel2K + dtx*sin(pi*(1+2*xc'/nx));
124 ta=repmat(ta_x,[1 ny nt]);
125 namf=['tair_',int2str(dtx),'x.bin'];
126 if kwr > 0,
127 fprintf('write to file: %s\n',namf);
128 fid=fopen(namf,'w','b'); fwrite(fid,ta,'real*8'); fclose(fid);
129 end;
130
131 cvapor_fac = 640380.000 ;
132 cvapor_exp = 5107.400 ;
133 atmrho = 1.200 ;
134 rh=70; %- specific humid <--> 70.% relative humid
135 tmpbulk = cvapor_fac*exp(-cvapor_exp./ta_x);
136 qa_x = (rh/100.)*tmpbulk/atmrho ;
137 qa=repmat(qa_x,[1 ny nt]);
138 namf=['qa',int2str(rh),'_',int2str(dtx),'x.bin'];
139 if kwr > 0,
140 fprintf('write to file: %s\n',namf);
141 fid=fopen(namf,'w','b'); fwrite(fid,qa,'real*8'); fclose(fid);
142 end;
143
144 %- salinity
145 sCst=30;
146 so=sCst*ones(nx,ny,nt);
147 namf='socn.bin';
148 %if kwr > 0,
149 % fprintf('write to file: %s\n',namf);
150 % fid=fopen(namf,'w','b'); fwrite(fid,so,'real*8'); fclose(fid);
151 %end;
152
153 muTf = 5.4e-2;
154 tfreeze=-muTf*sCst;
155 fprintf('T-freeze = %10.6f\n',tfreeze);
156 to_y=(yc-2)/ny;
157 to_y=tfreeze+0.5-to_y.*to_y;
158 to=repmat(to_y,[nx 1 nt]);
159 namf='tocn.bin';
160 if kwr > 0,
161 fprintf('write to file: %s\n',namf);
162 fid=fopen(namf,'w','b'); fwrite(fid,to,'real*8'); fclose(fid);
163 end;
164
165 %-- make some plots to check: ----------------
166
167 hScal=[-1.1 0.1]*abs(H0);
168 figure(1); clf;
169 subplot(211);
170 var=depth;
171 imagesc(xc,yc,var'); set(gca,'YDir','normal');
172 %caxis(hScal);
173 %change_colmap(-1);
174 colorbar;
175 title('Depth [m]');
176
177 subplot(413);
178 var=depth;
179 j1=2;
180 j2=ny/2;
181 j3=j2+1;
182 plot(xc,var(:,j1),'k-')
183 hold on; j=j+1;
184 plot(xc,var(:,j2),'ro-')
185 plot(xc,var(:,j3),'b-')
186 hold off;
187 axis([-nx/2 nx/2 hScal]);
188 grid
189 legend(int2str(j1),int2str(j2),int2str(j3));
190 title('Depth @ j= cst');
191
192 subplot(414);
193 i=nx/2;
194 plot(yc,var(i,:),'k-')
195 axis([0 ny H0*1.1 -H0*.1]);
196 grid
197 title(['Depth @ i=',int2str(i)]);
198
199 %--
200 dewPt=(qa_x*atmrho)/cvapor_fac;
201 dewPt=-cvapor_exp./log(dewPt);
202
203 figure(2);clf;
204 subplot(211)
205 plot(xc,ta_x-cel2K,'r-'); hold on;
206 plot(xc,dewPt-cel2K,'b-');
207 plot(xc,tfreeze*ones(nx,1),'k-');
208 hold off;
209 AA=axis; axis([-nx/2 nx/2 AA(3:4)]);
210 legend('ta','dew');
211 grid
212 title(['del-Temp-X= ',int2str(dtx),' ; RH= ',int2str(rh),' ; Air Temp (^oC)']);
213 subplot(212)
214 plot(yc,to_y,'b-'); hold on;
215 plot(yc,tfreeze*ones(ny,1),'k-');
216 hold off;
217 AA=axis; axis([0 ny AA(3:4)]);
218 grid
219 title('Ocean Temp ^oC');

  ViewVC Help
Powered by ViewVC 1.1.22