Warning, /verification/offline_exf_seaice/input/gendata.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 39c65383 on 2014-03-07 22:31:46 UTC
961c26b165 Jean*0001
e5bcb024b8 Jean*0002 kwr=1; kprt=0;
0fa1f2ea5a Jean*0003 nx=80; ny=42; nr=3; nt=1;
961c26b165 Jean*0004
0005 xc=[1:nx]; xc=xc-mean(xc);
0006 yc=[1:ny]-.5;
0007
0008 %------------------------------------------------------
0009
0010 windx=10.;
0011 H0=-100.;
0012
0013 namf='channel.bin';
0014 depth=H0*ones(nx,ny); depth(:,1)=0.;
f2108654b1 Gael*0015 if kwr > 0,
0016 fprintf('write to file: %s\n',namf);
0017 fid=fopen(namf,'w','b'); fwrite(fid,depth,'real*8'); fclose(fid);
0018 end
961c26b165 Jean*0019
0020 namf='bathy_3c.bin';
0021 msk=abs(xc)'*ones(1,ny)+ones(nx,1)*yc;
0022 depth=H0*ones(nx,ny); depth(:,1)=0.;
0023 depth(find(msk < 24))=0.;
0024 y2d=ones(nx,1)*yc;
0025 depth(find(y2d > ny/2))=H0;
0026 if kwr > 0,
0027 fprintf('write to file: %s\n',namf);
0028 fid=fopen(namf,'w','b'); fwrite(fid,depth,'real*8'); fclose(fid);
0029 end
0030
0031 namf='windx.bin';
0032 wnd=windx*ones(nx,ny,nt);
0033 if kwr > 0,
0034 fprintf('write to file: %s\n',namf);
0035 fid=fopen(namf,'w','b'); fwrite(fid,wnd,'real*8'); fclose(fid);
0036 end
0037
0038 %- file name convention: "const_{xx}.bin" <-> uniform value = xx (in percent)
0039 namf='const_00.bin';
0040 fld=0*ones(nx,ny,nt);
0041 if kwr > 0,
0042 fprintf('write to file: %s\n',namf);
0043 fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
0044 end
0045
0046 namf='const100.bin'; w0=1.;
0047 var=w0*ones(nx,ny);
0048 if kwr > 0,
0049 fprintf('write to file: %s\n',namf);
0050 fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
0051 end
0052
0053 namf='const+20.bin'; w0=0.2;
0054 var=w0*ones(nx,ny);
0055 if kwr > 0,
0056 fprintf('write to file: %s\n',namf);
0057 fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
0058 end
0059
0060 %namf='const+40.bin';
0061 %u0=0.4;
0062 %var=u0*ones(nx,ny);
0063 %if kwr > 0,
0064 % fprintf('write to file: %s\n',namf);
0065 % fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
0066 %end
0067
0068 %namf='const-10.bin';
0069 %v0=-0.1;
0070 %var=v0*ones(nx,ny);
0071 %if kwr > 0,
0072 % fprintf('write to file: %s\n',namf);
0073 % fid=fopen(namf,'w','b'); fwrite(fid,var,'real*8'); fclose(fid);
0074 %end
0075 %------------------------------------------------------
0076
0fa1f2ea5a Jean*0077 namf='ice0_area.bin'; iceC0=1.;
0078 iceConc=iceC0*ones(nx,ny); iceConc(:,1)=0;
0079 iceConc(:,2)=0.00*iceC0;
0080 iceConc(:,3)=0.10*iceC0;
0081 iceConc(:,end) =0.00*iceC0;
0082 iceConc(:,end-1)=0.01*iceC0;
0083 if kwr > 0,
f2108654b1 Gael*0084 fprintf('write to file: %s\n',namf);
0fa1f2ea5a Jean*0085 fid=fopen(namf,'w','b'); fwrite(fid,iceConc,'real*8'); fclose(fid);
f2108654b1 Gael*0086 end
0087
0fa1f2ea5a Jean*0088 namf='ice0_heff.bin'; iceH0=0.2;
0089 iceVol=iceConc*iceH0;
0090 if kwr > 0,
f2108654b1 Gael*0091 fprintf('write to file: %s\n',namf);
0fa1f2ea5a Jean*0092 fid=fopen(namf,'w','b'); fwrite(fid,iceVol,'real*8'); fclose(fid);
f2108654b1 Gael*0093 end
0094
0095 %------------------------------------------------------
0096
961c26b165 Jean*0097 dsw0=100;
0098 namf=['dsw_',int2str(dsw0),'.bin'];
0099 fld=dsw0*ones(nx,ny,nt);
0100 if kwr > 0,
0101 fprintf('write to file: %s\n',namf);
0102 fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
0103 end
0104
0fa1f2ea5a Jean*0105 dlw0=250;
961c26b165 Jean*0106 namf=['dlw_',int2str(dlw0),'.bin'];
0107 fld=dlw0*ones(nx,ny,nt);
0108 if kwr > 0,
0109 fprintf('write to file: %s\n',namf);
0110 fid=fopen(namf,'w','b'); fwrite(fid,fld,'real*8'); fclose(fid);
0111 end
0112
0113 cel2K=273.15; dtx=4; %- dtx = amplitude of air temp variations in X-dir
0114 ta_x=cel2K + dtx*sin(pi*(1+2*xc'/nx));
0115 ta=repmat(ta_x,[1 ny nt]);
0116 namf=['tair_',int2str(dtx),'x.bin'];
0117 if kwr > 0,
0118 fprintf('write to file: %s\n',namf);
0119 fid=fopen(namf,'w','b'); fwrite(fid,ta,'real*8'); fclose(fid);
0120 end;
0121
0122 cvapor_fac = 640380.000 ;
0123 cvapor_exp = 5107.400 ;
0124 atmrho = 1.200 ;
0125 rh=70; %- specific humid <--> 70.% relative humid
0126 tmpbulk = cvapor_fac*exp(-cvapor_exp./ta_x);
0127 qa_x = (rh/100.)*tmpbulk/atmrho ;
0128 qa=repmat(qa_x,[1 ny nt]);
0129 namf=['qa',int2str(rh),'_',int2str(dtx),'x.bin'];
0130 if kwr > 0,
0131 fprintf('write to file: %s\n',namf);
0132 fid=fopen(namf,'w','b'); fwrite(fid,qa,'real*8'); fclose(fid);
0133 end;
0134
0135 %- salinity
0136 sCst=30;
0137 so=sCst*ones(nx,ny,nt);
0138 namf='socn.bin';
0139 %if kwr > 0,
0140 % fprintf('write to file: %s\n',namf);
0141 % fid=fopen(namf,'w','b'); fwrite(fid,so,'real*8'); fclose(fid);
0142 %end;
0143
0144 muTf = 5.4e-2;
0fa1f2ea5a Jean*0145 tfreeze=-muTf*sCst;
961c26b165 Jean*0146 fprintf('T-freeze = %10.6f\n',tfreeze);
0fa1f2ea5a Jean*0147 %- parabolic profile in Y, max @ j=4, min @ j=ny, amplitude=1.K
0148 to_y=(yc-3.5)/(ny-4);
961c26b165 Jean*0149 to_y=tfreeze+0.5-to_y.*to_y;
0fa1f2ea5a Jean*0150 mnV=min(to_y); MxV=max(to_y); Avr=mean(to_y(2:end));
0151 fprintf(' SST* av,mn,Mx: %9.6f , %9.6f , %9.6f , %9.6f\n',Avr,mnV,MxV,MxV-mnV);
961c26b165 Jean*0152 to=repmat(to_y,[nx 1 nt]);
0153 namf='tocn.bin';
0154 if kwr > 0,
0155 fprintf('write to file: %s\n',namf);
0156 fid=fopen(namf,'w','b'); fwrite(fid,to,'real*8'); fclose(fid);
0157 end;
0158
0159 %-- make some plots to check: ----------------
0160
0161 hScal=[-1.1 0.1]*abs(H0);
0162 figure(1); clf;
0163 subplot(211);
0164 var=depth;
0165 imagesc(xc,yc,var'); set(gca,'YDir','normal');
0166 %caxis(hScal);
0167 %change_colmap(-1);
0168 colorbar;
0169 title('Depth [m]');
0170
0171 subplot(413);
0172 var=depth;
0173 j1=2;
0174 j2=ny/2;
0175 j3=j2+1;
0176 plot(xc,var(:,j1),'k-')
0177 hold on; j=j+1;
0178 plot(xc,var(:,j2),'ro-')
0179 plot(xc,var(:,j3),'b-')
0180 hold off;
0181 axis([-nx/2 nx/2 hScal]);
0182 grid
0183 legend(int2str(j1),int2str(j2),int2str(j3));
0184 title('Depth @ j= cst');
0185
0186 subplot(414);
0187 i=nx/2;
0188 plot(yc,var(i,:),'k-')
0189 axis([0 ny H0*1.1 -H0*.1]);
0190 grid
0191 title(['Depth @ i=',int2str(i)]);
0192
0193 %--
0194 dewPt=(qa_x*atmrho)/cvapor_fac;
0195 dewPt=-cvapor_exp./log(dewPt);
0196
0197 figure(2);clf;
e5bcb024b8 Jean*0198 subplot(311)
0199 P(1)=plot(xc,ta_x-cel2K,'r-'); hold on;
0200 P(2)=plot(xc,dewPt-cel2K,'b-');
0201 P(3)=plot(xc,tfreeze*ones(nx,1),'k-');
0202 set(P,'LineWidth',1);
0203 hold off; AA=axis;
0204 axis([-nx/2 nx/2 AA(3:4)]);
961c26b165 Jean*0205 legend('ta','dew');
0206 grid
e5bcb024b8 Jean*0207 xlabel('X')
0208 title(['Air Temp (^oC): del-Temp-X = ',int2str(dtx),' , RH= ',int2str(rh)]);
0209 subplot(312)
0210 P(1)=plot(yc,to_y,'b-'); hold on;
0211 P(2)=plot(yc,tfreeze*ones(ny,1),'k-');
0212 set(P,'LineWidth',1);
0213 hold off; AA=axis;
0214 L=line([1 1],AA(3:4)); set(L,'LineWidth',2.,'Color',[0 0 0]);
0215 axis([0 ny AA(3:4)]);
961c26b165 Jean*0216 grid
e5bcb024b8 Jean*0217 xlabel('Y')
961c26b165 Jean*0218 title('Ocean Temp ^oC');
0fa1f2ea5a Jean*0219
e5bcb024b8 Jean*0220 subplot(313)
0221 var=iceConc(1,:);
0222 P(1)=semilogy(yc,var,'b-x'); hold on;
0223 %plot(yc,var,'b-x'); hold on;
0224 var=iceVol(1,:);
0225 P(2)=semilogy(yc,var,'r-x');
0226 %plot(yc,var,'r-x');
0227 set(P,'LineWidth',1);
0228 hold off; AA=axis;
0229 L=line([1 1],AA(3:4)); set(L,'LineWidth',2.,'Color',[0 0 0]);
0230 axis([0 ny [0 2]*iceC0]);
0231 grid
0232 xlabel('Y')
0233 legend('iceC','hEff','Location','South');
0234 title('Initial ice in Channel : y-section');
0235 %-----
0236 if kprt == 1, f=2;
0237 namfig=sprintf('forcing_%2.2i',f);
0238 fprintf([' print fig= %2i to file: ',namfig,' '],f);
0239 set(f,'PaperOrientation','portrait')
0240 %set(f,'PaperPosition',[0.25 1.5 6. 8.]);
0241 set(f,'PaperPosition',[0.25 1.5 5.25 7.]);
0242 print(f,'-depsc2',namfig); fprintf('\n');
0243 end
0244 %-----
0fa1f2ea5a Jean*0245
0246 figure(3);clf;
0247 subplot(311)
0248 var=iceConc; ccB=[-1 12]/10;
0249 imagesc(xc,yc,var'); set(gca,'YDir','normal');
0250 caxis(ccB);
0251 %change_colmap(-1);
0252 colorbar;
0253 title('Ice Concentration in Channel');
0254
0255 subplot(312)
0256 var=iceVol; ccB=[-1 12]/50;
0257 imagesc(xc,yc,var'); set(gca,'YDir','normal');
0258 caxis(ccB);
0259 %change_colmap(-1);
0260 colorbar;
0261 title('Effective ice thickness in Channel');
0262
0263 subplot(313)
0264 var=iceConc(1,:);
0265 %plot(yc,var,'b-x'); hold on;
0266 semilogy(yc,var,'b-x'); hold on;
0267 var=iceVol(1,:);
0268 %plot(yc,var,'r-x'); hold off;
0269 semilogy(yc,var,'r-x'); hold on;
0270 AA=axis; axis([0 ny [0 2]*iceC0]);
0271 grid
0272 legend('iceC','hEff','Location','South');
0273 title('Initial ice in Channel : y-section');
e5bcb024b8 Jean*0274 %-----
0275
0276 return
