Back to home page

MITgcm

 
 

    


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