Back to home page

MITgcm

 
 

    


Warning, /verification/exp4/input/gendata.m is written in an unsupported language. File is not indexed.

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
688d11fba8 Alis*0001 % This is a matlab script that generates the input data
                0002 
                0003 % Dimensions of grid
                0004 nx=80;
                0005 ny=42;
                0006 nz=8;
                0007 % Nominal depth of model (meters)
                0008 H=4500;
                0009 % Scale of bump (m)
                0010 L=25e3;
                0011 % Height of bump (m)
                0012 dh=0.90*H;
                0013 % Horizontal resolution (m)
                0014 dx=5e3;
                0015 % Rotation
                0016 f=1e-4;
                0017 % Stratification
                0018 N=1.5 * f*L/H;
                0019 
                0020 % Gravity
                0021 g=9.81;
                0022 % E.O.S.
                0023 alpha=2.e-4;
                0024 
6c607fc2ca Jean*0025 Tz=N^2/(g*alpha);
                0026 fprintf(' Tz= %e ;',Tz);
688d11fba8 Alis*0027 
                0028 dz=H/nz;
6c607fc2ca Jean*0029 fprintf(' delZ = %d * %7.6g\n',nz,dz);
688d11fba8 Alis*0030 
                0031 x=(1:nx)*dx;x=x-mean(x);
                0032 y=(1:ny)*dx;y=y-mean(y);
                0033 z=-dz/2:-dz:-H;
                0034 
                0035 [Y,X]=meshgrid(y,x);
                0036 
                0037 % Temperature profile
6c607fc2ca Jean*0038 fprintf('Tref ='); fprintf(' %8.6g,',Tz*z-mean(Tz*z)); fprintf('\n');
688d11fba8 Alis*0039 
ab42872a05 Alis*0040 ieee='b';
6c607fc2ca Jean*0041 prec='real*8';
ab42872a05 Alis*0042 
688d11fba8 Alis*0043 % Gaussian bump
                0044 h=-H+dh*exp( -(X.^2+Y.^2)/(2*(L^2)) );
6c607fc2ca Jean*0045 fid=fopen('topog.bump','w',ieee); fwrite(fid,h,prec); fclose(fid);
688d11fba8 Alis*0046 
bd2e7b4f6b Mart*0047 % $$$ % Side walls + bump
                0048 % $$$ h(:,1)=0;
                0049 % $$$ h(:,ny)=0;
6c607fc2ca Jean*0050 % $$$ fid=fopen('topog.bumpchannel','w',ieee); fwrite(fid,h,prec); fclose(fid);
bd2e7b4f6b Mart*0051 
                0052 % $$$ % Simple channel
                0053 % $$$ h(:,1)=0;
                0054 % $$$ h(:,2:ny-1)=-H;
                0055 % $$$ h(:,ny)=0;
6c607fc2ca Jean*0056 % $$$ fid=fopen('topog.channel','w',ieee); fwrite(fid,h,prec); fclose(fid);
bd2e7b4f6b Mart*0057 
                0058 % initial fields for salinity
                0059 si = 35;
6c607fc2ca Jean*0060 %fid=fopen('S.init','w',ieee); fwrite(fid,si*ones(nx,ny,nz),prec); fclose(fid);
bd2e7b4f6b Mart*0061 
                0062 % open boundary conditions;
                0063 u0 = .25;
                0064 s0 = si+1;
c514ede061 Jean*0065 w0= 1.e-3;
bd2e7b4f6b Mart*0066 
                0067 % create two time slabs for testing
c514ede061 Jean*0068 uMerid = cat(3,u0*ones(nx,nz),u0*ones(nx,nz));
                0069 uZonal = cat(3,u0*ones(ny,nz),  zeros(ny,nz));
bd2e7b4f6b Mart*0070 sZonal = cat(3,s0*ones(ny,nz),s0*ones(ny,nz));
                0071 
c514ede061 Jean*0072 %- time varying fraction = 1 % of full velocity
                0073 du=u0*0.01;
                0074 uWest = cat(3,(u0+du)*ones(ny,nz),(u0-du)*ones(ny,nz));
                0075 uEast = cat(3,(u0-du)*ones(ny,nz),(u0+du)*ones(ny,nz));
                0076 
                0077 % to test Non-Hydrostatic OBCS:
                0078 w1=[0:nz-1]*pi/nz; w1=-w0*sin(w1);
                0079 wZonal = cat (3, ones(ny,1)*w1, zeros(ny,nz));
                0080 
0baf2dd287 Jean*0081 % to test prescribing Eta in NonLin-FreeSurf formulation:
                0082 et1=0.1;
                0083 etWest = cat(2,+et1*ones(ny,1),-et1*ones(ny,1));
                0084 etEast = cat(2,-et1*ones(ny,1),+et1*ones(ny,1));
                0085 
c514ede061 Jean*0086  fid=fopen('OBmeridU.bin','w',ieee); fwrite(fid,uMerid,prec); fclose(fid);
                0087 %fid=fopen('OBzonalU.bin','w',ieee); fwrite(fid,uZonal,prec); fclose(fid);
                0088  fid=fopen('OB_WestU.bin','w',ieee); fwrite(fid,uWest ,prec); fclose(fid);
                0089  fid=fopen('OB_EastU.bin','w',ieee); fwrite(fid,uEast ,prec); fclose(fid);
                0090  fid=fopen('OBzonalS.bin','w',ieee); fwrite(fid,sZonal,prec); fclose(fid);
                0091  fid=fopen('OBzonalW.bin','w',ieee); fwrite(fid,wZonal,prec); fclose(fid);
0baf2dd287 Jean*0092  fid=fopen('OB_WestH.bin','w',ieee); fwrite(fid,etWest,prec); fclose(fid);
                0093  fid=fopen('OB_EastH.bin','w',ieee); fwrite(fid,etEast,prec); fclose(fid);
6c607fc2ca Jean*0094 
                0095 %- rbcs mask & restauring tracer field:
                0096 msk=ones(nx,ny,nz);
                0097 xMx=max(x);
                0098 shapeX=(x-xMx)/dx;
                0099 shapeX=exp(shapeX*2/3);
                0100 
                0101 [I]=find(shapeX < 5.e-3); fprintf('zero out rbc-mask up to i= %i\n',max(I));
                0102 shapeX(I)=0.;
                0103 var=shapeX'*ones(1,ny*nz);
                0104 fid=fopen('rbcs_mask.bin','w',ieee); fwrite(fid,var,prec); fclose(fid);
                0105 
                0106 tr1=(si+s0)/2;
                0107 var=tr1*ones(nx,ny,nz);
                0108 fid=fopen('rbcs_Tr1_fld.bin','w',ieee); fwrite(fid,var,prec); fclose(fid);
                0109