Index: /issm/trunk/test/NightlyRun/test1301.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1301.m	(revision 5046)
+++ /issm/trunk/test/NightlyRun/test1301.m	(revision 5046)
@@ -0,0 +1,35 @@
+% This file can be run to check that the melting in simple conduction is correctly modeled.
+% There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature
+% at upper surface and an impose flux at its base. The result must be a linear temperature from the upper to the lower
+% surface with an imposed slope (Geothermal flux). if it is not the case, something is thermal modeling has been changed...
+% Just run this file in Matlab, with a properly setup Ice code. 
+
+md=model;
+md=mesh(md,'../Exp/Square.exp',100000);
+md=geography(md,'','');
+md=parameterize(md,'../Par/SquareThermal.par');
+md=extrude(md,3,2);
+md.verbose=2;
+md=setelementstype(md,'Pattyn','all');
+
+%Some conditions specific to melting test
+md.pressure=zeros(md.numberofgrids,1);
+md.observed_temperature=273.15*ones(md.numberofgrids,1);
+md.temperature=md.observed_temperature;
+md.spctemperature(:,2)=md.observed_temperature;
+md.rheology_B=paterson(md.temperature);
+
+%analytical results
+%melting heat = geothermal flux
+%Mb*L*rho=G   => Mb=G/L*rho
+melting=md.geothermalflux/(md.rho_ice*md.latentheat)*md.yts;
+
+%modeled  results
+md=solve(md,'analysis_type',ThermalSolutionEnum);
+
+%plot results
+comp_melting=zeros(md.numberofgrids,1);
+comp_melting(md.results.ThermalSolution.MeltingRate.index)=md.results.ThermalSolution.MeltingRate.value;
+relative=abs((comp_melting-melting)./melting)*100;
+relative(find(comp_melting==melting))=0;
+plotmodel(md,'data',comp_melting,'title','modeled melting','data',melting,'title','analytical melting','data',comp_melting-melting,'title','absolute error','data',relative,'title','relative error [%]','layer#all',1,'caxis#2',[1.02964 1.02966]*10^-4)
Index: /issm/trunk/test/NightlyRun/test1301_nightly.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1301_nightly.m	(revision 5046)
+++ /issm/trunk/test/NightlyRun/test1301_nightly.m	(revision 5046)
@@ -0,0 +1,3 @@
+field_names     ={'Melting'};
+field_tolerances={1e-13};
+field_values    ={comp_melting};
Index: /issm/trunk/test/NightlyRun/test1302.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1302.m	(revision 5046)
+++ /issm/trunk/test/NightlyRun/test1302.m	(revision 5046)
@@ -0,0 +1,35 @@
+% This file can be run to check that the advection-diffusion  is correctly modeled.
+% There is u=v=0 and w=cst everywhere the only thermal boundary conditions are an imposed temperature
+% at upper surface and an impose flux at its base.
+% Just run this file in Matlab, with a properly setup Ice code. 
+
+md=model;
+md=mesh(md,'../Exp/Square.exp',100000);
+md=geography(md,'','');
+md=parameterize(md,'../Par/SquareThermal.par');
+md=extrude(md,30,1);   %NB: the more one extrudes, the better (10-> relative~0.35%, 20->0.1%, 30->0.05%)
+md=setelementstype(md,'Pattyn','all'); 
+
+%Thermal boundary conditions
+pos1=find(md.elementonbed);     md.spctemperature(md.elements(pos1,1:3),1)=1; md.spctemperature(md.elements(pos1,1:3),2)=10;
+pos2=find(md.elementonsurface); md.spctemperature(md.elements(pos2,4:6),1)=1; md.spctemperature(md.elements(pos2,4:6),2)=0;
+md.vz=0.1*ones(md.numberofgrids,1);
+md.vel=sqrt( md.vx.^2+ md.vy.^2+ md.vz.^2);
+md.pressure=zeros(md.numberofgrids,1);
+
+%analytical results
+%d2T/dz2-w*rho_ice*c/k*dT/dz=0   T(surface)=0  T(bed)=10   => T=A exp(alpha z)+B
+alpha=0.1/md.yts*md.rho_ice*md.heatcapacity/md.thermalconductivity;   %alpha=w rho_ice c /k  and w=0.1m/an
+A=10/(exp(alpha*(-1000))-1);    %A=T(bed)/(exp(alpha*bed)-1)  with bed=-1000 T(bed)=10
+B=-A;
+md.observed_temperature=A*exp(alpha*md.z)+B;
+
+%modeled  results
+md=solve(md,'analysis_type',ThermalSolutionEnum);
+
+%plot results
+comp_temp=zeros(md.numberofgrids,1);
+comp_temp(md.results.ThermalSolution.Temperature.index)=md.results.ThermalSolution.Temperature.value;
+relative=abs((comp_temp-md.observed_temperature)./md.observed_temperature)*100;
+relative(find(comp_temp==md.observed_temperature))=0;
+plotmodel(md,'data',comp_temp,'title','modeled temperature','data','observed_temperature','view',3,'title','analytical temperature','view',3,'data',comp_temp-md.observed_temperature,'title','absolute error','view',3,'data',relative,'title','relative error [%]','view',3)
Index: /issm/trunk/test/NightlyRun/test1302_nightly.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1302_nightly.m	(revision 5046)
+++ /issm/trunk/test/NightlyRun/test1302_nightly.m	(revision 5046)
@@ -0,0 +1,3 @@
+field_names     ={'AdvectionTemperature'};
+field_tolerances={1e-13};
+field_values    ={comp_temp};
Index: /issm/trunk/test/NightlyRun/test1303.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1303.m	(revision 5046)
+++ /issm/trunk/test/NightlyRun/test1303.m	(revision 5046)
@@ -0,0 +1,30 @@
+% This file can be run to check that the conduction is correctly modeled.
+% There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature
+% at the lower and upper surface. The result must be a linear temperature from the upper to the lower
+% surface. if it is not the case, something is thermal modeling has been changed...
+% Just run this file in Matlab, with a properly setup Ice code. 
+
+md=model;
+md=mesh(md,'../Exp/Square.exp',100000);
+md=geography(md,'all','');
+md=parameterize(md,'../Par/SquareThermal.par');
+md=extrude(md,11,2);
+md=setelementstype(md,'Pattyn','all');
+pos1=find(md.elementonbed);     md.spctemperature(md.elements(pos1,1:3),1)=1; md.spctemperature(md.elements(pos1,1:3),2)=10;
+pos2=find(md.elementonsurface); md.spctemperature(md.elements(pos2,4:6),1)=1; md.spctemperature(md.elements(pos2,4:6),2)=0;
+md.pressure=zeros(md.numberofgrids,1);
+
+%analytical results
+%d2T/dz2=0 T(bed)=10 T(surface)=0  => T=0*(z-bed)/thickness+10*(surface-z)/thickness
+%each layer of the 3d mesh must have a constant value
+md.observed_temperature=10*(md.surface-md.z)./md.thickness;
+
+%modeled  results
+md=solve(md,'analysis_type',ThermalSolutionEnum);
+
+%plot results
+comp_temp=zeros(md.numberofgrids,1);
+comp_temp(md.results.ThermalSolution.Temperature.index)=md.results.ThermalSolution.Temperature.value;
+relative=abs((comp_temp-md.observed_temperature)./md.observed_temperature)*100;
+relative(find(comp_temp==md.observed_temperature))=0;
+plotmodel(md,'data',comp_temp,'title','modeled temperature','data','observed_temperature','view',3,'title','analytical temperature','view',3,'data',comp_temp-md.observed_temperature,'title','absolute error','view',3,'data',relative,'title','relative error [%]','view',3)
Index: /issm/trunk/test/NightlyRun/test1303_nightly.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1303_nightly.m	(revision 5046)
+++ /issm/trunk/test/NightlyRun/test1303_nightly.m	(revision 5046)
@@ -0,0 +1,3 @@
+field_names     ={'ConductionTemperature'};
+field_tolerances={1e-13};
+field_values    ={comp_temp};
Index: /issm/trunk/test/NightlyRun/test1304.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1304.m	(revision 5046)
+++ /issm/trunk/test/NightlyRun/test1304.m	(revision 5046)
@@ -0,0 +1,30 @@
+% This file can be run to check that the geothermal flux in simple conduction is correctly modeled.
+% There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature
+% at upper surface and an impose flux at its base. The result must be a linear temperature from the upper to the lower
+% surface with an imposed slope (Geothermal flux). if it is not the case, something is thermal modeling has been changed...
+
+md=model;
+md=mesh(md,'../Exp/Square.exp',100000);
+md=geography(md,'','');
+md=parameterize(md,'../Par/SquareThermal.par');
+md=extrude(md,11,1);
+md=setelementstype(md,'Pattyn','all');
+
+pos2=find(md.elementonsurface); md.spctemperature(md.elements(pos2,4:6),1)=1; md.spctemperature(md.elements(pos2,4:6),2)=0;
+md.pressure=zeros(md.numberofgrids,1);
+md.geothermalflux(:)=0.1; %100mW/m^2
+
+%analytical results
+%the result is linear with depth and is equal to 0 on the upper surface (See BC)
+%d2T/dz2=0  -k*dT/dz(bed)=G  T(surface)=0  => T=-G/k*(z-surface)
+md.observed_temperature=-0.1/md.thermalconductivity*(md.z-md.surface); %G=0.1 W/m2
+
+%modeled  results
+md=solve(md,'analysis_type',ThermalSolutionEnum);
+
+%plot results
+comp_temp=zeros(md.numberofgrids,1);
+comp_temp(md.results.ThermalSolution.Temperature.index)=md.results.ThermalSolution.Temperature.value;
+relative=abs((comp_temp-md.observed_temperature)./md.observed_temperature)*100;
+relative(find(comp_temp==md.observed_temperature))=0;
+plotmodel(md,'data',comp_temp,'title','modeled temperature','data','observed_temperature','view',3,'title','analytical temperature','view',3,'data',comp_temp-md.observed_temperature,'title','absolute error','view',3,'data',relative,'title','relative error [%]','view',3)
Index: /issm/trunk/test/NightlyRun/test1304_nightly.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1304_nightly.m	(revision 5046)
+++ /issm/trunk/test/NightlyRun/test1304_nightly.m	(revision 5046)
@@ -0,0 +1,3 @@
+field_names     ={'GeothermalFluxTemperature'};
+field_tolerances={1e-13};
+field_values    ={comp_temp};
Index: /issm/trunk/test/Par/SquareThermal.par
===================================================================
--- /issm/trunk/test/Par/SquareThermal.par	(revision 5046)
+++ /issm/trunk/test/Par/SquareThermal.par	(revision 5046)
@@ -0,0 +1,49 @@
+%Ok, start defining model parameters here
+
+%parallelization 
+md.cluster=oshostname();
+md.dt=0;
+
+disp('      creating thickness');
+h=1000;
+md.thickness=h*ones(md.numberofgrids,1);
+md.firn_layer=10*ones(md.numberofgrids,1);
+md.bed=-1000*ones(md.numberofgrids,1);
+md.surface=md.bed+md.thickness;
+
+disp('      creating velocities');
+md.vx=zeros(md.numberofgrids,1);
+md.vy=zeros(md.numberofgrids,1);
+md.vz=zeros(md.numberofgrids,1);
+
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag_coefficient=200*ones(md.numberofgrids,1); %q=1.
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag_coefficient(md.elements(pos,:))=0;
+md.drag_p=ones(md.numberofelements,1);
+md.drag_q=ones(md.numberofelements,1);
+
+disp('      creating temperatures');
+md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
+
+disp('      creating flow law paramter');
+md.rheology_B=paterson(md.observed_temperature);
+md.rheology_n=3*ones(md.numberofelements,1);
+
+disp('      creating accumulation rates');
+md.accumulation_rate=ones(md.numberofgrids,1)/md.yts; %1m/a
+md.melting_rate=0*ones(md.numberofgrids,1)/md.yts; %1m/a
+
+%Deal with boundary conditions:
+
+disp('      boundary conditions for diagnostic model');
+md=SetMarineIceSheetBC(md,'../Exp/SquareFront.exp');
+
+disp('      boundary conditions for thermal model');
+md.temperature=md.observed_temperature;
+md.spctemperature(:,1)=ones(md.numberofgrids,1); %surface temperature
+md.spctemperature(:,2)=md.observed_temperature;
+md.geothermalflux=zeros(md.numberofgrids,1); 
+pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=1*10^-3; %1 mW/m^2