Warning, /verification/front_relax/input.top/gendata.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 8233d0ce on 2022-02-13 17:14:26 UTC
8233d0ceb9 Jean*0001
0002 dzL=20;
0003 hup=160; kL=round(hup/dzL);
0004 N=15;
0005 k=[1:N];
0006 for k=1:N,
0007 %dd(k)=1.20436^k;
0008 dd(k)=1.1997^k;
0009 end
0010 dd=dd*dzL;
0011 dd=round(dd);
0012 tot=cumsum(dd);
0013 fprintf(' zF - %i :',hup);
0014 %fprintf(' %5.0f',tot(end-5:end));fprintf('\n');
0015 fprintf(' %4.0f',tot(1:end));fprintf('\n');
0016 fprintf(' variable dz:');
0017 fprintf(' %3.0f,',dd);fprintf('\n');
0018
0019 %- uniform dz=20 m down to "hup", then increasing exponentially
0020 % + set remaining (up to nr=25) thickness (near bottom) uniformly
0021 dM=dd(end); kD=10-kL;
0022 dz=[ones(1,kL)*dzL dd ones(1,kD)*dM]; nz=size(dz,2);
0023 zF=[0 -cumsum(dz)]; zC=(zF(1:nz)+zF(2:nz+1))/2;
0024 fprintf(' nz= %i , Full depth: %6.1f\n',nz,-zF(end));
0025
0026 dy=1.e+5; ny=32;
0027 yF=[0:ny]-ny/2; yF=yF*dy; yC=(yF(1:ny)+yF(2:ny+1))/2;
0028
0029 %- set initial temperature:
0030 N2=4.e-6;
0031 alphaT=2.e-4;
0032 gravity=9.81;
0033 fo=1.e-4;
0034 %- N2= db/dz = g alphaT dT/dz
0035 dTdz=N2/gravity/alphaT;
0036
0037 slope=1e-3;
0038 %- jmc: use a smaller slope and slower decay with depth (increase Ho)
0039 slope=6e-4;
0040 dTdy=-slope*dTdz;
0041
0042 fprintf(' N= %9.3e ; g= %5.3f ; alphaT= %9.3e ; dTdz= %e\n', ...
0043 sqrt(N2),gravity,alphaT,dTdz);
0044
0045 Ly=24.*dy; yM=0; Ho=2205;
0046 Ho=2400;
0047 T=zeros(ny,nz);
0048 km=1;
0049 for k=km:nz,
0050 T(:,k)=10 +dTdz*zC(k) ...
0051 +dTdy*Ly*(tanh((yC-yM-200.e3)/200.e3)-tanh((yC-yM+200.e3)/200.e3))/2 ...
0052 *exp(-(3*zC(k)/Ho)^2);
0053 end
0054 %- Create a ML with very weak stratif
0055 N2ML=1.e-12;
0056 %- N2= db/dz = g alphaT dT/dz
0057 dTdzML=N2ML/gravity/alphaT;
0058 kmx=3;
0059 if kmx > 1,
0060 TML=T(:,kmx);
0061 for k=kmx-1:-1:1
0062 TML=TML+dzL*dTdzML;
0063 T(:,k)=TML;
0064 end
0065 end
0066
0067 fprintf('- horiz mean T profile:\n');
0068 var=mean(T); i2=0;
0069 i1=i2+1; i2=floor(nz/2);
0070 fprintf(' %4.1f,',var(i1:i2));
0071 fprintf('\n');
0072 i1=i2+1; i2=nz;
0073 fprintf(' %4.1f,',var(i1:i2));
0074 fprintf('\n');
0075
0076 %- Thermal wind
0077 alphaT=2.e-4;
0078 dbdy=zeros(ny+1,nz);
0079 dbdy(2:ny,:)=gravity*alphaT*(T([2:ny],:)-T([1:ny-1],:))/dy;
0080 ug=zeros(ny+1,nz+1);
0081 for k=nz:-1:1;
0082 ug(:,k) = ug(:,k+1) - dbdy(:,k)*dz(k)/fo ;
0083 end
0084 uCg=ug([1:ny],[1:nz])+ug([2:ny+1],[2:nz+1]) ...
0085 +ug([2:ny+1],[1:nz])+ug([1:ny],[2:nz+1]);
0086 uCg=uCg/4;
0087 %size(uCg)
0088
0089 %- Salt (passive tracer)
0090 S=zeros(ny,nz);
0091 for k=1:nz,
0092 %S(:,k)=ones(ny,1)+zC(k)/Ho; % Linear with z
0093 %S(:,k)=exp(2*zC(k)/Ho); % Exponential with z
0094 S(:,k)=exp(-(2*yC/Ly).^2); % Exponential with y
0095 end
0096 S=10+2*S;
0097
0098 %- add a bump at the top: Exponential shape (width: 2*nyD*dy), height: Hbump
0099 Hbump=50; nyD=5;
0100 ztop=exp(-(yC/dy/nyD).^2); % Exponential width: 2*nyD*dy
0101 zB=ztop(2); ztop=-Hbump*(ztop-zB)/(max(ztop)-zB);
0102 ztop=min(ztop,0);
0103
0104 %- plot to check:
0105
0106 figure(1);clf;
0107 subplot(211);
0108 var=squeeze(T); var(find(var==0))=NaN;
0109 %pcolor(yC,zC,var'); colorbar;
0110 [cs,h]=contour(yC,zC,var',[6:.5:13 13.3]);clabel(cs);
0111 grid
0112 title('t\_ini');
0113
0114 subplot(223);
0115 var=T(1,:);
0116 plot(var,zC,'b-');
0117 var=T(ny/2,:);
0118 hold on;
0119 plot(var,zC,'r-');
0120 i1=ny/2-3; i2=1+ny-i1;
0121 var=T(i1,:);
0122 plot(var,zC,'g-');
0123 var=T(i2,:);
0124 plot(var,zC,'c-');
0125 hold off;
0126 grid
0127 axis([9.4 11 -300 0]);
0128
0129 subplot(224);
0130 var=T(1,:);
0131 plot(var,zC,'b-');
0132 var=T(ny/2,:);
0133 hold on;
0134 plot(var,zC,'r-');
0135 i1=ny/2-3; i2=1+ny-i1;
0136 var=T(i1,:);
0137 plot(var,zC,'g-');
0138 var=T(i2,:);
0139 plot(var,zC,'c-');
0140 hold off;
0141 grid;
0142
0143 figure(2);clf;
0144 subplot(211);
0145 var=squeeze(uCg);
0146 %pcolor(yc,zc,var'); colorbar;
0147 cl=[-10:10]/100;
0148 [cs,h]=contour(yC(1:ny),zC,var',cl);clabel(cs);
0149 grid
0150 title('u\_ini');
0151
0152 subplot(223);
0153 plot(yC,S(:,1),'k-'); axis([yF(1) yF(end) 10 12.1]);
0154 %var=squeeze(S); cl=[-10:10]/10;
0155 %[cs,h]=contour(yC(1:ny),zC,var',cl);clabel(cs);
0156 grid
0157 title('s\_ini');
0158
0159 subplot(426);
0160 plot(yC,ztop,'r-'); axis([yF(1) yF(end) [-1.04 0.04]*50]);
0161 grid
0162 title('top\_bump');
0163
0164 hSlp=[1:ny]*8/(ny-2); hSlp=hSlp-mean(hSlp);
0165 hSlp=zF(end)+300*(1+tanh(hSlp)); hSlp(1)=0; hSlp(ny)=0;
0166 subplot(428);
0167 plot(yC,hSlp,'b-');
0168 AA=axis;
0169 %axis([yF(1) yF(end) AA(3:4)]);
0170 yBnd=[hSlp(2) hSlp(ny-1)]+[-1 1]*100;
0171 axis([yF(1) yF(end) [-26 -18]*100]);
0172 grid
0173 title('topo\_slp');
0174
0175 %return
0176
0177 H=zF(end)*ones(ny,1);
0178 H(1)=0; H(ny)=0;
0179 %fid=fopen('topog.bin','w','b'); fwrite(fid,H,'real*8'); fclose(fid);
0180
0181 fid=fopen('topo_slp.bin','w','b'); fwrite(fid,hSlp,'real*8'); fclose(fid);
0182
0183 fid=fopen('top_bump.bin','w','b'); fwrite(fid,ztop,'real*8'); fclose(fid);
0184
0185 fid=fopen('t_ini.bin','w','b'); fwrite(fid,T,'real*8'); fclose(fid);
0186
0187 fid=fopen('u_ini.bin','w','b'); fwrite(fid,uCg,'real*8'); fclose(fid);
0188
0189 fid=fopen('s_ini.bin','w','b'); fwrite(fid,S,'real*8'); fclose(fid);
0190
0191 return
