front_relax/input/gendata.m

fcf1b050b0 Alis*0001 % This is a matlab script that generates the input data
                0002 prec='real*8';
                0003 ieee='b';
5b172de0d2 Jean*0004 kwr=1;
                0005
                0006 %- pick simple number (easier to compare P & Z coords):
                0007 gravity=10; rhoC=1000.;
fcf1b050b0 Alis*0008
                0009 % Dimensions/resolution of grid (resolution/depth in m)
                0010 nx=1; dx=10e3;
                0011 ny=32; dy=10e3;
5b172de0d2 Jean*0012 %- set original level thickness:
fcf1b050b0 Alis*0013 dz=[50 50 55 60 65 70 80 95 120 155 200 260 320 400 480];
5b172de0d2 Jean*0014
fcf1b050b0 Alis*0015 nz=prod(size(dz));
5b172de0d2 Jean*0016 %- add 10 more levels:
                0017 nr=nz+10;
fcf1b050b0 Alis*0018
                0019 % Nominal depth of model (meters)
                0020 Ho=sum(dz);
                0021
5b172de0d2 Jean*0022 % Vertical Coordinates
                0023 zf=[0 -cumsum(dz)];
                0024 zc=(zf(1:end-1)+zf(2:end))/2;
                0025 drF=dz;
                0026 %- after defining centered position (zc), get vertical spacing in between (to use if
                0027 %  we want the interface to be @ the middle beetween 2 center --> specifying delRc=)
                0028 delRc=zeros(1,nz+1);
                0029 delRc(2:nz)=zc(1:nz-1)-zc(2:nz);
                0030 delRc(1)=zf(1)-zc(1); delRc(nz+1)=zc(nz)-zf(nz+1);
                0031 % need to adjust nz+1 value since we will add more levels there:
                0032  delRc(nz+1)=2*delRc(nz+1);
                0033 fprintf(' delRc='); fprintf(' %4.1f,',delRc); fprintf('\n');
                0034 %- set "zfc" assuming interface are @ middle:
                0035 zfc=zf; zfc(2:nz)=(zc(1:nz-1)+zc(2:nz))*0.5;
                0036 drF=zfc(1:nz)-zfc(2:nz+1);
                0037
                0038 % Topography (channel):
fcf1b050b0 Alis*0039 H=-Ho*ones(nx,ny);
                0040 H(:,end)=0;  % Solid wall in North
5b172de0d2 Jean*0041 namf='bathy_inZ.bin';
                0042 if kwr > 0,
                0043  fprintf('write to file: %s\n',namf);
                0044  fid=fopen(namf,'w',ieee); fwrite(fid,H,prec); fclose(fid);
                0045 end
fcf1b050b0 Alis*0046
                0047 % Size of domain (m)
                0048 Lx=dx*nx;
                0049 Ly=dy*(ny-1);  % Solid wall in North
                0050
74b90434ef Alis*0051 % Variable resolution
                0052 y=(1:ny)/(ny-0)-0.5;
5b172de0d2 Jean*0053 dyF=1-0.3*exp( -(5*y).^2 );
                0054 %dyF=ones(1,ny);              % Constant resolution
                0055 dyF=dyF/sum(dyF(1:ny-1))*Ly;
                0056 namf='dy.bin';
                0057 if kwr > 0,
                0058  fprintf('write to file: %s\n',namf);
                0059  fid=fopen(namf,'w',ieee); fwrite(fid,dyF,prec); fclose(fid);
                0060 end
74b90434ef Alis*0061
fcf1b050b0 Alis*0062 % Coordinates
                0063 xc=((1:nx)-0.5)*dx;
5b172de0d2 Jean*0064 yf=-Ly/2+[0 cumsum(dyF)];
74b90434ef Alis*0065 yc=(yf(1:end-1)+yf(2:end))/2;
fcf1b050b0 Alis*0066 [X,Y]=ndgrid(xc,yc);
                0067
                0068 % Stratification
                0069 fo=1e-4;
                0070 alpha=2e-4;
5b172de0d2 Jean*0071 %gravity=9.81; has been set @ the top
fcf1b050b0 Alis*0072 N_over_f=20;
                0073 N2=(N_over_f*fo)^2;
                0074 dTdz=N2/(alpha*gravity);
                0075
                0076 slope=1e-3;
                0077 dTdy=-slope*dTdz;
                0078
                0079 for k=1:nz,
                0080  T(:,:,k)=20 +dTdz*zc(k) ...
                0081              +dTdy*Ly*sin(pi*Y/Ly)...
                0082              *exp(-(3*zc(k)/Ho)^2);
                0083 end
5b172de0d2 Jean*0084 dTy=(T(1,[2:ny],:)-T(1,[1:ny-1],:));
fcf1b050b0 Alis*0085 T(:,end,:)=0;
0837799444 Jean*0086 %- add 10 more levels:
5b172de0d2 Jean*0087 t2d=zeros(nx,ny,nr); t2d(:,:,[1:nz])=T;
                0088 namf='Tini_+10l.bin';
                0089 if kwr > 0,
                0090  fprintf('write to file: %s\n',namf);
                0091  fid=fopen(namf,'w',ieee); fwrite(fid,t2d,prec); fclose(fid);
                0092 end
                0093
                0094 %- compute local N2:
                0095 dTz=(T(1,:,[1:nz-1])-T(1,:,[2:nz]));
                0096 dTz=squeeze(dTz)./repmat(delRc(2:nz),[ny 1]);
                0097 locN2=dTz(1:ny-1,:)*(alpha*gravity); MxV=max(locN2(:));
                0098 fprintf(' max(N2) = %6.4e ; max(c1) = %6.4f\n',MxV,Ho*sqrt(MxV)/pi);
                0099
                0100 %- Thermal wind
                0101 dyC=zeros(ny+1,1); dyC(1)=dyF(1); dyC(ny+1)=dyF(ny);
                0102 dyC(2:ny,1)=(dyF(2:ny)+dyF(1:ny-1))*0.5;
                0103 dbdy=zeros(ny+1,nz);
                0104 dbdy(2:ny,:)=gravity*alpha*squeeze(dTy);
                0105 dbdy=dbdy./repmat(dyC,[1 nz]);
                0106 dbdy(1,:)=dbdy(2,:); dbdy(ny+1,:)=dbdy(ny,:);
                0107 ug=zeros(ny+1,nz+1);
                0108 for k=nz:-1:1;
                0109  ug(:,k) = ug(:,k+1) - dbdy(:,k)*drF(k)/fo ;
                0110 end
                0111 uCg=ug([1:ny],[1:nz])+ug([2:ny+1],[2:nz+1])  ...
                0112    +ug([2:ny+1],[1:nz])+ug([1:ny],[2:nz+1]);
                0113 uCg=uCg/4; uCg(end,:)=0;
                0114 %- add 10 more levels:
                0115 u2d=zeros(nx,ny,nr); u2d(:,:,[1:nz])=reshape(uCg,[1 ny nz]);
                0116 namf='Uini_geos.bin';
                0117 if kwr > 1,
                0118  fprintf('write to file: %s\n',namf);
                0119  fid=fopen(namf,'w',ieee); fwrite(fid,u2d,prec); fclose(fid);
                0120 end
fcf1b050b0 Alis*0121
5b172de0d2 Jean*0122 % Salt (passive tracer), only fct of Y:
                0123 S=zeros(nx,ny,nz);
fcf1b050b0 Alis*0124 for k=1:nz,
                0125 %S(:,:,k)=ones(nx,ny)+zc(k)/Ho;  % Linear with z
                0126 %S(:,:,k)=exp(2*zc(k)/Ho);       % Exponential with z
5b172de0d2 Jean*0127 %S(1,:,k)=exp(2*yc/Ly);           % Exponential with y
                0128  S(1,:,k)=exp(-(2*yc/Ly).^2);     % Exponential with y
                0129 end
                0130 S(:,end,:)=0;
                0131 %- add 10 more levels:
                0132 s2d1=zeros(nx,ny,nr); s2d1(:,:,[1:nz])=S;
                0133 namf='Sini_Ydir.bin';
                0134 if kwr > 0,
                0135  fprintf('write to file: %s\n',namf);
                0136  fid=fopen(namf,'w',ieee); fwrite(fid,s2d1,prec); fclose(fid);
                0137 end
                0138
                0139 % Salt (passive tracer), centered patch:
                0140 % centered location: (yS,zS) ; patch width: (Lh,Lv)
                0141 S=zeros(nx,ny,nz);
                0142  yS=yc(16); Lh=Ly/8;
                0143  zS=zc(10); Lv=500;
                0144 for k=1:nz,
                0145  S(1,:,k)=exp(-((yc-yS)/Lh).^2)*exp(-((zc(k)-zS)/Lv).^2);  % Exponential with y^2 + z^2
fcf1b050b0 Alis*0146 end
                0147 S(:,end,:)=0;
0837799444 Jean*0148 %- add 10 more levels:
5b172de0d2 Jean*0149 s2d2=zeros(nx,ny,nr); s2d2(:,:,[1:nz])=S;
                0150 namf='Sini_Patch.bin';
                0151 if kwr > 0,
                0152  fprintf('write to file: %s\n',namf);
                0153  fid=fopen(namf,'w',ieee); fwrite(fid,s2d2,prec); fclose(fid);
                0154 end
                0155
                0156 figure(1);clf; colormap('jet'); ccB=[0 0];
                0157  subplot(211);
                0158  yax=yc*1.e-3; %- in km
                0159  var=squeeze(T); zax=zc; ccB=[15 21]*1; cInt=[15:0.5:21];
                0160  var=squeeze(t2d(:,:,1:nz)); zax=zc;
                0161 %imagesc(yax,zax,var'); if zax(12) < 0, set(gca,'YDir','normal'); end
                0162 %pcolor(yax,zax,var'); if zax(12) > 0, set(gca,'YDir','reverse'); end
                0163  contourf(yax,zax,var',cInt); if zax(12) > 0, set(gca,'YDir','reverse'); end
                0164  if ccB(2) > ccB(1), caxis(ccB); end
                0165  colorbar;
                0166
                0167  subplot(212);
                0168 %var=squeeze(S); zax=zc; ccB=[-1 11]*0.1; cInt=[0:1:10]*0.1;
                0169 %var=squeeze(s2d1(:,:,1:nz)); ccB=[-1 10]*0.1; cInt=[0:1:20]*5.e-2;
                0170  var=squeeze(s2d2(:,:,1:nz)); ccB=[-1 10]*0.1; cInt=[0:1:20]*5.e-2;
                0171 %yax=yf(1:ny)*1.e-3; %- in km
                0172 %var=squeeze(uCg); zax=zc; ccB=[-1 10]*8.e-3; cInt=[0:1:20]*4.e-3;
                0173 %ccB=[0 0];
                0174 %var(:,:)=NaN; var(1:ny-1,2:nz)=log10(locN2); var(1:ny-1,2:nz)=locN2;
                0175 %imagesc(yax,zax,var'); if zax(12) < 0, set(gca,'YDir','normal'); end
                0176 %pcolor(yax,zax,var'); if zax(12) > 0, set(gca,'YDir','reverse'); end
                0177  contourf(yax,zax,var',cInt); if zax(12) > 0, set(gca,'YDir','reverse'); end
                0178  if ccB(2) > ccB(1), caxis(ccB); end
                0179  colorbar;
                0180
                0181 figure(2); clf;
                0182  yax=[1:ny]; yax=yax-mean(yax);
                0183  subplot(311);
                0184  plot(yax,dyF,'b-');
                0185  AA=axis; AA(1:2)=[-ny/2 ny/2]; axis(AA);
                0186  grid on;
                0187  title('dyF');
                0188
                0189  subplot(613);
                0190  plot(yax,H,'b-');
                0191  AA=axis; AA(1:2)=[-ny/2 ny/2]; axis(AA);
                0192  grid on;
                0193  title('H');
                0194
                0195  subplot(223);
                0196  zax=-[1:nz];
                0197  plot(zc,zax,'kx');
                0198  zax=0.5-[1:nz+1];
                0199  hold on;
                0200  plot(zf,zax,'b-');
                0201  plot(zfc,zax,'r-');
                0202  hold off;
                0203  grid
                0204  title('zc , zf , zfc');
                0205
                0206 %- plot relative difference:
                0207  var=(zfc-zf)./delRc;
                0208  subplot(224);
                0209  plot(var,zax,'b-');
                0210  grid
                0211  title(' (zfc - zf)/dzC');
                0212