Back to home page

MITgcm

 
 

    


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