%% Generation of input data for ice shelf
%
% Requires an entry in Modelrun.csv for this specific run with setup parameters.
% 
% 1. Create bed topography
% 2. Create ice shelf topography
% 3. oceanographic setting and ice shelf anomaly file
% 4. Create velocity boundary
% 5. Temperature and salinity initial fields
% 6. plot input data
%
% Nadine Steiger 09.01.2020 
%-------------------------------------------------------

clear all; close all


model_path = '/media/nst092/Seagate_Nadine/Modelling/';
table_path = [model_path,'Modelruns.csv'];

%% read setup for experiment
% ----------------------------------------

Exp = readtable(table_path);                       % read table with experiment variables
runnum = 205;                                                                   % experiment number

run('../../setgrid_univers.m')                                                 % run function to set the grid of the domain

L           = Exp.L(ind);
dH          = Exp.dH(ind);
topo        = Exp.topo{ind};        % the bed topography between 'box', 'vshape', 'vshapelab' and 'sinshape';
flow        = 'sideflow';           % Exp.inflowBC{ind};   % chose the inflow at x=0 between 'noflow', 'sinflow', 'boxflow'; 
eos         = 'linear';             % Exp.eos{ind};     % the equation of state between 'linear' and nonlinear: 'jmd95z' or 'mdjwf'; 
iceshelf    = 'tilted';             % Exp.iceshelf{ind};
strat       = Exp.strat{ind};       % 'const' for no stratification, 'linear' for linear stratification


% set parameters

acc = 'real*8';
prec='real*8';
ieee='b';

use grid
[Y,X]=meshgrid(y,x);




%% 1. Create bottom topography (topo)
% -----------------------------------------------

% Read the dimensions from table
w_ch    = Exp.w(ind);                               % width of channel center
alpha   = Exp.alpha(ind);                  % angle of channel sides 
d_ch    = H;                                        % depth in channel center
d_shelf = Exp.d_shelf(ind);                         % depth of shelf (flanks)
w       = round((d_ch-d_shelf)/alpha*2/500)*500;    % width of channel slopes
dh      = d_ch-d_shelf;                        

beta    = Exp.beta(ind);                            % slope of the whole channel



% calculate topography
switch topo
    case 'box'                      
        h = -ones(Nx,Ny)*H;
        h(1,:) = 0;
        h(:,1) = 0;
        
    case 'vshape'
        h = dh*ones(length(x),length(y));
        ih_min = length(y)/2-w/(2*dx);
        ih_max = length(y)/2+w/(2*dx);
        h(:,ih_min+1:ih_max) = abs(alpha*Y(:,ih_min+1:ih_max));
        h = h-H;

        % add slope of channel
        for ii = 1:length(x)
        h(ii,:) = h(ii,:) + beta*(max(x) - x(ii));
        end
 

    case 'sinshape'
        w_slope = dh/alpha;                          % width of the sloping walls

        h = dh*ones(length(x),length(y));
        for jj = find(y>=w/2-w_slope,1,'first'):find(y>=w/2,1,'first')
        h(:,jj) = dh/2*sin(pi/(2*w_slope/2)*(y(jj)-w_slope)) +dh/2;
        end
        for jj = find(y<=-(w/2-w_slope),1,'last')+1:find(y>=w/2-w_slope,1,'first')-1
        h(:,jj) = 0;
        end
        for jj = find(y<=-w/2,1,'last'):find(y<=-(w/2-w_slope),1,'last')
        h(:,jj) = dh/2*sin(pi/(2*w_slope/2)*(abs(y(jj))-w_slope)) +dh/2;
        end
        h = h-H;
      
        % add slope of channel
        for ii = 1:length(x)
        h(ii,:) = h(ii,:) + beta*(max(x) - x(ii));
        end
        
      case 'shelf'
        % add channel
        h= abs(alpha*(Y));
        %h(h>alpha*x(w/2/dy)+ beta*(max(x))- x(w/2/dy)) = alpha*x(w/2/dy)+ beta*(max(x)- x(w/2/dy));
        
        % add shelf slope
        h(1:w/2/dy,:) = alpha*X(1:w/2/dy,:);
        
        % add round corner      
        xc = 1:w/dx; 
        yc = 1:w/dy; 
        hc = cone(xc,yc,dh);
        
        ih_min = length(y)/2-w/(2*dx);
        ih_max = length(y)/2+w/(2*dx);
        h(1:w/2/dy,ih_min+1:Ny/2) = hc(1:w/2/dy,yc(end/2+1:end));
        h(1:w/2/dy,Ny/2+1:ih_max) = hc(1:w/2/dy,yc(1:end/2));
        
         % add slope of channel
        for ii = 1:length(x)
        h(ii,:) = h(ii,:) + beta*(max(x) -x(ii));
        end
        
        % lower channel to channel depth
        h = h-d_ch;
        
        % add flanks
        h(h>-d_shelf) = -d_shelf;
       
        iw_ch = (w_ch)/2/dy;
        h(:,1:Ny/2-iw_ch) = h(:,iw_ch+1:Ny/2);
        h(:,Ny/2+iw_ch+1:Ny) = h(:,Ny/2+1:Ny-iw_ch);
        h(:,Ny/2-iw_ch:Ny/2+iw_ch+1) = h(100,Ny/2);

end 

figure;
contourf(h,20);colorbar;
caxis([z(end) -250])

fid=fopen(['topog.channel'],'w',ieee); fwrite(fid,h,prec); fclose(fid);




%% 2. Create ice shelf topography (iceshelf)
% ----------------------------------------------
nxi         = find(x==Exp.xis(ind));                       % first grid of ice shelf
icefront    = Exp.icefront(ind);    % tilt of ice shelf front
dHdx        = Exp.dHdx(ind);

Hmax        = -dH;
Hmin        = -(dHdx*(x(end)-x(nxi))-Hmax);
xc          = x;


icetopo = min(Hmax,Hmin + dHdx*(max(x)-xc))'*ones(1,Ny);
icetopo(1:nxi-1,:) = 0;
  if icefront>0
      for i = 1:(icefront/dx)
          icetopo(nxi-icefront/dx+i,:) = i*Hmax/(icefront/dx);
      end
  end
icetopo(icetopo<-H+10) = -H+10;

figure;plot(icetopo(:,Ny/2))
fid=fopen(['icetopo'],'w','b'); fwrite(fid,icetopo,acc);fclose(fid);




%% 3. oceanographic setting
% ---------------------------------------------

% set grid config
dz = deltaZ*ones(1,Nz);
zgp1 = [0,cumsum(dz)];
zc = .5*(zgp1(1:end-1)+zgp1(2:end));
zg = zgp1(1:end-1);
dz = diff(zgp1);



% Set stratification
gravity  = 9.81;                        % Gravity
% talpha   = 2e-4;                     % in 1/oC; value equal to default in MITgcm
% sbeta    = 7.4e-04;                  % in 1/psu; value equal to default in MITgcm
t        = Exp.T(ind)*ones(Nz,1);
if strcmp(strat,'linear')
    s = Exp.S(ind):0.5/(Nz-1):Exp.S(ind)+0.5;
    s = s';
elseif strcmp(strat,'const')
    s = Exp.S(ind)*ones(Nz,1);
end

% parameter for linear eof
tref     = t(1);                      % use surface value for reference
sref     = s(1);                      % use surface value for reference
pref     = 0;
rhoConst = round(sw_dens(sref,tref,pref)*10)/10;      % calculate the density at tref and sref with seawater package

talpha = round(sw_alpha(sref, tref,pref)*10^7)*10^-7;
sbeta = round(sw_beta(sref,tref,pref)*10^6)*10^-6;
p = abs(zc)*gravity*rhoConst*1e-4;    % pressure in dbar

switch eos
    case 'linear'
      r = rhoConst*(1-talpha*(t-tref)+sbeta*(s-sref));      % drho is 0 for no stratification
     case 'jmd95z'
      r = densjmd95(s,t,p');                                    % drho is increasing without stratification, due to pressure dependency of rho
     case 'mdjwf'
      r = densmdjwf(s,t,p'); 
end



k=1;
dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
dp = p;
kp = 0;
while std(dp) > 1e-13
  phiHydF(k) = 0;
  p0 = p;
  kp = kp+1;
  for k = 1:Nz
    switch eos
     case 'linear'
      drho = rhoConst*(1-talpha*(t(k)-tref)+sbeta*(s(k)-sref))-rhoConst;      % drho is constant for no stratification
     case 'jmd95z'
      drho = densjmd95(s(k),t(k),p(k))-rhoConst;                                    % drho is increasing without stratification, due to pressure dependency of rho
     case 'mdjwf'
      drho = densmdjwf(s(k),t(k),p(k))-rhoConst;                                    % drho is increasing without stratification, due to pressure dependency of rho
     otherwise
      error(sprintf('unknown EOS: %s',eos))
    end
    phiHydC(k)   = phiHydF(k) + dzm(k)*gravity*drho/rhoConst;
    phiHydF(k+1) = phiHydC(k) + dzp(k)*gravity*drho/rhoConst;
  end
  switch eos
   case 'mdjwf'
    p = (gravity*rhoConst*abs(zc) + phiHydC*rhoConst)/gravity/rhoConst;
  end
  dp = p-p0;
end




%% pressure field under ice shelf
% --------------------------------------------------

% After modifying the code in calc_phi_hyd.F on Apr26,2012 this is the
% consistent way of computing phi0surf. For this, we need the grid
% information (hFacC's). For convenience, it's taken from a previous model
% run.
%
% The way of computing phi0surf consistent with code prior to Apr26,2012
% is recovered by setting drloc*dphi=0


hf=rdmds(['../run00',num2str(2),'/run/hFacC']);    % get the grid information from previous model run or...
% hf is a matrix that it zero in all grid points where there is topography
% or ice shelf and it is 1 whereever there is ocean.

% if there is no previous run from which hFacC can be extracted, run one
% time step without SHELFICEloadAnomalyFile to get the hFacC.
 


msk=sum(hf,3); msk(msk>0)=1;
phi0surf = zeros(Nx,Ny);

if ~isempty(hf)
for ix=1:Nx
  for iy=1:Ny
      k=max(find(abs(zg)<abs(icetopo(ix,iy))));
    % k=max(abs(zg)<abs(icetopo(ix,iy)));
    if isempty(k)
      k=0;
    end
    if k>0
      kp1=min(k+1,Nz);
      drloc=1-hf(ix,iy,k);
      %drloc=(abs(icetopo(ix,iy))-abs(zg(k)))/dz(k);
      dphi = phiHydF(kp1)-phiHydF(k);
      phi0surf(ix,iy) = (phiHydF(k)+drloc*dphi)*rhoConst*msk(ix,iy);   % pressure load ANOMALY in Pa
    end
  end
end
end
fid=fopen(['phi0surf'],'w','b'); fwrite(fid,phi0surf,acc);fclose(fid);


%% 4. Create velocity boundary conditions
% -------------------------------------------------


% velocity field at boundary;
u0 = Exp.uBCmax(ind);           % max velocity for the inflow
y0 = 5000;


switch flow
    case 'noflow'
        u0 = 0;
        
    case 'sinflow'
%         for i=1:Nz
%             umat(:,i) = u0*sin((y-y0)/(w/2-y0)*pi);
%         end
%         umat(abs(y-y0)>w/2-y0,:) = 0;
        umat = zeros(Ny,Nz);

        for i=1:Nz
            umat(1:Ny/2,i) = u0*sin((y(1:Ny/2)+y0)/(w/2-y0)*pi);
            umat(Ny/2+1:end,i) = u0*sin((y(Ny/2+1:end)-y0)/(w/2-y0)*pi);
        end
        umat(abs(y)<y0,:) = 0;        
        
        uWest = cat(3,umat,umat);
        uEast = cat(3,zeros(Ny,Nz),zeros(Ny,Nz)); %uEast = cat(3,(u0-du)*ones(Ny,Nz),(u0+du)*ones(Ny,Nz));
        vNorth = cat(3,zeros(Nx,Nz),zeros(Nx,Nz)); %uEast = cat(3,(u0-du)*ones(Ny,Nz),(u0+du)*ones(Ny,Nz));
        vSouth = cat(3,zeros(Nx,Nz),zeros(Nx,Nz)); %uEast = cat(3,(u0-du)*ones(Ny,Nz),(u0+du)*ones(Ny,Nz));
       
        
    case 'boxflow'
        umat = zeros(Ny,Nz);

        infl1 = 5000;          % at which y should inflow start (in m)
        infl2 = 15000;         % at which y should inflow stop (in m)
        ind_infl = Ny/2+infl1/dy+1:Ny/2+infl2/dy;
        umat(ind_infl,:) = u0*ones(length(ind_infl),Nz);            % positive inflow
        umat(Ny-ind_infl+1,:) = -u0*ones(length(ind_infl),Nz);         % return flow
        
        uWest = cat(3,umat,umat);
        uEast = cat(3,zeros(Ny,Nz),zeros(Ny,Nz)); %uEast = cat(3,(u0-du)*ones(Ny,Nz),(u0+du)*ones(Ny,Nz));
        vNorth = cat(3,zeros(Nx,Nz),zeros(Nx,Nz)); %uEast = cat(3,(u0-du)*ones(Ny,Nz),(u0+du)*ones(Ny,Nz));
        vSouth = cat(3,zeros(Nx,Nz),zeros(Nx,Nz)); %uEast = cat(3,(u0-du)*ones(Ny,Nz),(u0+du)*ones(Ny,Nz));

        
    case 'sideflow'
        vmat = zeros(Nx,Nz);
        for i=1:Nz
            vmat(abs(x)<w/2,i) = u0*sin((x(abs(x)<w/2))/(w/2)*pi);
        end
        vNorth = cat(3,-vmat,-vmat);
        vSouth = cat(3,-vmat,-vmat);
        
        uWest = cat(3,zeros(Ny,Nz),zeros(Ny,Nz)); %uEast = cat(3,(u0-du)*ones(Ny,Nz),(u0+du)*ones(Ny,Nz));
        uEast = cat(3,zeros(Ny,Nz),zeros(Ny,Nz)); %uEast = cat(3,(u0-du)*ones(Ny,Nz),(u0+du)*ones(Ny,Nz));

end


 
 % to test prescribing Eta in NonLin-FreeSurf formulation:
%  et1=0;                      % set Eta for the beginning
%  etWest = cat(2,+et1*ones(Ny,1),-et1*ones(Ny,1));
%  etEast = cat(2,-et1*ones(Ny,1),+et1*ones(Ny,1));

 
 

 fid=fopen('OB_WestU.bin','w',ieee); fwrite(fid,uWest ,prec); fclose(fid);
 fid=fopen('OB_EastU.bin','w',ieee); fwrite(fid,uEast ,prec); fclose(fid);
 fid=fopen('OB_NorthV.bin','w',ieee); fwrite(fid,vNorth ,prec); fclose(fid);
 fid=fopen('OB_SouthV.bin','w',ieee); fwrite(fid,vSouth ,prec); fclose(fid); 
 
%  fid=fopen('input/OBzonalS.bin','w',ieee); fwrite(fid,sZonal,prec); fclose(fid);
%  fid=fopen('input/OBzonalW.bin','w',ieee); fwrite(fid,wZonal,prec); fclose(fid);
%  fid=fopen('input/OB_WestH.bin','w',ieee); fwrite(fid,etWest,prec); fclose(fid);
%  fid=fopen('input/OB_EastH.bin','w',ieee); fwrite(fid,etEast,prec); fclose(fid);
%  



%% 5. Temperature field and salinity

% initial temperature field
T_ref   = Exp.T(ind)*ones(1,Nz);

Tin     = repmat(t,1,Nx,Ny);
T0      = permute(Tin,[2 3 1]);
fid     = fopen('T0.bin','w',ieee); fwrite(fid,T0,prec); fclose(fid);

% initial fields for salinity
Sin     =  repmat(s,1,Nx,Ny);
S0      = permute(Sin,[2 3 1]);
fid     = fopen('S0.bin','w',ieee); fwrite(fid,S0,prec); fclose(fid);

% prescribe temperature and salinity at Boundary
 
 tWest = cat(3,ones(Ny,1)*t',ones(Ny,1)*t');
 sWest = cat(3,ones(Ny,1)*s',ones(Ny,1)*s');
 tNorth = cat(3,ones(Nx,1)*t',ones(Nx,1)*t');
 sNorth = cat(3,ones(Nx,1)*s',ones(Nx,1)*s');

 fid=fopen('OBW_T.bin','w',ieee); fwrite(fid,tWest,prec); fclose(fid);
 fid=fopen('OBW_S.bin','w',ieee); fwrite(fid,sWest,prec); fclose(fid);
 fid=fopen('OBE_T.bin','w',ieee); fwrite(fid,tWest,prec); fclose(fid);
 fid=fopen('OBE_S.bin','w',ieee); fwrite(fid,sWest,prec); fclose(fid);
 fid=fopen('OBN_T.bin','w',ieee); fwrite(fid,tNorth,prec); fclose(fid);
 fid=fopen('OBN_S.bin','w',ieee); fwrite(fid,sNorth,prec); fclose(fid);
 fid=fopen('OBS_T.bin','w',ieee); fwrite(fid,tNorth,prec); fclose(fid);
 fid=fopen('OBS_S.bin','w',ieee); fwrite(fid,sNorth,prec); fclose(fid); 



%% plot input data
% -------------------------------------

% plot topography and ice shelf
icetopo2 = icetopo;
icetopo2(icetopo2==0) = NaN;
icetopo2(icetopo2<h) = NaN;
h1 = h;
%h1(h1>icetopo) = NaN;

f1 = myfigure(28,9);
plot3(X/1000,Y/1000,h1,'Color',[0.4 0.4 0.4]); hold on;
plot3(X/1000,Y/1000,icetopo2,'Color',[0.2 0 0.5]); box on; 
view([-24 16]); ylim([-Ny*dy/2/1000 Ny*dy/2/1000]); zlim([min(z) 0])
xlabel('y (km)'); ylabel('x (km)'); zlabel('depth (m)')
saveas(f1,'Topography_channel','jpeg'); 

% plot water depth
f2 = myfigure(28,10);
h2 = h - icetopo;
h2(h2>0) = 0;
contourf(X/1000,Y/1000,h2,40);cb = colorbar;
xlabel('x (km)'); ylabel('y (km)'); cb.Label.String = 'water depth [m]';
saveas(f2,'Waterdepth','jpeg'); 

% plot phi0surf   
f3 = myfigure(28,9);
contourf(X/1000,Y/1000,phi0surf); cb = colorbar;
xlabel('x (km)'); ylabel('y (km)'); cb.Label.String = 'phi0surf';
saveas(f3,'phi0surf','jpeg'); 
 
 
% plot velocity field at boundary
f4 = myfigure(28,9);
if strcmp(flow,'sideflow')
    contourf(x/1000,z,vmat',20);cb = colorbar;xlabel('x (km)');
    
else
contourf(y/1000,z,umat',20);cb = colorbar;xlabel('y (km)');
end
ylabel('depth (m)'); cb.Label.String = 'velocity (m/s)';
saveas(f4,'uWest','jpeg'); 

% plot stratification
long = -128; lat = -74;
p = gsw_p_from_z(z,lat);
SA = gsw_SA_from_SP(s',p,long,lat);
CT = gsw_CT_from_t(SA,t',p);
Rho = gsw_rho(SA,CT,p);

f5 = myfigure(28,18);
subplot(231)
plot(t,z); xlabel('T (^o C)'); ylabel('z (m)'); xlim([-0.5 0.5])
subplot(232)
plot(s,z); xlabel('s (psu)'); ylabel('z (m)'); xlim([34.3 35])
subplot(233)
plot(r,z); xlabel('\rho (kg/m^3)'); ylabel('z (m)'); xlim([1027.5 1031])
text(r(10),z(4/5*Nz),{['\rho_0 = ',num2str(rhoConst)],['\alpha = ',num2str(talpha)],['\beta =',num2str(sbeta)],['T0 = ',num2str(tref)],['S0 =',num2str(sref)]})
title(['eos = ',eos])

subplot(234)
plot(CT,p); xlabel('CT (^o C)'); ylabel('pressure'); set(gca,'YDir','reverse'); xlim([-0.5 0.5])
subplot(235)
plot(SA,p); xlabel('SA (g/kg)'); ylabel('pressure'); set(gca,'YDir','reverse'); xlim([34.3 35])
subplot(236)
plot(Rho,p); xlabel('\rho (kg/m^3)'); ylabel('pressure'); set(gca,'YDir','reverse'); xlim([1027.5 1031])



saveas(f5,'Stratification','jpeg'); 
 

 
