Index: /issm/trunk/src/m/classes/public/modeldefault/defaultparams.m
===================================================================
--- /issm/trunk/src/m/classes/public/modeldefault/defaultparams.m	(revision 2216)
+++ /issm/trunk/src/m/classes/public/modeldefault/defaultparams.m	(revision 2216)
@@ -0,0 +1,113 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Some hardcoded parameters for this deck%%%%%%%%%%
+
+%some parameterization for this parameter file :)
+thicknesspath=options.thicknesspath;
+firnpath=options.firnpath;
+surfacepath=options.surfacepath;
+mosaicpath=options.mosaicpath;
+temperaturepath=options.temperaturepath;
+heatfluxpath=options.heatfluxpath;
+
+%Solution parameters
+
+	%parallelization 
+	md.cluster=oshostname();
+	md.np=8;
+	md.time=60;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	disp('      reading bedmap thicknesses');
+	md.thickness=InterpFromFile(md.x,md.y,thicknesspath,10);
+
+	disp('      reading firn layer');
+	md.firn_layer=InterpFromFile(md.x,md.y,firnpath,0);
+
+	disp('      reading Bamber dem');
+	md.surface=InterpFromFile(md.x,md.y,surfacepath,10);
+
+	%correct the surface by taking into account the firn layer
+	rho_ice=917;
+	rho_firn=830;
+	firn_layer_correction=md.firn_layer*(1-rho_firn/rho_ice);
+	md.surface=md.surface-firn_layer_correction;
+	
+	%Correct thickness by assuming hydrostatic equilibrium 10km down the grounding line
+	md=ThicknessCorrection(md);
+
+	%some corrections
+	minsurf=1;
+	pos=find(isnan(md.surface) | (md.surface<=0));
+	md.surface(pos)=minsurf;
+	pos=find(isnan(md.thickness) | (md.thickness<=0));
+	md.thickness(pos)=minsurf/(1-md.rho_ice/md.rho_water);
+	md.bed=md.surface-md.thickness;
+
+	disp('      reading velocities from Rignot');
+	md=plugvelocities(md,mosaicpath,0);
+	
+	%drag md.drag or stress
+	md.drag_type=2; %0 none 1 plastic 2 viscous
+	md.drag=300*ones(md.numberofgrids,1); %q=1.
+		
+	%zones of high md.drag
+	%[rhighmd.drag]=ArgusContourToMesh(md.elements,md.x,md.y,expread('HighDrag.exp',1),'node');
+	%pos=find(highmd.drag);md.drag(pos)=10^3;
+
+	%Take care of iceshelves: no drag md.drag
+	pos=find(md.elementoniceshelf);
+	md.drag(md.elements(pos,:))=0;
+	md.p=ones(md.numberofelements,1);
+	md.q=ones(md.numberofelements,1);
+
+	%Load md.temperature from Giovinetto:
+	disp('      loading temperature');
+	md.temperature=InterpFromFile(md.x,md.y,temperaturepath,253);
+	while ~isempty(find(isnan(md.temperature))),
+		pos=find(isnan(md.temperature));
+		if ((pos+1)<=length(md.temperature)),
+			md.temperature(pos)=md.temperature(pos+1);
+		else
+			md.temperature(pos)=md.temperature(pos-1);
+		end
+	end
+
+	%flow law
+	disp('      creating flow law paramters');
+	md.n=3*ones(md.numberofelements,1);
+	md.B=paterson(md.temperature);
+
+	%zones of shear margin softening
+	%[rweakb]=ArgusContourToMesh(md.elements,md.x,md.y,expread('Weakmd.BPIG.exp',1),'node');
+	%pos=find(weakb);md.B(pos)=.3*md.B(pos);
+		
+	%rifts: none for now.
+	%if isstruct(rifts), %we have rifts, choose whether they are full of water or air.
+	%	for i=1:length(rifts), 
+	%		rifts(i).fill='water'; %choice:'air','water','ice'
+	%		rifts(i).friction=10^11;
+	%	end
+	%end
+
+	disp('      creating accumulation rates');
+	md.accumulation=ones(md.numberofgrids,1); %1m/a
+	
+	%Deal with boundary conditions:
+
+	disp('      thermal model');
+	md.melting=zeros(md.numberofgrids,1);
+	md.observed_temperature=md.temperature;
+	
+	disp('      reading geothermal flux');
+	load(heatfluxpath);
+	md.geothermalflux=InterpFromGrid(x_m,y_m,heatflux_Antarctica,md.x,md.y,80);
+	pos=find(md.geothermalflux==0);md.geothermalflux(pos)=80;
+	md.geothermalflux=md.geothermalflux/1000; %map is given in mW/m^2, we need it in W/m^2
+
+	disp('      boundary conditions');
+	md=SetMarineIceSheetBC(md);
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+md.counter=3;
Index: /issm/trunk/src/m/classes/public/modeldefault/modeldefault.m
===================================================================
--- /issm/trunk/src/m/classes/public/modeldefault/modeldefault.m	(revision 2216)
+++ /issm/trunk/src/m/classes/public/modeldefault/modeldefault.m	(revision 2216)
@@ -0,0 +1,66 @@
+function md=modeldefault(varargin),
+%MODELDEFAULT - build a model with default parameters
+%
+%   This routine is used to create and parameterize a model.
+%   Several options are available, but only the domain outline
+%   is mandatory.
+%   Default options will be applied when not specified.
+%
+%   Available options:
+%      o 'domainoutline' followed by the name of an ARGUS file
+%        describing the domain to be meshed
+%      o 'resolution': mesh resolution
+%      o 'iceshelves': ice shelves contours in an ARGUS file
+%      o 'islands': islands contours in an ARGUS file
+%      o 'thicknesspath': path to the thickness data file
+%      o 'firnpath': path to the thickness data file
+%      o 'surfacepath': path to the thickness data file
+%      o 'velocitypath': path to the thickness data file
+%      o 'temperaturepath': path to the thickness data file
+%      o 'heatfluxpath': path to the thickness data file
+%      o 'numlayers': number of layers in the mesh vertical extrusion
+%      o 'extrusionexponent': exponent used in the vertical extrusion
+%      o 'elementstype': type of elements for diagnostic runs
+%
+%   Usage:
+%      md=modeldefault(varargin)
+%
+%   Examples:
+%      md=modeldefault('domainoutline','Pig.exp');
+
+%some checks
+if nargout~=1
+	help modeldefault
+	error('modeldefault error message: bad usage')
+end
+
+%recover options
+disp(sprintf('\nRecovering options...\n'));
+options=optionlist2cell(varargin{:});
+
+%process options
+options=modeldefault_process(options);
+
+%initialize output:
+md=model;
+
+%Mesh the model
+disp(sprintf('\nMesh generation...'));
+md=mesh(md,options.domainoutline,options.resolution);
+
+%geography
+disp(sprintf('Geography...'));
+md=geography(md,options.iceshelves,options.islands);
+
+%parameterization
+disp(sprintf('Model parameterization...'));
+defaultparams;
+
+%extrusion
+disp(sprintf('Mesh extrusion...'));
+md=extrude(md,options.numlayers,options.extrusionexponent);
+
+%elements type
+disp(sprintf('Setting elemtens type...'));
+md=setelementstype(md,options.elementstype,'all');
+disp(sprintf('done'));
Index: /issm/trunk/src/m/classes/public/modeldefault/modeldefault_process.m
===================================================================
--- /issm/trunk/src/m/classes/public/modeldefault/modeldefault_process.m	(revision 2216)
+++ /issm/trunk/src/m/classes/public/modeldefault/modeldefault_process.m	(revision 2216)
@@ -0,0 +1,148 @@
+function optionstruct=modeldefault_process(options),
+%MODELDEFAULT_PROCESS - process modeldefault options
+%
+%   Usage:
+%      optionstruct=modeldefault_process(options)
+
+%retrieve some paths
+global ISSM_DIR
+modeldatapath='/u/wilkes-r1b/larour/ModelData';
+
+%initialize output
+optionstruct=struct();
+
+%domainoutline
+optionstruct.domainoutline=find_option(options,'domainoutline');
+if isempty(optionstruct.domainoutline),
+	error('modeldefault_process error pessage: at least the option ''domainoutline'' must be specified')
+end
+
+%resolution
+optionstruct.resolution=find_option(options,'resolution');
+if isempty(optionstruct.resolution),
+	optionstruct.resolution=10000;
+end
+disp(sprintf('   %-18s: %g','resolution',optionstruct.resolution));
+
+%iceshelves
+optionstruct.iceshelves=find_option(options,'iceshelves');
+if isempty(optionstruct.iceshelves),
+	iceshelves=[ISSM_DIR '/../models/Antarctica/Exp_Par/Iceshelves.exp'];
+	if exist(iceshelves),
+		optionstruct.iceshelves=iceshelves;
+	else
+		optionstruct.iceshelves='';
+	end
+end
+disp(sprintf('   %-18s: %s','iceshelves',optionstruct.iceshelves));
+
+%islands
+optionstruct.islands=find_option(options,'islands');
+if isempty(optionstruct.islands),
+	islands=[ISSM_DIR '/../models/Antarctica/Exp_Par/Islands.exp'];
+	if exist(islands),
+		optionstruct.islands=islands;
+	else
+		optionstruct.islands='';
+	end
+end
+disp(sprintf('   %-18s: %s','islands',optionstruct.islands));
+
+%thicknesspath
+optionstruct.thicknesspath=find_option(options,'thicknesspath');
+if isempty(optionstruct.thicknesspath),
+	thicknesspath=[modeldatapath '/BedMap/gridded/thickness.mat'];
+	if exist(thicknesspath),
+		optionstruct.thicknesspath=thicknesspath;
+	else
+		optionstruct.thicknesspath='';
+	end
+end
+disp(sprintf('   %-18s: %s','thicknesspath',optionstruct.thicknesspath));
+
+%firnpath
+optionstruct.firnpath=find_option(options,'firnpath');
+if isempty(optionstruct.firnpath),
+	firnpath=[modeldatapath '/BroekeFirn1km/firn.mat'];
+	if exist(firnpath),
+		optionstruct.firnpath=firnpath;
+	else
+		optionstruct.firnpath='';
+	end
+end
+disp(sprintf('   %-18s: %s','firnpath',optionstruct.firnpath));
+
+%surfacepath
+optionstruct.surfacepath=find_option(options,'surfacepath');
+if isempty(optionstruct.surfacepath),
+	surfacepath=[modeldatapath '/BamberDEMAntarctica1km/surface_smooth30_lowslope.mat'];
+	if exist(surfacepath),
+		optionstruct.surfacepath=surfacepath;
+	else
+		optionstruct.surfacepath='';
+	end
+end
+disp(sprintf('   %-18s: %s','surfacepath',optionstruct.surfacepath));
+
+%mosaicpath
+optionstruct.mosaicpath=find_option(options,'mosaicpath');
+if isempty(optionstruct.mosaicpath),
+	mosaicpath=[modeldatapath '/RignotAntarcticaVelMosaicRampErsAlos/RignotAntVel.mat'];
+	if exist(mosaicpath),
+		optionstruct.mosaicpath=mosaicpath;
+	else
+		optionstruct.mosaicpath='';
+	end
+end
+disp(sprintf('   %-18s: %s','mosaicpath',optionstruct.mosaicpath));
+
+%temperaturepath
+optionstruct.temperaturepath=find_option(options,'temperaturepath');
+if isempty(optionstruct.temperaturepath),
+	temperaturepath=[modeldatapath '/GiovinettoZwallyTemperatures92/Giovinetto_Temperatures.mat'];
+	if exist(temperaturepath),
+		optionstruct.temperaturepath=temperaturepath;
+	else
+		optionstruct.temperaturepath='';
+	end
+end
+disp(sprintf('   %-18s: %s','temperaturepath',optionstruct.temperaturepath));
+
+%heatfluxpath
+optionstruct.heatfluxpath=find_option(options,'heatfluxpath');
+if isempty(optionstruct.heatfluxpath),
+	heatfluxpath=[modeldatapath '/HeatfluxAntarctica/RignotHeatFlux.mat'];
+	if exist(heatfluxpath),
+		optionstruct.heatfluxpath=heatfluxpath;
+	else
+		optionstruct.heatfluxpath='';
+	end
+end
+disp(sprintf('   %-18s: %s','heatfluxpath',optionstruct.heatfluxpath));
+
+%numlayers
+optionstruct.numlayers=find_option(options,'numlayers');
+if isempty(optionstruct.numlayers),
+	optionstruct.numlayers=0;
+end
+disp(sprintf('   %-18s: %g','numlayers',optionstruct.numlayers));
+
+%extrusionexponent
+optionstruct.extrusionexponent=find_option(options,'extrusionexponent');
+if isempty(optionstruct.extrusionexponent),
+	optionstruct.extrusionexponent=3;
+end
+disp(sprintf('   %-18s: %g','extrusionexponent',optionstruct.extrusionexponent));
+
+
+%elementstype
+optionstruct.elementstype=find_option(options,'elementstype');
+if isempty(optionstruct.elementstype),
+	if optionstruct.numlayers>0,
+		optionstruct.elementstype='pattyn';
+	else
+		optionstruct.elementstype='macayeal';
+	end
+end
+disp(sprintf('   %-18s: %s','elementstype',optionstruct.elementstype));
+