%function geninit()
%  This function generates initial condition files for the mitgcm.
% if init_vel=1, then the front problem is used
% if init_vel=0, then resting (other than small symmetry-breaking
% random noise) initial conditions are used.

clear all;  

%-------------------------
%outdir='/Users/sanjiv/CODES/MITgcm/CASES/MLAdjust/BoB/IC/notopo/'; 
outdir='/Users/sanjiv/CODES/MITgcm/CASES/MLAdjust/BoB/IC/NSwall/'; 
%-------------------------
nx=100;
ny=300;
nz=40;
init_vel=1;
dxspacing=400;
dyspacing=dxspacing;

Lx=dxspacing*nx;
Ly=dyspacing*ny;

fprintf(' nx= %i , ny= %i , nz= %i ; dx=%6.1f , dy=%6.1f\n', ...
          nx,ny,nz,dxspacing,dxspacing)

%-- Params
g=9.81;
tAlpha = -3.16e-4; % 16N, BoB (2013 radiator survey) 
sBeta= abs(tAlpha)*2; 
f0=4e-5;
rho=1027;
day=24*60^2;
prec='real*8';
ieee='b';

H=100;            %Max depth

fprintf(' Lx=%6.1f , Ly=%6.1f , H=%6.1f\n',Lx,Ly,H);

%-- Grid: x
dx=ones(1,nx);                                  % uniform resolution
dx=dx*Lx/sum(dx); 
xf=cumsum([0 dx]); % Face x points
xc=(xf(1:end-1)+xf(2:end))/2; % Centered x points

%-- Grid: y
dy=ones(1,ny);                                  % uniform resolution
dy=dy*Ly/sum(dy); 
yf=cumsum([0 dy]);  % Face y-points
yc=(yf(1:end-1)+yf(2:end))/2;  %Centered y-points
L=yc(end)-yc(1);	% this takes into account the wall of topography!!!

%-- Grid: z
dh=H/nz*ones(1,nz);
zf=-round(cumsum([0 dh])/1)*1;   % Face z points
dh=-diff(zf);
zc=(zf(1:end-1)+zf(2:end))/2;  % centered z points
nz=length(dh);
H=sum(dh);

[XT,YT,ZT]=ndgrid(xc,yc,zc); % This is the centered, temperature grid.
[XU,YU,ZU]=ndgrid(xf,yc,zc); % This is the U grid.
[XV,YV,ZV]=ndgrid(xc,yf,zc); % This is the V grid.
[XW,YW,ZW]=ndgrid(xc,yc,zf); % This is the W grid.
[XB,YB]=ndgrid(xc,yc); % This is the Bathymetry grid.

% Stratification Params
Dml=20;             % Mixed-layer Depth
D=Dml;
Lf=2*1e3;          % Half-Width of Front
y0=yc(round(length(yc)/2));     % centered location of Front

Ms=-8e-7;     %db/dy  2e-4 is about 1 K/50 km
%Ns=10*Ms^2/f0^2 ;  %db/dz
Ns=3e-4;

deltatheta= 0.25*Lf*Ms/g/tAlpha;
deltaS    = 0.75*Lf*Ms/g/sBeta; 
T0=20;S0=34.5; 
dthetadz = -0.8*Ns/g/tAlpha;
dSdz     = -0.2*Ns/g/sBeta; 

fprintf(' Ms=%15.6e , Ns=%15.6e , T0= %6.3f, S0=%6.3f\n',Ms,Ns,T0,S0);
fprintf(' deltatheta=%15.6e , dthetadz=%15.6e\n',deltatheta,dthetadz);
fprintf(' deltaS=%15.6e , dSdz=%15.6e\n',deltaS,dSdz);
%-- Bathymetry
% Bathymetry is on Xc, Yc
hh=ones(size(XB));
%--uncomment the two lines to have a wall at the southern/northern boundaries-----
hh(:,end)=0*hh(:,end);
hh(:,1)=0*hh(:,1);
%--------------------------------------------------------------------------------
hh=-H*hh;

figure(1); clf
subplot(221)
pcolor(XB/1e3,YB/1e3,hh); shading interp 
axis equal

subplot(223)
plot(xc/1e3,dx/1e3,'.'); xlabel('X (km)');ylabel('\Delta x (km)')
title('\Delta x')

subplot(222)
plot(dy/1e3,yc/1e3,'.'); ylabel('Y (km)');xlabel('\Delta y (km)')
title('\Delta y')

subplot(224)
plot(dh,zc,'.'); ylabel('Z (m)');xlabel('\Delta z (m)')
title('\Delta z')

fid=fopen(strcat(outdir,'topo_BoB.bin'),'w',ieee); fwrite(fid,hh,prec); fclose(fid);

% Initial Temp Profile
figure(2); clf

zdecay=10; 
theta=T0+dthetadz*(ZT-ZT(1,1,end)) + deltatheta*exp((ZT+Dml)/zdecay).*tanh((YT-y0)/Lf)/2;
salt =S0+    dSdz*(ZT-ZT(1,1,end)) +     deltaS*exp((ZT+Dml)/zdecay).*tanh((YT-y0)/Lf)/2;

% Impose a strong initial and restoring mixed layer to Dml
i=1;
iml=1;

while (ZT(1,1,iml)>-Dml)&(iml<length(ZT(1,1,:)))
  iml=iml+1;
end

while (i<iml)
  theta(:,:,i)=theta(:,:,iml);
  salt(:,:,i) = salt(:,:,iml);
  i=i+1;
end


fprintf('      Z    , Theta: j=1 ,  j=end\n');
fprintf(' %10.4f %10.4f %10.4f\n', ...
[ZT(1,1,1) theta(1,1,1) theta(1,end,1) ; ...
 ZT(1,1,i) theta(1,1,i) theta(1,end,i);  ...
 ZT(1,1,end) theta(1,1,end) theta(1,end,end)]' ...
);

subplot(3,2,1)
[h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)));
axis([(yc(round(length(yc)/2))-1.5*Lf)/1e3 (yc(round(length(yc)/2))+1.5*Lf)/1e3 -Dml 0])
colorbar
title('Potl Temp')
xlabel('y (km)');ylabel('z (m)')

subplot(3,2,2)
[h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)),20);
colorbar
title('Potl Temp')
xlabel('y (km)');ylabel('z (m)')


subplot(3,2,3)
dens=theta*rho*tAlpha+salt*rho*sBeta;
[h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)));
title('Density')
xlabel('y (km)');ylabel('z (m)')
axis([(yc(round(length(yc)/2))-1.5*Lf)/1e3 (yc(round(length(yc)/2))+1.5*Lf)/1e3 -Dml 0])

%Spice
x1=xc(round(length(xc)/4));     % centered location of Front #1
x2=xc(round(3*length(xc)/4));     % centered location of Front #2
%spice=T0+dthetadz*(ZT-ZT(1,1,end))+deltatheta*(tanh((XT-x1)/Lf)-tanh((XT-x2)/Lf)-1)/2;
spice=T0+dthetadz*(ZT-ZT(1,1,end))+deltatheta*exp((ZT+Dml)/zdecay).*(tanh((XT-x1)/Lf)-tanh((XT-x2)/Lf)-1)/2;

i=1;
while (i<iml)
  spice(:,:,i)=spice(:,:,iml);
  i=i+1;
end

subplot(3,2,4)
[h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)));
title('Density')
xlabel('y (km)');ylabel('z (m)')

subplot(3,2,5)
[h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(salt(1,:,:)));
title('Salt')
xlabel('y (km)');ylabel('z (m)')

subplot(3,2,6)
plot(YB(1,:)/1e3,hh(1,:))
title('topography')
xlabel('y (km)');ylabel('value')


%-- Perturb Initial temperature
pert=rand(size(theta(:,:,1)));
pert=1e-5*(pert-0.5);
for i=1:length(theta(1,1,:))
  theta(:,:,i)=theta(:,:,i)+pert;
  salt (:,:,i)=salt(:,:,i)+pert;
end

fid=fopen(strcat(outdir,'theta_BoB.bin'),'w',ieee);
fwrite(fid,theta,prec); fclose(fid);

fid=fopen(strcat(outdir,'salt_BoB.bin'),'w',ieee);
fwrite(fid,salt,prec); fclose(fid);

tauMax=0.05;
tau=tauMax*sin(pi*YB);
fid=fopen(strcat(outdir,'windx_sin.y'),'w',ieee); 
fwrite(fid,tau,prec); fclose(fid);