diff --git a/fdelmodc3D/boundaries3D.c b/fdelmodc3D/boundaries3D.c
index 9ede67a286a0d47e97eeaef60b78a36be2b1f754..5a235970a14367854be4b5209631362d9cbbd401 100644
--- a/fdelmodc3D/boundaries3D.c
+++ b/fdelmodc3D/boundaries3D.c
@@ -1303,13 +1303,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
 					for (iz=izo; iz<ize; iz++) {
-						vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+						vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+							c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+								tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+							c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+								tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 
 						vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[ib-iz];
 					}
@@ -1326,13 +1326,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+							vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+								c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+									tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+								c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+									tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(ibz-iz)];
 						}
@@ -1348,13 +1348,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
-								vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+								vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+									c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+									c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+										tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 			
 								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
 							}
@@ -1371,13 +1371,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
-								vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+								vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+									c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+									c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+										tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 			
 								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
 							}
@@ -1404,13 +1404,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+								vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+									c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+									c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+										tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 							
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(ibz-iz)];
 						}
@@ -1426,13 +1426,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
-								vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+								vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+									c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+									c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+										tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 			
 								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
 							}
@@ -1449,13 +1449,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
-								vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+								vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+									c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+									c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+										tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 			
 								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
 							}
@@ -1476,13 +1476,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+							vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+								c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+									tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+								c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+									tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibz-iz)];
 						}
@@ -1500,13 +1500,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+							vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+								c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+									tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+								c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+									tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibz-iz)];
 						}
@@ -2570,13 +2570,14 @@ MID 	left 	mid 	mid
 #pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
 					for (iz=izo; iz<ize; iz++) {
-						vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+						vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+							c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+								tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+							c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+								tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+
 						vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[iz-ib];
 					}
 				}
@@ -2594,13 +2595,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+							vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+								c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+									tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+								c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+									tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 							
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
 						}
@@ -2617,13 +2618,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
-								vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+								vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+									c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+									c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+										tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 			
 								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
 							}
@@ -2643,13 +2644,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
-								vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+								vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+									c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+									c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+										tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 			
 								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
 							}
@@ -2677,13 +2678,13 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+							vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+								c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+									tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+								c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+									tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 							
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
 						}
@@ -2699,13 +2700,13 @@ MID 	left 	mid 	mid
 						for (iy=iyo; iy<iye; iy++) {
 							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
-								vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+								vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+									c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+									c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+										tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 			
 								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
 							}
@@ -2722,13 +2723,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
-								vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+								vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+									c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+									c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+										tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 			
 								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
 							}
@@ -2749,13 +2750,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+							vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+								c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+									tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+								c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+									tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
 						}
@@ -2773,13 +2774,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+							vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+								c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+									tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+								c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+									tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
 						}
@@ -3113,7 +3114,7 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*(
+						vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*(
 										c1*(tyy[iy*n2*n1+ix*n1+iz]     - tyy[(iy-1)*n2*n1+ix*n1+iz] +
 											tyz[iy*n2*n1+ix*n1+iz+1]   - tyz[iy*n2*n1+ix*n1+iz] +
 											txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz])  +
@@ -3167,13 +3168,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
 					for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+						vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+							c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+								tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+							c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+								tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 						
 						vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[ib-ix];
 					}
@@ -3190,13 +3191,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+							vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+								c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+									tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+								c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+									tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
 						}
@@ -3214,13 +3215,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+							vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+								c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+									tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+								c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+									tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
 						}
@@ -3619,13 +3620,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
 					for (iz=izo; iz<ize; iz++) {
-						vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+						vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+							c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+								tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+							c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+								tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 						vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[ix-ib]; 
 					}
@@ -3643,13 +3644,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+							vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+								c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+									tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+								c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+									tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)];
 						}
@@ -3668,13 +3669,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
-							vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-										c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-										c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-											tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-											txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+							vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+								c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+									tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+								c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+									tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+									txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
 						}
@@ -3839,13 +3840,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
 					for (iz=izo; iz<ize; iz++) {
-						vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-									c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-										txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-									c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-										txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+						vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+							c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+								tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+							c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+								tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 						vz[iy*n1*n2+ix*n1+iz] *= bnd.tapy[iby-iy]; 
 					}
@@ -4008,13 +4009,13 @@ MID 	left 	mid 	mid
 #pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
 					for (iz=izo; iz<ize; iz++) {
-						vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
-									c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
-										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-										txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
-									c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-										tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
-										txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
+						vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
+							c1*(tzz[iy*n2*n1+ix*n1+iz]     - tzz[iy*n2*n1+ix*n1+iz-1] +
+								tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
+							c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
+								tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
+								txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 						vz[iy*n1*n2+ix*n1+iz] *= bnd.tapy[iy-iby];
 					}
diff --git a/fdelmodc3D/elastic4dc_3D.c b/fdelmodc3D/elastic4dc_3D.c
index 8aa04a4be8d19586b9c333c54cf5684277a37e08..8c9f4ba48d43c78edabe704811e3196abd511a34 100644
--- a/fdelmodc3D/elastic4dc_3D.c
+++ b/fdelmodc3D/elastic4dc_3D.c
@@ -129,7 +129,7 @@ long elastic4dc_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l
                         tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
                         txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz])  +
                     c2*(tzz[iy*n2*n1+ix*n1+iz+1]   - tzz[iy*n2*n1+ix*n1+iz-2] +
-                        tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+                        tyz[(iy+2)*n2*n1+ix*n1+iz] - tyz[(iy-1)*n2*n1+ix*n1+iz] +
                         txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz])  );
             }
         } 
diff --git a/fdelmodc3D/getParameters3D.c b/fdelmodc3D/getParameters3D.c
index eb81ea79929b49977d42baa796613cfc7431b78d..ea9c79a05579e3638956e24bbe71460c53ae77d7 100644
--- a/fdelmodc3D/getParameters3D.c
+++ b/fdelmodc3D/getParameters3D.c
@@ -57,6 +57,7 @@ long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar
 	float xsrc1, xsrc2, ysrc1, ysrc2, zsrc1, zsrc2, tsrc1, tsrc2, tlength, tactive;
 	float src_anglex, src_angley, src_velox, src_veloy, px, py, grad2rad, rdelay, scaledt;
 	float *xsrca, *ysrca, *zsrca, rrcv;
+	float Mxx, Myy, Mzz, Mxy, Myz, Mxz;
 	float rsrc, oxsrc, oysrc, ozsrc, dphisrc, ncsrc;
 	size_t nsamp;
 	long i, j, l;
@@ -244,7 +245,7 @@ long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar
 	else {
 		dispfactor = 6;
 		//stabfactor = 0.606; /* courant number */
-		stabfactor = 0.496; //Joeri check: number changes in 3D case. Sei&Simes, 1995.
+		stabfactor = 0.495; //Joeri check: number changes in 3D case. Sei&Simes, 1995.
 		/*However as we will see in the sequel, the stability condition
 has only an indicative role. We will have to choose p = c. A t/Ax much less
 than the maximum allowed by the stability condition to fulfill the precision
@@ -1040,6 +1041,17 @@ criteria we have imposed.*/
 			case 10 : fprintf(stderr,"Moment tensor"); break;
 		}
 		fprintf(stderr,"\n");
+		if (src->type==9) {
+			Mxx = -1.0*(sin(src->dip)*cos(src->rake)*sin(2.0*src->strike)+sin(src->dip*2.0)*sin(src->rake)*sin(src->strike)*sin(src->strike));
+			Myy = sin(src->dip)*cos(src->rake)*sin(2.0*src->strike)-sin(src->dip*2.0)*sin(src->rake)*cos(src->strike)*cos(src->strike);
+			Mzz = sin(src->dip*2.0)*sin(src->rake);
+			Mxz = -1.0*(cos(src->dip)*cos(src->rake)*cos(src->strike)+cos(src->dip*2.0)*sin(src->rake)*sin(src->strike));
+			Mxy = sin(src->dip)*cos(src->rake)*cos(src->strike*2.0)+0.5*(sin(src->dip*2.0)*sin(src->rake)*sin(src->strike*2.0));
+			Myz = -1.0*(cos(src->dip)*cos(src->rake)*sin(src->strike)-cos(src->dip*2.0)*sin(src->rake)*cos(src->strike));
+			vmess("Strike %.3f (%.2f degrees) Rake %.3f (%.2f degrees) Dip %.3f (%.2f degrees)",src->strike,180.0*src->strike/M_PI,src->rake,180.0*src->rake/M_PI,src->dip,180.0*src->dip/M_PI);
+			vmess("Mxx %.2f Myy %.2f Mzz %.2f",Mxx,Myy,Mzz);
+			vmess("Mxy %.2f Mxz %.2f Myz %.2f",Mxy,Mxz,Myz);
+		}
 		if (wav->random) vmess("Wavelet has a random signature with fmax=%.2f", wav->fmax);
 		if (src->n>1) {
 			vmess("*******************************************");
diff --git a/utils/HomG.c b/utils/HomG.c
index c6da0bd4d761e0eef2f26cccea258c112daed855..d16fb269f65dcf7c43f9c28860825b1bcfd67159 100755
--- a/utils/HomG.c
+++ b/utils/HomG.c
@@ -87,23 +87,34 @@ char *sdoc[] = {
 " Optional parameters: ",
 " ",
 "   file_out= ................ Filename of the output",
+"   direction=z .............. The direction over which the data is stacked, can be x, y or z",
 "   numb= .................... integer number of first snapshot file",
 "   dnumb= ................... integer number of increment in snapshot files",
 "   zmax= .................... Integer number of maximum depth level",
-"   inx= ..................... Number of sources per depth level",
 "   zrcv= .................... z-coordinate of first receiver location",
 "   xrcv= .................... x-coordinate of first receiver location",
 "   zfps=0 ................... virtual source data are in SU format (=0) or zfp compressed (=1)",
 "   zfpr=0 ................... virtual receiver data are in SU format (=0) or zfp compressed (=1)",
-"   shift=0.0 ................ shift per shot",
-"   scheme=0 ................. Scheme used for retrieval. 0=Marchenko,",
-"                              1=Marchenko with multiple sources, 2=classical",
+"   cp=1000.0 ................ Velocity at the top of the medium in m/s",
+"   rho=1000.0 ............... Density at the top of the medium in kg/m^3",
+"   scheme=0 ................. Scheme for the retrieval",
+"   .......................... scheme=0 Marchenko homogeneous Green's function retrieval with G source",
+"   .......................... scheme=1 Marchenko homogeneous Green's function retrieval with f2 source",
+"   .......................... scheme=2 Marchenko Green's function retrieval with source depending on virtual receiver location",
+"   .......................... scheme=3 Marchenko Green's function retrieval with G source",
+"   .......................... scheme=4 Marchenko Green's function retrieval with f2 source",
+"   .......................... scheme=5 Classical homogeneous Green's function retrieval",
+"   .......................... scheme=6 Marchenko homogeneous Green's function retrieval with multiple G sources",
+"   .......................... scheme=7 Marchenko Green's function retrieval with multiple G sources",
+"   .......................... scheme=8 f1+ redatuming",
+"   .......................... scheme=9 f1- redatuming",
+"   .......................... scheme=10 2i IM(f1) redatuming",
 NULL};
 
 int main (int argc, char **argv)
 {
 	FILE    *fp_in, *fp_shot, *fp_out;
-	char    *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100];
+	char    *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], *direction;
 	float   *rcvdata, *Ghom, *shotdata, *shotdata_jkz, rho, fmin, fmax;
 	float   dt, dy, dx, t0, y0, x0, xmin, xmax1, sclsxgx, dxrcv, dyrcv, dzrcv;
 	float   *conv, *conv2, *tmp1, *tmp2, cp, shift;
@@ -132,28 +143,40 @@ int main (int argc, char **argv)
 	if (!getparlong("numb", &numb)) numb=0;
     if (!getparlong("dnumb", &dnumb)) dnumb=1;
 	if (!getparlong("scheme", &scheme)) scheme = 0;
-	if (!getparlong("ntmax", &ntmax)) ntmax = 0;
 	if (!getparlong("verbose", &verbose)) verbose = 0;
 	if (!getparlong("zfps", &zfps)) zfps = 0;
 	if (!getparlong("zfpr", &zfpr)) zfpr = 0;
+    if (!getparstring("direction", &direction)) direction = "z";
 	if (fin == NULL) verr("Incorrect vr input");
 	if (fshot == NULL) verr("Incorrect vs input");
 
     /*----------------------------------------------------------------------------*
     *   Split the filename so the number can be changed
     *----------------------------------------------------------------------------*/
-    count = dignum(numb);
+    // count = dignum(numb);
+	// if (dnumb == 0) dnumb = 1;
+	// sprintf(fins,"z%li",numb);
+	// fp_in = fopen(fin, "r");
+	// if (fp_in == NULL) {
+	// 	verr("error on opening basefile=%s", fin);
+	// }
+	// fclose(fp_in);
+	// ptr  = strstr(fin,fins);
+	// pos1 = ptr - fin;
+   	// sprintf(fbegin,"%*.*s", pos1, pos1, fin);
+   	// sprintf(fend,"%s", fin+pos1+count+1);
+
 	if (dnumb == 0) dnumb = 1;
-	sprintf(fins,"z%li",numb);
+	sprintf(fins,"%s%li",direction,numb);
 	fp_in = fopen(fin, "r");
 	if (fp_in == NULL) {
 		verr("error on opening basefile=%s", fin);
 	}
 	fclose(fp_in);
 	ptr  = strstr(fin,fins);
-	pos1 = ptr - fin;
-   	sprintf(fbegin,"%*.*s", pos1, pos1, fin);
-   	sprintf(fend,"%s", fin+pos1+count+1);
+	pos1 = ptr - fin + 1;
+   	sprintf(fbegin,"%*.*s", pos1-1, pos1-1, fin);
+   	sprintf(fend,"%s", fin+pos1+1);
 
     /*----------------------------------------------------------------------------*
     *   Determine the amount of files to be read
@@ -161,7 +184,7 @@ int main (int argc, char **argv)
 	file_det = 1;
 	nzvr=0;
 	while (file_det) {
-        sprintf(fins,"z%li",nzvr*dnumb+numb);
+        sprintf(fins,"%s%li",direction,nzvr*dnumb+numb);
         sprintf(fin,"%s%s%s",fbegin,fins,fend);
         fp_in = fopen(fin, "r");
         if (fp_in == NULL) {
@@ -189,7 +212,7 @@ int main (int argc, char **argv)
     /*----------------------------------------------------------------------------*
     *   Determine the other sizes of the files
     *----------------------------------------------------------------------------*/
-    sprintf(fins,"z%li",numb);
+    sprintf(fins,"%s%li",direction,numb);
     sprintf(fin,"%s%s%s",fbegin,fins,fend);
     if (zfpr) getVirReczfp(fin, &nxvr, &nyvr, &nxr, &nyr, &ntr);
     else getVirRec(fin, &nxvr, &nyvr, &nxr, &nyr, &ntr);
@@ -227,7 +250,7 @@ int main (int argc, char **argv)
         vmess("Number of samples for each source   : x=%li y=%li t=%li",nxvs,nyvs,ntvs);
         vmess("Sampling distance is                : x=%.3f y=%.3f t=%.3f",dx,dy,dt);
         vmess("Scaling of the transforms           : %.3f",scl);
-        vmess("Transform operators                 : fmin=%.1f fmax=%.1f cp=%.1f",fmin,fmax,cp);
+        vmess("Transform operators                 : fmin=%.1f fmax=%.1f cp=%.1f rho=%.1f",fmin,fmax,cp,rho);
     }
 
     if (ntr!=ntvs) verr("number of t-samples between virtual source (%li) and virtual receivers (%li) is not equal",ntvs,ntr);
@@ -336,7 +359,7 @@ int main (int argc, char **argv)
         }
         if (scheme==6 || scheme==8 || scheme==9 || scheme==10) tmp1 = (float *)calloc(nyr*nxr*ntr,sizeof(float));
 
-        sprintf(fins,"z%li",ir*dnumb+numb);
+        sprintf(fins,"%s%li",,direction,ir*dnumb+numb);
 		sprintf(fin2,"%s%s%s",fbegin,fins,fend);
         fp_in = fopen(fin2, "r");
 		if (fp_in == NULL) {