diff --git a/src/DMSwarm2mesh.cpp b/src/DMSwarm2mesh.cpp
index 8c5cde4..28cba09 100644
--- a/src/DMSwarm2mesh.cpp
+++ b/src/DMSwarm2mesh.cpp
@@ -37,6 +37,9 @@ extern PetscScalar *conductivity;
 
 extern PetscInt periodic_boundary;
 
+extern PetscInt WITH_SHEAR_H;
+extern PetscReal Xi_shear;
+
 extern double visc_MIN;
 
 extern PetscReal air_threshold_density;
@@ -45,6 +48,8 @@ extern PetscReal rho0_scaled;
 
 extern PetscReal epsilon_x;
 
+//Shear heating only in 2D!!! Must be implemented for 3D!!!
+
 PetscErrorCode Swarm2Mesh_2d(){
 
 	PetscErrorCode ierr;
@@ -103,6 +108,8 @@ PetscErrorCode Swarm2Mesh_2d(){
 	PetscReal *strain_rate_fac;
 	PetscInt *layer_array;
 
+	PetscReal inter_H_aux;
+
 	ierr = DMSwarmGetLocalSize(dms,&nlocal);CHKERRQ(ierr);
 
 	ierr = DMSwarmGetField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);CHKERRQ(ierr);
@@ -152,33 +159,38 @@ PetscErrorCode Swarm2Mesh_2d(){
 		if (rx<0 || rx>1) {printf("weird rx=%f , Swarm2Mesh\n",rx); exit(1);}
 		if (rz<0 || rz>1) {printf("weird rz=%f , Swarm2Mesh\n",rz); exit(1);}
 
-
+		if (WITH_SHEAR_H == 1){
+			inter_H_aux = inter_H[layer_array[p]] + Xi_shear*4*geoq_fac[p]*strain_rate_fac[p]*strain_rate_fac[p]/inter_rho[layer_array[p]];
+		}
+		else {
+			inter_H_aux = inter_H[layer_array[p]];
+		}
 
 
 		rfac = (1.0-rx)*(1.0-rz);
 		qq_rho	[k][i] += rfac*inter_rho[layer_array[p]];
-		qq_H	[k][i] += rfac*inter_H[layer_array[p]];
+		qq_H	[k][i] += rfac*inter_H_aux;
 		qq_strain[k][i] += rfac*strain_fac[p];
 		qq_strain_rate[k][i] += rfac*strain_rate_fac[p];
 		qq_cont	[k][i] += rfac;
 
 		rfac = (rx)*(1.0-rz);
 		qq_rho	[k][i+1] += rfac*inter_rho[layer_array[p]];
-		qq_H	[k][i+1] += rfac*inter_H[layer_array[p]];
+		qq_H	[k][i+1] += rfac*inter_H_aux;
 		qq_strain[k][i+1] += rfac*strain_fac[p];
 		qq_strain_rate[k][i+1] += rfac*strain_rate_fac[p];
 		qq_cont	[k][i+1] += rfac;
 
 		rfac = (1.0-rx)*(rz);
 		qq_rho	[k+1][i] += rfac*inter_rho[layer_array[p]];
-		qq_H	[k+1][i] += rfac*inter_H[layer_array[p]];
+		qq_H	[k+1][i] += rfac*inter_H_aux;
 		qq_strain[k+1][i] += rfac*strain_fac[p];
 		qq_strain_rate[k+1][i] += rfac*strain_rate_fac[p];
 		qq_cont	[k+1][i] += rfac;
 
 		rfac = (rx)*(rz);
 		qq_rho	[k+1][i+1] += rfac*inter_rho[layer_array[p]];
-		qq_H	[k+1][i+1] += rfac*inter_H[layer_array[p]];
+		qq_H	[k+1][i+1] += rfac*inter_H_aux;
 		qq_strain[k+1][i+1] += rfac*strain_fac[p];
 		qq_strain_rate[k+1][i+1] += rfac*strain_rate_fac[p];
 		qq_cont	[k+1][i+1] += rfac;
diff --git a/src/header.h b/src/header.h
index 176eff1..2936e73 100644
--- a/src/header.h
+++ b/src/header.h
@@ -27,6 +27,8 @@ PetscScalar sp_d_c = 0.0;
 // Parameter file boolean variables
 PetscInt WITH_NON_LINEAR = 0; // 1=True, 0=False
 PetscInt WITH_ADIABATIC_H = 0; // 1=True, 0=False
+PetscInt WITH_SHEAR_H = 0; // 1=True, 0=False
+PetscReal Xi_shear = 1; //scale factor
 PetscInt WITH_RADIOGENIC_H = 0; // 1=True, 0=False
 PetscInt PLASTICITY = 1; // 1=True, 0=False
 PetscInt direct_solver = 1; // 1=direct, 0=iterative
diff --git a/src/reader.cpp b/src/reader.cpp
index 528bfaa..788923f 100644
--- a/src/reader.cpp
+++ b/src/reader.cpp
@@ -53,6 +53,8 @@ extern int bcT_front;
 extern int bcT_back;
 extern PetscInt WITH_NON_LINEAR;
 extern PetscInt WITH_ADIABATIC_H;
+extern PetscInt WITH_SHEAR_H;
+extern PetscReal Xi_shear;
 extern PetscInt WITH_RADIOGENIC_H;
 extern PetscInt PLASTICITY;
 extern PetscReal denok_min;
@@ -268,6 +270,8 @@ PetscErrorCode reader(int rank, const char fName[]){
 			else if (strcmp(tkn_w, "geoq") == 0) {geoq_on = check_a_b(tkn_w, tkn_v, "on", "off");}
 			else if (strcmp(tkn_w, "non_linear_method") == 0) {WITH_NON_LINEAR = check_a_b(tkn_w, tkn_v, "on", "off");}
 			else if (strcmp(tkn_w, "adiabatic_component") == 0) {WITH_ADIABATIC_H = check_a_b(tkn_w, tkn_v, "on", "off");}
+			else if (strcmp(tkn_w, "shear_heating") == 0) {WITH_SHEAR_H = check_a_b(tkn_w, tkn_v, "on", "off");}
+			else if (strcmp(tkn_w, "Xi_shear") == 0) {Xi_shear = atof(tkn_v);}
 			else if (strcmp(tkn_w, "radiogenic_component") == 0) {WITH_RADIOGENIC_H = check_a_b(tkn_w, tkn_v, "on", "off");}
 			else if (strcmp(tkn_w, "plasticity") == 0) {PLASTICITY = check_a_b(tkn_w, tkn_v, "on", "off");}
 			else if (strcmp(tkn_w, "top_normal_velocity") == 0) {bcv_top_normal = check_a_b(tkn_w, tkn_v, "fixed", "free");}
@@ -409,6 +413,8 @@ PetscErrorCode reader(int rank, const char fName[]){
 	MPI_Bcast(&c_heat_capacity,1,MPI_DOUBLE,0,PETSC_COMM_WORLD);
 	MPI_Bcast(&WITH_NON_LINEAR,1,MPI_INT,0,PETSC_COMM_WORLD);
 	MPI_Bcast(&WITH_ADIABATIC_H,1,MPI_INT,0,PETSC_COMM_WORLD);
+	MPI_Bcast(&WITH_SHEAR_H,1,MPI_INT,0,PETSC_COMM_WORLD);
+	MPI_Bcast(&Xi_shear,1,MPIU_REAL,0,PETSC_COMM_WORLD);
 	MPI_Bcast(&WITH_RADIOGENIC_H,1,MPI_INT,0,PETSC_COMM_WORLD);
 	MPI_Bcast(&PLASTICITY,1,MPI_INT,0,PETSC_COMM_WORLD);
 	MPI_Bcast(&bcv_top_normal,1,MPI_INT,0,PETSC_COMM_WORLD);