Index: /issm/trunk-jpl/src/m/parameterization/interpISMIP6AntarcticaOcn.m
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/interpISMIP6AntarcticaOcn.m	(revision 28040)
+++ /issm/trunk-jpl/src/m/parameterization/interpISMIP6AntarcticaOcn.m	(revision 28040)
@@ -0,0 +1,83 @@
+function basalforcings = interpISMIP6AntarcticaOcn(md,model_name,scenario)
+%interpISMIP6AntarcticaOcn - interpolate chosen ISMIP6 atmospheric forcing to model
+%
+%   Input:
+%     - md (model object)
+%     - model_name  (string): name of the climate model 
+%                             Examples: cnrm-esm2-1 ccsm4cesm2 cnrm-cm6-1 csiro-mk3-6-0
+%                                       hadgem2-es ipsl-cm5a-mr miroc-esm-chem noresm1-m ukesm1-0-ll
+%
+%     - scenario    (string): name of the climate scenario
+%                             Examples: rcp2.6, rcp8.5, ssp126, ssp585
+%   Output:
+%     - basalforcings: prepared to be input directly into md.smb
+%
+%   Examples:
+%      md.basalforcings = interpISMIP6AntarcticaOcn(md,'miroc-esm-chem','rcp8.5');
+
+% Find appropriate directory
+switch oshostname(),
+	case {'totten'}
+		path='/totten_1/ModelData/ISMIP6/Projections/AIS/Ocean_Forcing/';
+	otherwise
+		error('machine not supported yet, please provide your own path');
+end
+
+rootname = [path model_name '_' scenario '/'];
+if ~exist(rootname,'dir')
+   error(['this path does not exist or the ' model_name ' and ' scenario ' are not available in this combination.']);
+end
+
+%load TF data
+tfnc    = [rootname '/1995-2100/' upper(model_name) '_thermal_forcing_8km_x_60m.nc'];
+x_n     = double(ncread(tfnc,'x'));
+y_n     = double(ncread(tfnc,'y'));
+tf_data = double(ncread(tfnc,'thermal_forcing'));
+z_data  = double(ncread(tfnc,'z'));
+
+%Build tf cell array
+t = 1:size(tf_data,4);
+tf = cell(1,1,size(tf_data,3));
+for i=1:size(tf_data,3)  %Iterate over depths
+	disp(['   == Interpolating over depth ' num2str(i) '/' num2str(size(tf_data,3))]);
+	
+	temp_matrix=[];
+	for ii=1:size(tf_data,4) %Iterate over time steps
+		%temp_tfdata=InterpFromGridToMesh(x_n,y_n,tf_data(:,:,i,ii)',md.mesh.x,md.mesh.y,0);
+		temp_tfdata=InterpFromGrid(x_n,y_n,tf_data(:,:,i,ii)',md.mesh.x,md.mesh.y);
+		temp_matrix = [temp_matrix temp_tfdata];
+	end
+	tf{:,:,i} = [temp_matrix ; t];
+end
+
+%load Delta and gamma data
+deltatnc_median = [path '/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_median.nc'];
+basin_datanc    = [path '/imbie2/imbie2_basin_numbers_8km.nc'];
+deltaT_median   = double(ncread(deltatnc_median,'deltaT_basin'));
+gamma0_median   = double(ncread(deltatnc_median,'gamma0'));
+basinid_data    = double(ncread(basin_datanc,'basinNumber'));
+
+disp('   == Interpolating basin Id');
+num_basins = length(unique(basinid_data));
+deltat_median = NaN(length(unique(basinid_data)),1);
+
+for i=0:num_basins-1
+	pos = find(basinid_data==i);
+	deltat_temp = deltaT_median(pos);
+	deltat_temp = deltat_temp(1);
+	deltat_median(i+1) = deltat_temp;
+end
+
+%Deal with basins ID
+x_el = mean(md.mesh.x(md.mesh.elements),2);
+y_el = mean(md.mesh.y(md.mesh.elements),2);
+basinid = InterpFromGrid(x_n,y_n,basinid_data',x_el, y_el, 'nearest')+1;
+
+%Set ISMIP6 basal melt rate parameters
+basalforcings            = basalforcingsismip6(md.basalforcings);
+basalforcings.basin_id   = basinid;
+basalforcings.num_basins = num_basins;
+basalforcings.delta_t    = deltat_median;
+basalforcings.tf_depths  = z_data';
+basalforcings.gamma_0    = gamma0_median;
+basalforcings.tf         = tf;
Index: /issm/trunk-jpl/src/m/parameterization/interpISMIP6AntarcticaSMB.m
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/interpISMIP6AntarcticaSMB.m	(revision 28040)
+++ /issm/trunk-jpl/src/m/parameterization/interpISMIP6AntarcticaSMB.m	(revision 28040)
@@ -0,0 +1,72 @@
+function smb = interpISMIP6AntarcticaSMB(md,model_name,scenario)
+%interpISMIP6AntarcticaSMB - interpolate chosen ISMIP6 atmospheric forcing to model
+%
+%   Input:
+%     - md (model object)
+%     - model_name  (string): name of the climate model 
+%                             Examples: CESM2 CNRM_CM6 CNRM_ESM2 CSIRO-Mk3-6-0 HadGEM2-ES IPSL-CM5A-MR 
+%                                       ccsm4 miroc-esm-chem noresm1-m
+%     - scenario    (string): name of the climate scenario
+%                             Examples: rcp26, rcp85, ssp126, ssp585
+%   Output:
+%     - smb: prepared to be input directly into md.smb
+%
+%   Examples:
+%      md.smb = interpISMIP6AntarcticaSMB(md,'miroc-esm-chem','rcp85');
+
+% Find appropriate directory
+switch oshostname(),
+	case {'totten'}
+		path='/totten_1/ModelData/ISMIP6/Projections/AIS/Atmosphere_Forcing/';
+	otherwise
+		error('machine not supported yet, please provide your own path');
+end
+
+rootname = [path model_name '_' scenario '/'];
+if ~exist(rootname,'dir')
+   error(['this path does not exist or the ' model_name ' and ' scenario ' are not available in this combination.']);
+end
+
+disp('   == loading TS and SMB climatology data');
+smbclimnc           = [rootname '/Regridded_2km/' upper(model_name) '_2km_clim_1995-2014.nc'];
+lat                 = double(ncread(smbclimnc,'lat'));
+lon                 = double(ncread(smbclimnc,'lon'));
+smb_clim_data       = double(ncread(smbclimnc,'smb_clim'));
+ts_clim_data        = double(ncread(smbclimnc,'ts_clim'));
+
+disp('   == loading TS and SMB anomoly data');
+smbanomnc           = [rootname '/Regridded_2km/' upper(model_name) '_2km_anomaly_1995-2100.nc'];
+smb_anomaly_data    = double(ncread(smbanomnc,'smb_anomaly'));
+ts_anomaly_data     = double(ncread(smbanomnc,'ts_anomaly'));
+
+%Create SMB and TS matrix
+disp('   == Interpolating on model');
+t=[1:size(smb_anomaly_data,3)];
+[x_n y_n]=ll2xy(lat(:,1),lon(:,1),-1);
+y_n = x_n;
+smb_clim = InterpFromGridToMesh(x_n,y_n,smb_clim_data',md.mesh.x,md.mesh.y,0);
+ts_clim  = InterpFromGridToMesh(x_n,y_n,ts_clim_data',md.mesh.x,md.mesh.y,0);
+temp_matrix_smb = []; temp_matrix_ts = [];
+for i = 1:size(smb_anomaly_data,3)
+	%SMB
+	temp_smb        = InterpFromGridToMesh(x_n,y_n,smb_anomaly_data(:,:,i)',md.mesh.x,md.mesh.y,0);
+	temp_smb        = temp_smb+smb_clim;
+	temp_matrix_smb = [temp_matrix_smb temp_smb];
+	%TS
+	temp_ts         = InterpFromGridToMesh(x_n,y_n,ts_anomaly_data(:,:,i)',md.mesh.x,md.mesh.y,0);
+	temp_ts         = temp_smb+smb_clim;
+	temp_matrix_ts  = [temp_matrix_ts temp_ts];
+	clear temp_smb; clear temp_ts;
+end
+
+% convert to m/yr
+rhoi = md.materials.rho_ice;
+temp_matrix_smb = temp_matrix_smb*(31556926/1000)*(1000/rhoi);
+
+%Save Data (1995-2100)
+smb = SMBforcing();
+smb.mass_balance = [temp_matrix_smb ; t];
+
+%What do we do with surface temp?
+%md.miscellaneous.dummy.ts = [temp_matrix_ts ; t];
+end