tutorial_plume_on_slope/input/gendata.m

bf8b271b6c Jean*0001 % This is a matlab script that generates the input data
                0002 % variable x resolution
                0003
                0004 prec='real*8';
                0005 ieee='b';
                0006
                0007 % Dimensions of grid
                0008 %nx=80;
                0009 nx=320;
                0010 ny=1;
                0011 nz=60;
                0012 % Nominal depth of model (meters)
                0013 H=200.0;
                0014 % Size of domain
                0015 Lx=6.40e3;
                0016
                0017 % Horizontal resolution (m)
                0018 %Variable resolution
                0019 res1=2*Lx/(3*nx);
                0020 L1=Lx/2;
                0021 L2=Lx - L1;
                0022 n2=nx - (L1/res1);
                0023 res2=L2/n2;
                0024 A=res2 - res1;
                0025 iswitch1=L1/res1;
                0026 width=40;
                0027 dx=zeros(nx,1);
                0028 for i=1:nx
                0029 dx(i) = res1 + 0.5*A*( tanh( (i-iswitch1)/width) + 1);
                0030 %dx(i) = Lx/nx;
                0031 end
                0032
                0033 dy = Lx/nx
                0034
                0035 % Flux
                0036 Qo=200;
                0037
                0038 % Stratification
                0039 gravity=9.81;
                0040 talpha=2.0e-4;
                0041 N2=0.0;
                0042 Tz=N2/(gravity*talpha);
                0043
                0044 dz=H/nz;
                0045 sprintf('delZ = %d * %7.6g,',nz,dz)
                0046
                0047 x=zeros(nx,1);
                0048 x(1) = dx(1);
                0049 for i=2:nx
                0050 x(i)=x(i-1) + dx(i);
                0051 end
                0052 z=-dz/2:-dz:-H;
                0053
                0054 %Tanh function for cooling
                0055 xswitch = 2.50e3 + Lx/2.0;
                0056 qwidth = 0.1e3;
                0057 Q=0.0*rand([nx,ny]);
                0058 for i=1:nx
                0059 Q(i,:) = Q(i,:) + Qo*0.5*(tanh((Lx-x(i)-xswitch)/qwidth) + 1);
                0060 %Q(i,:) = Q(i,:) + Qo*0.5*(tanh((x(i)-xswitch)/qwidth) + 1);
                0061 end
                0062 fid=fopen('Qnet.forcing','w',ieee); fwrite(fid,Q,prec); fclose(fid);
                0063
                0064 % Temperature profile
                0065 Tref=Tz*z-mean(Tz*z);
                0066 [sprintf('Tref =') sprintf(' %8.6g,',Tref)]
                0067 %t = 0.25*rand([nx,ny,nz]);
                0068 t=0.01*rand([nx,ny,nz]);
                0069 for k=1:nz
                0070  t(:,:,k) = t(:,:,k) + Tref(k);
                0071 end
                0072 fid=fopen('T.init','w',ieee); fwrite(fid,t,prec); fclose(fid);
                0073
                0074 % Sloping channel
                0075 % tanh function for slope
                0076 slope=0.15
                0077 offset=1.5e3 + Lx/2.0;
                0078 dmax=-40.0;
                0079 h1=dmax;
                0080 h2=-H;
                0081 hdiff=(h1-h2);
                0082 xwidth=hdiff/(2.0*slope);
                0083 d=0.0*rand([nx,ny]);
                0084 for i=1:nx
                0085 for j=1:ny
                0086 %d(i,j) = hdiff/2*( exp((x(i)-offset)/xwidth) - exp(-(x(i)-offset)/xwidth))/ ( exp((x(i)-offset)/xwidth) + exp(-(x(i)-offset)/xwidth)) + hdiff/2 -H;
                0087 d(i,j) = hdiff/2*(tanh((Lx-x(i)-offset)/xwidth) + 1) - H;
                0088 end
                0089 end
                0090 d(1,:)=0.0;
                0091 fid=fopen('topog.slope','w',ieee); fwrite(fid,d,prec); fclose(fid);
                0092 plot(x,d(:,1))
                0093 fid=fopen('dx.bin','w',ieee); fwrite(fid,dx,prec); fclose(fid);
                0094