% Create constants.mat that lists all constants used in the M-AMOC1.0 AMOC-emulator code

%% GCM
modelname='UVic';

%% Time constants
eq_length =500; 			% length of equilibrium simulation in years
trans_length_hist =156;  		% length of historical simulation in years
trans_length_rcp45nomelt =295; 		% length of RCP4.5 simulation in years
trans_length_rcp45gGISmelt = 295;  	% length of RCP4.5 GISmelt simulation in years
trans_length_rcp85nomelt = 295; 	% length of RCP8.5 simulation in years
trans_length_rcp85gGISmelt = 295;  	% length of RCP8.5 GISmelt simulation in years

ebm_res = 52; 				% EBM time-resolution (timesteps in every year)
amoc_res = 13; 				% AMOC time-resolution (timesteps in every year) % Note that ebm_res/amoc_res has to be an integer
daysinyr = 364; 			% Number of days in a year

eq_length = eq_length(model)*ebm_res;	% length of equilibrium simulation in EBM time-steps
trans_length = (trans_length_hist+trans_length_rcp45nomelt+trans_length_rcp45gGISmelt+trans_length_rcp85nomelt+trans_length_rcp85gGISmelt)*ebm_res; % Totoal length of transient simulations in EBM time-steps
dt_ebm = daysinyr/ebm_res*86400.0; 	% EBM time step length in seconds
dt_amoc = 1/amoc_res; 			% AMOC time step length in years
yr2sec=60*60*24*daysinyr; 		% number of seconds in year

%% Fitting procedure constants
numraninit=100; 			% number of iterations with different random initial parameter values for fitting procedure
maxlength=300;				% maximum number of simulations in fitting procedure before iteration is cancelled.
numpara=7; 				% number of parameters tuned in fitting procedure
parinci=2; 				% perturbation applied to initial Zickfeld et al. values to find initial parameter values (2=+-200%)
raninitpar=1; 				% Set the calculations of the initial parameter values in fitting procedure to random (1) or prescribed (0)
Tstati=0.2; 				% Initial value for statistical 'Temperature'
uphti=0.6; 				% Initial value for uphill climbing threshold
R2gradStopcrit=-0.05; 			% Stopping criterion for fitting procedure. -0.05 means that at least a 5% decrease in R2 has to be made over the last 'R2gradlength' experiments, otherwise iteration is canceled.
R2gradlength=50; 			% Number of experiments over which the R2 gradient is calculated for the stopping criterion.
R2exept=3e4; 				% Maximum R2-value of a parameter set to be exepted as a successfull fit.
numbestmod=10; 				% Number of best emulator fits used as AMOC-emulators

%% Geometry constants
deg2rad  = pi/180.0;        		% degrees to radians conversion factor
deg2dist = 111194.9;        		% distance that corresponds to one degree of latitude (m)

% EBM
jmt = 7;                    		% = number of grid points in the latitudinal (north-south) direction (calculated points are from 2 through jmt-1, end points are boundaries)

% Initialize arrays
dytdeg=zeros(1,jmt);  
dyt=zeros(1,jmt);  
yt=zeros(1,jmt);
sint=zeros(1,jmt); 
cst=zeros(1,jmt);
dyu=zeros(1,jmt);
dyudeg=zeros(1,jmt);
yu=zeros(1,jmt);
csu=zeros(1,jmt);

% Define temperature grid cells
dxt = 1; 				% T grid cell zonal width at equator (m)
yt = [0 -80 -50 8 58 80 0]; 		% midpoints T grid cell (deg N)
dytdeg = [0 20 40 75 25 20 0]; 		% T grid cell meridional height (deg)
dyt = dytdeg*deg2dist;     		% T grid cell meridional height (m)
sint = sin(yt*deg2rad); 		% sine of T grid cell latitude
cst  = cos(yt*deg2rad); 		% cosine of T grid cell latitude

% Define velocity grid cells
dxu         = 1;    			% U grid cell zonal width at equator (m)
yu = [0 -90 -70 -30 45 70 90 0]; 	% midpoints U grid cell (deg N)
dyudeg(2:6)   = dytdeg(1:5)/2+dytdeg(2:6)/2; % U grid cell meridional height (deg)
dyu = dyudeg*deg2dist;     		% U grid cell meridional height (m)
csu(2:jmt-1)  = cos(yu(2:jmt-1)*deg2rad); % cosine of U grid cell latitude

% Define height of atmosphere
hatm = 8400.0;              		% height of atmosphere (m)

% AMOC
V=[1.1*10^(17) 0.68*10^(17) 0.4*10^(17) 0.05*10^(17)]; % Volumes of AMOC basins (m3)
z=[3000 1000 3000]; 			% Depth of AMOC basins (m)

%% EBM constants
% Physical constants
rhoair  = 1.005;           		% density of air (kg m^(-3))
cpair   = 1000.0;          		% specific heat of air (J kg^(-1) K^(-1))
tzero   = 273.15;          		% freezing point (K)
stebol  = 5.67E-08;        		% Stefan Boltzmann's constant (W m^(-2) K ^(-4))
CtoK    = 272.15;          		% Conversion from Celsius to Kelvin

% Model parameters
apln(1:jmt)      = 0.0;                 % planetary albedo (-)
apln(2:6)        = 0.6-0.4*cst(2:6);    % Parametrization of planetary albedo From lecture notes Thomas Stocker (2008)
diff_cnt(1:jmt)  = 0.0;                 % Diffusion coefficient for northern face of T grid cell (latitude-dependent, m^2 s^(-1)
diff_cnt(2:6)    = 2*10^5*cos(deg2rad*yu(3:7)); % From lecture notes Andreas Schmittner
ceff             = rhoair*cpair*hatm;   % Atmospheric heat capacity (J m^(-2) K^(-1))
solin(1:jmt)     = 0.0;                 % solar insolation (W m^(-2))
for j=2:jmt-1
    solin(j)=0;
    for lat=yu(j):yu(j+1)
        solin(j) = solin(j)+295+125*cos(2*lat*deg2rad); % solar insolation on T grid
    end
    solin(j)=solin(j)/(yu(j+1)-yu(j));
end

%% AMOC constants
rhowat  = 1025;            		% density of sea water (kg m^(-3)) 
cpwat   = 4000.0;          		% specific heat of water (J kg^(-1) K^(-1))
S0      = 35;              		% Global reference salinity (psu)

Teq = [7 11.4 4.7 4.7]+CtoK; 		% Initial ocean temperature (K) from Zickfeld et al.
Seq= [35 35.08 35.02 35.02]; 		% Initial ocean salinity (psu) from Zickfeld et al. + average of 35psu

ki=2.54*10^18; 				% Flow constant (m^3 yr^(-1)) from Zickfeld et al.
thaui=7.3*10^8;				% thermal coupling constant (J yr^(-1) m^(-2) K^(-1))  from Zickfeld et al.

F1i=0.014;				% Initial freshwater flux F1 (Sv) from Zickfeld et al.
F2i=0.065;				% Initial freshwater flux F2 (Sv) from Zickfeld et al.

h1i = -0.0047; 				% Initial Freshwater flux F1 to global temperature change relation constant (Sv/K) from from Zickfeld et al.
h2i = 0.0139; 				% Initial Freshwater flux F2 to global temperature change relation constant (Sv/K) from from Zickfeld et al.

lambdai=thaui./(cpwat*rhowat*z); 	% lambda values (yr^(-1))

save('constants.mat')
