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