internal_wave/input/gendata.m

fc0c9d8762 Alis*0001 % This is a matlab script that generates the input data
                0002 % variable x resolution
                0003 prec='real*8';
                0004 ieee='b';
                0005
                0006
                0007 % Dimensions of grid
                0008 nx=60;
dbad296d97 Alis*0009 ny=1;
fc0c9d8762 Alis*0010 nz=20;
                0011 % Nominal depth of model (meters)
                0012 H=200.0;
                0013 % Size of domain
                0014 Lx=13.3e3;
                0015
                0016 % Horizontal resolution (m)
                0017 dx=zeros(nx,1);
                0018 for i=1:nx
                0019 dx(i) = Lx/(nx+1);
                0020 end
                0021
                0022 dy = Lx/nx
                0023 % Stratification
                0024 gravity=9.81;
                0025 talpha=2.0e-4;
                0026 N2=1e-6;
                0027
                0028 Tz=N2/(gravity*talpha);
                0029
                0030 dz=H/nz;
                0031 sprintf('delZ = %d * %7.6g,',nz,dz)
                0032
                0033 x=zeros(nx,1);
                0034 x(1) = dx(1);
                0035 for i=2:nx
                0036 x(i)=x(i-1) + dx(i);
                0037 end
                0038 z=-dz/2:-dz:-H;
                0039
                0040 %[Y,X]=meshgrid(y,x);
                0041
                0042 % Temperature profile
                0043 Tref=Tz*z-mean(Tz*z);
                0044 [sprintf('Tref =') sprintf(' %8.6g,',Tref)]
                0045 t=0.0*rand([nx,ny,nz]);
                0046 for k=1:nz
                0047 t(:,:,k) = t(:,:,k) + Tref(k);
                0048 end
                0049 fid=fopen('T.init','w',ieee); fwrite(fid,t,prec); fclose(fid);
                0050
                0051 % Sloping channel
dbad296d97 Alis*0052 slope=0.03
fc0c9d8762 Alis*0053 offset=2.5e3;
                0054 dmax=-40.0;
                0055 d=0.0*rand([nx,ny]);
                0056 for i=1:nx
                0057 for j=1:ny
                0058 d(i,j) = -H;
                0059 d(i,j) = d(i,j) + slope*(x(i) - offset);
                0060 if d(i,j) < -H ;
                0061 d(i,j) = -H;
                0062 end
                0063 if d(i,j) > dmax;
                0064 d(i,j) = dmax;
                0065 end
                0066 end
                0067 end
dbad296d97 Alis*0068 d(nx,:)=0.0;
fc0c9d8762 Alis*0069 fid=fopen('topog.slope','w',ieee); fwrite(fid,d,prec); fclose(fid);
dbad296d97 Alis*0070 plot(x,d(:,1))
fc0c9d8762 Alis*0071
                0072 fid=fopen('delXvar','w',ieee); fwrite(fid,dx,prec); fclose(fid);
                0073
                0074 %convex slope
                0075 nxslope=(dmax + H)/(slope)
                0076 d1=zeros(nx,ny);
                0077 hamp=(H-dmax)/5.0
                0078 pi=4.0*atan(1.0)
                0079 for i=1:nx
                0080 for j=1:ny
                0081 if x(i) < (offset + nxslope)
                0082 if x(i) < offset
                0083 d1(i,j) = d(i,j);
                0084 else
                0085 d1(i,j) = -H;
                0086 d1(i,j) = d(i,j) + hamp*sin(pi*(x(i)-offset)/nxslope);
                0087 if d1(i,j) < -H ;
                0088 d1(i,j) = -H;
                0089 end
                0090 if d1(i,j) > dmax;
                0091 d1(i,j) = dmax;
                0092 end
                0093 end
                0094 else
                0095 d1(i,j) = d(i,j);
                0096 end
                0097 end
                0098 end
                0099 %d1(end-1:end,:)=d1(1:2,:); % debug by aja
                0100 fid=fopen('topog.convex','w',ieee); fwrite(fid,d1,prec); fclose(fid);
dbad296d97 Alis*0101 plot(x,d1(:,1),'g')
fc0c9d8762 Alis*0102
                0103 %convex slope
                0104 d2=zeros(nx,ny);
                0105 for i=1:nx
                0106 for j=1:ny
                0107 if x(i) < (offset + nxslope)
                0108 if x(i) < offset
                0109 d2(i,j) = d(i,j);
                0110 else
                0111 d2(i,j) = -H;
                0112 d2(i,j) = d(i,j) - hamp*sin(pi*(x(i)-offset)/nxslope);
                0113 if d2(i,j) < -H ;
                0114 d2(i,j) = -H;
                0115 end
                0116 if d2(i,j) > dmax;
                0117 d2(i,j) = dmax;
                0118 end
                0119 end
                0120 else
                0121 d2(i,j) = d(i,j);
                0122 end
                0123 end
                0124 end
                0125 %d2(end-1:end,:)=d2(1:2,:); % debug by aja
                0126 fid=fopen('topog.concave','w',ieee); fwrite(fid,d2,prec); fclose(fid);
                0127 hold on
dbad296d97 Alis*0128 plot(x,d2(:,1),'r')
                0129 hold off
fc0c9d8762 Alis*0130
                0131
                0132 fid=fopen('delXvar','w',ieee); fwrite(fid,dx,prec); fclose(fid);