front_relax/input.mxl/gendata.m

463e482235 Jean*0001
                0002 dzL=20;
                0003 hup=200; kL=round(hup/dzL);
                0004 N=15;
                0005 k=[1:N];
                0006 for k=1:N,
                0007  dd(k)=1.20436^k;
                0008 end
                0009 dd=dd*dzL;
                0010 dd=round(dd);
                0011 tot=cumsum(dd)/100;
                0012 fprintf(' %5.2f',tot(end-5:end));fprintf('\n');
                0013 fprintf(' %2.0f.,',dd);fprintf('\n');
                0014
                0015 %- uniform dz=20 m down to 200m, then increasing exponentially
                0016 dz=[ones(1,kL)*dzL dd]; nz=size(dz,2);
                0017 zF=[0 -cumsum(dz)]; zC=(zF(1:nz)+zF(2:nz+1))/2;
5eacf8adec Jean*0018 dy=1.e+5; ny=32;
                0019 yF=[0:ny]-ny/2; yF=yF*dy; yC=(yF(1:ny)+yF(2:ny+1))/2;
463e482235 Jean*0020
                0021 %- set initial temperature:
                0022 N2=4.e-6;
                0023 alphaT=2.e-4;
                0024 gravity=9.81;
                0025 fo=1.e-4;
                0026 %- N2= db/dz = g alphaT dT/dz
                0027 dTdz=N2/gravity/alphaT;
                0028
                0029 slope=1e-3;
                0030 dTdy=-slope*dTdz;
                0031
                0032 fprintf(' N= %9.3e ; g= %5.3f ; alphaT= %9.3e ; dTdz= %e\n', ...
                0033          sqrt(N2),gravity,alphaT,dTdz);
                0034
5eacf8adec Jean*0035 Ly=24.*dy; yM=0; Ho=2205;
                0036 T=zeros(ny,nz);
463e482235 Jean*0037 km=6; kmx=6;
                0038 for k=km:nz,
                0039  T(:,k)=10 +dTdz*zC(k) ...
5eacf8adec Jean*0040            +dTdy*Ly*(tanh((yC-yM-200.e3)/200.e3)-tanh((yC-yM+200.e3)/200.e3))/2 ...
463e482235 Jean*0041            *exp(-(3*zC(k)/Ho)^2);
                0042 end
                0043 %Create a ML
                0044 for k=1:km
                0045  T(:,k)=T(:,km);
                0046 end
                0047 TML=mean(T(:,1:kmx),2);
                0048 for k=1:kmx
                0049  T(:,k)=TML;
                0050 end
                0051
                0052 %-- Thermal wind
                0053 alphaT=2.e-4;
5eacf8adec Jean*0054 dbdy=zeros(ny+1,nz);
                0055 dbdy(2:ny,:)=gravity*alphaT*(T([2:ny],:)-T([1:ny-1],:))/dy;
                0056 ug=zeros(ny+1,nz+1);
463e482235 Jean*0057 for k=nz:-1:1;
5eacf8adec Jean*0058  ug(:,k) = ug(:,k+1) - dbdy(:,k)*dz(k)/fo ;
463e482235 Jean*0059 end
5eacf8adec Jean*0060 uCg=ug([1:ny],[1:nz])+ug([2:ny+1],[2:nz+1])  ...
                0061    +ug([2:ny+1],[1:nz])+ug([1:ny],[2:nz+1]);
                0062 uCg=uCg/4;
                0063 %size(uCg)
463e482235 Jean*0064
                0065 % Salt (passive tracer)
5eacf8adec Jean*0066 S=zeros(ny,nz);
463e482235 Jean*0067 for k=1:nz,
5eacf8adec Jean*0068 %S(:,k)=ones(ny,1)+zC(k)/Ho;  % Linear with z
463e482235 Jean*0069 %S(:,k)=exp(2*zC(k)/Ho);      % Exponential with z
5eacf8adec Jean*0070  S(:,k)=exp(-(2*yC/Ly).^2);   % Exponential with y
463e482235 Jean*0071 end
                0072 S=10+2*S;
                0073
                0074 %- plot to check:
                0075
                0076 figure(1);clf;
                0077 subplot(211);
                0078 var=squeeze(T); var(find(var==0))=NaN;
5eacf8adec Jean*0079 %pcolor(yC,zC,var'); colorbar;
                0080 [cs,h]=contour(yC,zC,var',[6:.5:13 13.3]);clabel(cs);
463e482235 Jean*0081 grid
5eacf8adec Jean*0082 title('t\_ini');
463e482235 Jean*0083
                0084 subplot(223);
                0085 var=T(1,:);
                0086 plot(var,zC,'b-');
                0087 grid
                0088 axis([9.4 10 -300 0]);
                0089 subplot(224);
                0090 plot(var,zC,'b-');
5eacf8adec Jean*0091 var=T(ny/2,:);
463e482235 Jean*0092 hold on;
                0093 plot(var,zC,'r-');
5eacf8adec Jean*0094 i1=ny/2-3; i2=1+ny-i1;
463e482235 Jean*0095 var=T(i1,:);
                0096 plot(var,zC,'g-');
                0097 var=T(i2,:);
                0098 plot(var,zC,'c-');
                0099 hold off;
                0100 grid;
                0101
                0102 figure(2);clf;
                0103 subplot(211);
5eacf8adec Jean*0104 var=squeeze(uCg);
463e482235 Jean*0105 %pcolor(yc,zc,var'); colorbar;
                0106 cl=[-10:10]/100;
5eacf8adec Jean*0107 [cs,h]=contour(yC(1:ny),zC,var',cl);clabel(cs);
463e482235 Jean*0108 grid
5eacf8adec Jean*0109 title('u\_ini');
                0110
                0111 subplot(223);
                0112 plot(yC,S(:,1),'k-'); axis([yF(1) yF(end) 10 12.1]);
                0113 %var=squeeze(S); cl=[-10:10]/10;
                0114 %[cs,h]=contour(yC(1:ny),zC,var',cl);clabel(cs);
463e482235 Jean*0115 grid
5eacf8adec Jean*0116 title('s\_ini');
463e482235 Jean*0117
5eacf8adec Jean*0118 hSlp=[1:ny]*8/(ny-2); hSlp=hSlp-mean(hSlp);
                0119 hSlp=zF(end)+300*(1+tanh(hSlp)); hSlp(1)=0; hSlp(ny)=0;
                0120 subplot(224);
                0121 plot(yC,hSlp,'b-');
                0122 AA=axis;
                0123 axis([yF(1) yF(end) AA(3:4)]);
                0124 grid
                0125 title('topo\_slp');
463e482235 Jean*0126
5eacf8adec Jean*0127 %return
                0128
                0129 H=zF(end)*ones(ny,1);
                0130 H(1)=0; H(ny)=0;
463e482235 Jean*0131 fid=fopen('topog.bin','w','b'); fwrite(fid,H,'real*8'); fclose(fid);
                0132
5eacf8adec Jean*0133 fid=fopen('topo_slp.bin','w','b'); fwrite(fid,hSlp,'real*8'); fclose(fid);
                0134
463e482235 Jean*0135 fid=fopen('t_ini.bin','w','b'); fwrite(fid,T,'real*8'); fclose(fid);
                0136
5eacf8adec Jean*0137 fid=fopen('u_ini.bin','w','b'); fwrite(fid,uCg,'real*8'); fclose(fid);
463e482235 Jean*0138
                0139 fid=fopen('s_ini.bin','w','b'); fwrite(fid,S,'real*8'); fclose(fid);
                0140
                0141 return