Back to home page

MITgcm

 
 

    


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