Warning, /verification/MLAdjust/input/gendata.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 36741444 on 2014-02-12 23:15:51 UTC
367414446b Jean*0001 %function geninit()
0002 % This function generates initial condition files for the mitgcm.
0003 % if init_vel=1, then the front problem is used
0004 % if init_vel=0, then resting (other than small symmetry-breaking
0005 % random noise) initial conditions are used.
0006
0007 nx=50;
0008 ny=nx/2+1;
0009 nz=40;
0010 init_vel=1;
0011 dxspacing=1000;
0012 dyspacing=dxspacing;
0013
0014 Lx=dxspacing*nx;
0015 Ly=dyspacing*ny;
0016
0017 fprintf(' nx= %i , ny= %i , nz= %i ; dx=%6.1f , dy=%6.1f\n', ...
0018 nx,ny,nz,dxspacing,dxspacing)
0019
0020 %-- Params
0021 g=9.81;
0022 tAlpha=-2e-4;
0023 f0=7.29e-5;
0024 rho=1035;
0025 day=24*60^2;
0026 prec='real*8';
0027 ieee='b';
0028
0029 H=200; %Max depth
0030
0031 fprintf(' Lx=%6.1f , Ly=%6.1f , H=%6.1f\n',Lx,Ly,H);
0032
0033 %-- Grid: x
0034 dx=ones(1,nx); % uniform resolution
0035 dx=dx*Lx/sum(dx);
0036 xf=cumsum([0 dx]); % Face x points
0037 xc=(xf(1:end-1)+xf(2:end))/2; % Centered x points
0038
0039 %-- Grid: y
0040 dy=ones(1,ny); % uniform resolution
0041 dy=dy*Ly/sum(dy);
0042 yf=cumsum([0 dy]); % Face y-points
0043 yc=(yf(1:end-1)+yf(2:end))/2; %Centered y-points
0044 L=yc(end)-yc(1); % this takes into account the wall of topography!!!
0045
0046 %-- Grid: z
0047 dh=H/nz*ones(1,nz);
0048 zf=-round(cumsum([0 dh])/1)*1; % Face z points
0049 dh=-diff(zf);
0050 zc=(zf(1:end-1)+zf(2:end))/2; % centered z points
0051 nz=length(dh);
0052 H=sum(dh);
0053
0054 [XT,YT,ZT]=ndgrid(xc,yc,zc); % This is the centered, temperature grid.
0055 [XU,YU,ZU]=ndgrid(xf,yc,zc); % This is the U grid.
0056 [XV,YV,ZV]=ndgrid(xc,yf,zc); % This is the V grid.
0057 [XW,YW,ZW]=ndgrid(xc,yc,zf); % This is the W grid.
0058 [XB,YB]=ndgrid(xc,yc); % This is the Bathymetry grid.
0059
0060 % Stratification Params
0061 Dml=200; % Mixed-layer Depth
0062 D=Dml;
0063 Lf=2*1e3; % Half-Width of Front
0064 y0=yc(round(length(yc)/2)); % centered location of Front
0065
0066 Ms=(2e-4)^2 ; %db/dy 2e-4 is about 1 K/50 km
0067 Ns=9*Ms^2/f0^2 ; %db/dz
0068
0069 deltatheta=Lf*Ms/g/tAlpha;
0070 T0=17;
0071 dthetadz=-Ns/g/tAlpha;
0072
0073 fprintf(' Ms=%15.6e , Ns=%15.6e , T0= %6.3f\n',Ms,Ns,T0);
0074 fprintf(' deltatheta=%15.6e , dthetadz=%15.6e\n',deltatheta,dthetadz);
0075
0076 %-- Bathymetry
0077 % Bathymetry is on Xc, Yc
0078 hh=ones(size(XB));
0079 hh(:,end)=0*hh(:,end);
0080 hh=-H*hh;
0081
0082 figure(1); clf
0083 subplot(221)
0084 pcolor(XB/1e3,YB/1e3,hh)
0085 axis equal
0086
0087 subplot(223)
0088 plot(xc/1e3,dx/1e3,'.'); xlabel('X (km)');ylabel('\Delta x (km)')
0089 title('\Delta x')
0090
0091 subplot(222)
0092 plot(dy/1e3,yc/1e3,'.'); ylabel('Y (km)');xlabel('\Delta y (km)')
0093 title('\Delta y')
0094
0095 subplot(224)
0096 plot(dh,zc,'.'); ylabel('Z (m)');xlabel('\Delta z (m)')
0097 title('\Delta z')
0098
0099 fid=fopen('topo_sl.bin','w',ieee); fwrite(fid,hh,prec); fclose(fid);
0100
0101 % Initial Temp Profile
0102 figure(2); clf
0103
0104 theta=T0+dthetadz*(ZT-ZT(1,1,end))+deltatheta*tanh((YT-y0)/Lf)/2;
0105
0106 % Impose a strong initial and restoring mixed layer to Dml
0107 i=1;
0108 iml=1;
0109
0110 while ((ZT(1,1,iml)>-Dml)&(iml<length(ZT(1,1,:))))
0111 iml=iml+1;
0112 end
0113
0114 while (i<iml)
0115 theta(:,:,i)=theta(:,:,iml);
0116 i=i+1;
0117 end
0118
0119 fprintf(' Z , Theta: j=1 , j=end\n');
0120 fprintf(' %10.4f %10.4f %10.4f\n', ...
0121 [ZT(1,1,1) theta(1,1,1) theta(1,end,1) ; ...
0122 ZT(1,1,i) theta(1,1,i) theta(1,end,i); ...
0123 ZT(1,1,end) theta(1,1,end) theta(1,end,end)]' ...
0124 );
0125
0126 subplot(3,2,1)
0127 [h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)));
0128 axis([(yc(round(length(yc)/2))-1.5*Lf)/1e3 (yc(round(length(yc)/2))+1.5*Lf)/1e3 -Dml 0])
0129 colorbar
0130 title('Potl Temp')
0131 xlabel('y (km)');ylabel('z (m)')
0132
0133 subplot(3,2,2)
0134 [h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)));
0135 colorbar
0136 title('Potl Temp')
0137 xlabel('y (km)');ylabel('z (m)')
0138
0139
0140 subplot(3,2,3)
0141 dens=theta*rho*tAlpha+rho;
0142 [h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)));
0143 title('Density')
0144 xlabel('y (km)');ylabel('z (m)')
0145 axis([(yc(round(length(yc)/2))-1.5*Lf)/1e3 (yc(round(length(yc)/2))+1.5*Lf)/1e3 -Dml 0])
0146
0147 %Spice
0148
0149 x1=xc(round(length(xc)/4)); % centered location of Front #1
0150 x2=xc(round(3*length(xc)/4)); % centered location of Front #2
0151 spice=T0+dthetadz*(ZT-ZT(1,1,end))+deltatheta*(tanh((XT-x1)/Lf)-tanh((XT-x2)/Lf)-1)/2;
0152
0153 i=1;
0154 while (i<iml)
0155 spice(:,:,i)=spice(:,:,iml);
0156 i=i+1;
0157 end
0158
0159 subplot(3,2,4)
0160 [h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)));
0161 title('Density')
0162 xlabel('y (km)');ylabel('z (m)')
0163
0164 subplot(3,2,5)
0165 [h,c]=contourf(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(spice(:,1,:)));
0166 title('Spice')
0167 xlabel('x (km)');ylabel('z (m)')
0168
0169 subplot(3,2,6)
0170 plot(YB(1,:)/1e3,hh(1,:))
0171 title('topography')
0172 xlabel('x (km)');ylabel('value')
0173
0174 %-- Perturb Initial temperature
0175 pert=rand(size(theta(:,:,1)));
0176 pert=1e-5*(pert-0.5);
0177 for i=1:length(theta(1,1,:))
0178 theta(:,:,i)=theta(:,:,i)+pert;
0179 end
0180
0181 fid=fopen('thetaInitial.bin','w',ieee);
0182 fwrite(fid,theta,prec); fclose(fid);
0183
0184 fid=fopen('spiceInitial.bin','w',ieee);
0185 fwrite(fid,spice,prec); fclose(fid);
0186
