diff --git a/fdelmodc3D/boundaries3D.c b/fdelmodc3D/boundaries3D.c
index 3a3f2cb2f7a7cc245c654ea4163175ab5b0f38d5..1496fbb7388ce82d83eb4a0913ee289ac8c1b417 100644
--- a/fdelmodc3D/boundaries3D.c
+++ b/fdelmodc3D/boundaries3D.c
@@ -79,7 +79,7 @@ MID 	left 	mid 	mid
     dy  = mod.dy;
     dt  = mod.dt;
     fac = dt/dx;
-    if ( (bnd.top==2) || (bnd.bot==2) || (bnd.lef==2) || (bnd.rig==2) ) pml=1;
+    if ( (bnd.top==2) || (bnd.bot==2) || (bnd.lef==2) || (bnd.rig==2) || (bnd.fro==2) || (bnd.bac==2) ) pml=1;
 	else pml=0;
 
 	ibnd = mod.iorder/2-1;
@@ -104,7 +104,6 @@ MID 	left 	mid 	mid
 
 	if (bnd.top==3) { /* rigid surface at top */
 #pragma omp for private (ix, iy, iz) nowait
-//#pragma ivdep
         for (iy=1; iy<=ny; iy++) {
         	#pragma ivdep
             for (ix=1; ix<=nx; ix++) {
@@ -118,7 +117,6 @@ MID 	left 	mid 	mid
 	}
 	if (bnd.rig==3) { /* rigid surface at right */
 #pragma omp for private (ix, iy, iz) nowait
-//#pragma ivdep
         for (iy=1; iy<=ny; iy++) {
         	#pragma ivdep
             for (iz=1; iz<=nz; iz++) {
@@ -222,7 +220,6 @@ MID 	left 	mid 	mid
 			ibz = (bnd.ntap+izo-1);
 #pragma omp for private(ix,iy,iz)
 			for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
 					#pragma ivdep
 					for (iz=izo; iz<ize; iz++) {
@@ -244,7 +241,6 @@ MID 	left 	mid 	mid
 				ibx = (ixo);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
@@ -265,7 +261,6 @@ MID 	left 	mid 	mid
 					iby = (bnd.ntap+iyo-1);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
@@ -287,7 +282,6 @@ MID 	left 	mid 	mid
 					iby = (iyo);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
@@ -317,7 +311,6 @@ MID 	left 	mid 	mid
 				ibx = (bnd.ntap+ixo-1);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
@@ -337,7 +330,6 @@ MID 	left 	mid 	mid
 					iby = (bnd.ntap+iyo-1);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
@@ -359,7 +351,6 @@ MID 	left 	mid 	mid
 					iby = (iyo);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
@@ -385,7 +376,6 @@ MID 	left 	mid 	mid
 				iby = (bnd.ntap+iyo-1);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
@@ -408,7 +398,6 @@ MID 	left 	mid 	mid
 				iby = (iyo);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
@@ -435,7 +424,6 @@ MID 	left 	mid 	mid
 			ib = (bnd.ntap+izo-1);
 #pragma omp for private(ix,iy,iz)
 			for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
 					#pragma ivdep
 					for (iz=izo; iz<ize; iz++) {
@@ -457,7 +445,6 @@ MID 	left 	mid 	mid
 				ibx = (ixo);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
@@ -476,7 +463,6 @@ MID 	left 	mid 	mid
 					iby = (bnd.ntap+iyo-1);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
@@ -499,7 +485,6 @@ MID 	left 	mid 	mid
 					iby = (iyo);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
 							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
@@ -815,92 +800,109 @@ MID 	left 	mid 	mid
 				}
 			}
 		}
-		
-	}
 
-	/*********/
-	/* Bottom */
-	/*********/
-	if (bnd.bot==4) {
-		
-		if (mod.ischeme <= 2) { /* Acoustic scheme */
+		else { /* Elastic scheme */
 			
 			/* Vx field */
-			/* mid bottom mid vx*/
+			/* mid top mid vx */
 			ixo = mod.ioXx;
 			ixe = mod.ieXx;
 			iyo = mod.ioXy;
 			iye = mod.ieXy;
-			izo = mod.ieXz;
-			ize = mod.ieXz+bnd.ntap;
-			
-			ib = (ize-bnd.ntap);
+			izo = mod.ioXz-bnd.ntap;
+			ize = mod.ioXz;
+	
+			ibz = (bnd.ntap+izo-1);
 #pragma omp for private(ix,iy,iz)
 			for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
+					#pragma ivdep
 					for (iz=izo; iz<ize; iz++) {
 						vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-									c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-									c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 
-						vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[iz-ib];
+						vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[ibz-iz];
 					}
 				}
 			}
-			/* right bottom mid corner vx */
+
+
+			/* right top mid corner vx */
 			if (bnd.rig==4) {
 				ixo = mod.ieXx;
 				ixe = mod.ieXx+bnd.ntap;
-				ibz = (izo);
+				ibz = (bnd.ntap+izo-1);
 				ibx = (ixo);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
+						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
 							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 		
-							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
-				/* right bottom front corner vx */
+
+
+				/* right top front corner vx */
 				if (bnd.fro==4) {
 					iyo = mod.ioXy-bnd.ntap;
 					iye = mod.ioXy;
 					iby = (bnd.ntap+iyo-1);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
+							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
 								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-											c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 			
-								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
 				}
-				/* right bottom back corner vx */
+		
+		
+				/* right top back corner vx */
 				if (bnd.bac==4) {
 					iyo = mod.ieXy;
 					iye = mod.ieXy+bnd.ntap;
 					iby = (iyo);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
+							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
 								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-											c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 			
-								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
@@ -911,187 +913,244 @@ MID 	left 	mid 	mid
 			ixe = mod.ieXx;
 			iyo = mod.ioXy;
 			iye = mod.ieXy;
-			izo = mod.ieXz;
-			ize = mod.ieXz+bnd.ntap;
+			izo = mod.ioXz-bnd.ntap;
+			ize = mod.ioXz;
 
-			/* left bottom mid corner vx*/
+			/* left top mid corner vx */
 			if (bnd.lef==4) {
 				ixo = mod.ioXx-bnd.ntap;
 				ixe = mod.ioXx;
-				ibz = (izo);
+				ibz = (bnd.ntap+izo-1);
 				ibx = (bnd.ntap+ixo-1);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
+						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
 							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 							
-							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
-				/* left bottom front corner vx */
+
+				/* left top front corner vx */
 				if (bnd.fro==4) {
 					iyo = mod.ioXy-bnd.ntap;
 					iye = mod.ioXy;
 					iby = (bnd.ntap+iyo-1);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
+							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
 								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-											c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 			
-								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
 				}
-				/* left bottom back corner vx*/
+
+
+				/* left top back corner vx */
 				if (bnd.bac==4) {
 					iyo = mod.ieXy;
 					iye = mod.ieXy+bnd.ntap;
 					iby = (iyo);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
+							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
 								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-											c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 			
-								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
 				}
 			}
-			/* mid bottom front corner vx */
+
+
+			/* mid top front corner vx */
 			if (bnd.fro==4) {
 				ixo = mod.ioXx;
 				ixe = mod.ieXx;
 				iyo = mod.ioXy-bnd.ntap;
 				iye = mod.ioXy;
-				ibz = (izo);
+				ibz = (bnd.ntap+izo-1);
 				iby = (bnd.ntap+iyo-1);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
+						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
 							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 		
-							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
 			}
-			/*  mid bottom back corner vx */
+
+
+			/* mid top back corner vx */
 			if (bnd.bac==4) {
 				iyo = mod.ieXy;
 				iye = mod.ieXy+bnd.ntap;
-				ibz = (izo);
+				ibz = (bnd.ntap+izo-1);
 				iby = (iyo);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
+						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
 							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 		
-							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
-			}
-
+			}		
+	
 
 			/* Vy field */
-			/* mid bottom mid vy */
+			/* mid top mid vy */
 			ixo = mod.ioYx;
 			ixe = mod.ieYx;
 			iyo = mod.ioYy;
 			iye = mod.ieYy;
-			izo = mod.ieYz;
-			ize = mod.ieYz+bnd.ntap;
+			izo = mod.ioYz-bnd.ntap;
+			ize = mod.ioYz;
 	
-			ib = (ize-bnd.ntap);
+			ib = (bnd.ntap+izo-1);
 #pragma omp for private(ix,iy,iz)
 			for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
+					#pragma ivdep
 					for (iz=izo; iz<ize; iz++) {
 						vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 
-						vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapy[iz-ib];
+						vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapy[ib-iz];
 					}
 				}
 			}
-			/* right bottom mid corner vy */
+	
+
+			/* right top mid corner vy */
 			if (bnd.rig==4) {
 				ixo = mod.ieYx;
-				ixe = ixo+bnd.ntap;
-				ibz = (izo);
+				ixe = mod.ieYx+bnd.ntap;
+				ibz = (bnd.ntap+izo-1);
 				ibx = (ixo);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
+						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
 							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 		
-							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
-				/* right bottom front corner vy */
+				/* right top front corner vy */
 				if (bnd.fro==4) {
 					iyo = mod.ioYy-bnd.ntap;
 					iye = mod.ioYy;
 					iby = (bnd.ntap+iyo-1);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
+							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
 								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-											c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 			
-								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
 				}
-				/* right bottom back corner vy */
+
+	
+
+				/* right top back corner vy */
 				if (bnd.bac==4) {
 					iyo = mod.ieYy;
 					iye = mod.ieYy+bnd.ntap;
 					iby = (iyo);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
+							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
 								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-											c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 			
-								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
@@ -1101,14 +1160,15 @@ MID 	left 	mid 	mid
 			ixe = mod.ieYx;
 			iyo = mod.ioYy;
 			iye = mod.ieYy;
-			izo = mod.ieYz;
-			ize = mod.ieYz+bnd.ntap;
+			izo = mod.ioYz-bnd.ntap;
+			ize = mod.ioYz;
 
-			/* left bottom mid corner vy */
+
+			/* left top mid corner vy */
 			if (bnd.lef==4) {
 				ixo = mod.ioYx-bnd.ntap;
 				ixe = mod.ioYx;
-				ibz = (izo);
+				ibz = (bnd.ntap+izo-1);
 				ibx = (bnd.ntap+ixo-1);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
@@ -1116,14 +1176,19 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 							
-							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
-				/* left bottom front corner vy */
+
+				/* left top front corner vy */
 				if (bnd.fro==4) {
 					iyo = mod.ioYy-bnd.ntap;
 					iye = mod.ioYy;
@@ -1134,15 +1199,20 @@ MID 	left 	mid 	mid
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
 								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-											c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 			
-								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
 				}
-				/* left bottom back corner vy */
+
+				/* left top back corner vy */
 				if (bnd.bac==4) {
 					iyo = mod.ieYy;
 					iye = mod.ieYy+bnd.ntap;
@@ -1153,22 +1223,26 @@ MID 	left 	mid 	mid
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
 								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-											c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 			
-								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
 				}
 			}
-			/* mid bottom front corner vy */
+			/* front top mid corner vy */
 			if (bnd.fro==4) {
 				ixo = mod.ioYx;
 				ixe = mod.ieYx;
 				iyo = mod.ioYy-bnd.ntap;
 				iye = mod.ioYy;
-				ibz = (izo);
+				ibz = (bnd.ntap+izo-1);
 				iby = (bnd.ntap+iyo-1);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
@@ -1176,19 +1250,23 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 		
-							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
 			}
-			/* mid bottom back corner vy */
+			/* back top mid corner vy */
 			if (bnd.bac==4) {
 				iyo = mod.ieYy;
 				iye = mod.ieYy+bnd.ntap;
-				ibz = (izo);
+				ibz = (bnd.ntap+izo-1);
 				iby = (iyo);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
@@ -1196,10 +1274,14 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 		
-							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
@@ -1207,49 +1289,56 @@ MID 	left 	mid 	mid
 
 
 			/* Vz field */
+			/* mid top mid vz*/
 			ixo = mod.ioZx;
 			ixe = mod.ieZx;
 			iyo = mod.ioZy;
 			iye = mod.ieZy;
-			izo = mod.ieZz;
-			ize = mod.ieZz+bnd.ntap;
-			
-			ib = (ize-bnd.ntap);
+			izo = mod.ioZz-bnd.ntap;
+			ize = mod.ioZz;
+	
+			ib = (bnd.ntap+izo-1);
 #pragma omp for private (ix, iy, iz) 
 			for (ix=ixo; ix<ixe; ix++) {
 #pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
 					for (iz=izo; iz<ize; iz++) {
 						vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-									c1*(tzz[iy*n1*n2+ix*n1+iz]   - tzz[iy*n1*n2+ix*n1+iz-1]) +
-									c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
-						vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[iz-ib];
+										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*n1*n2+ix*n1+iz] *= bnd.tapz[ib-iz];
 					}
 				}
 			}
-
-
-			/* right bottom corner */
+			/* right top mid corner vz */
 			if (bnd.rig==4) {
 				ixo = mod.ieZx;
 				ixe = ixo+bnd.ntap;
-				ibz = (izo);
+				ibz = (bnd.ntap+izo-1);
 				ibx = (ixo);
-#pragma omp for private(ix,iy,iz)
+#pragma omp for private(ix,iz)
 				for (ix=ixo; ix<ixe; ix++) {
 #pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]   - tzz[iy*n1*n2+ix*n1+iz-1]) +
-										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
-							
-							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+										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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
-
-				/* right bottom front corner */
+				/* right top front corner vz */
 				if (bnd.fro==4) {
 					iyo = mod.ioZy-bnd.ntap;
 					iye = mod.ioZy;
@@ -1260,18 +1349,19 @@ MID 	left 	mid 	mid
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
 								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
-											c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+										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*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
 				}
-
-
-
-				/* right bottom back corner */
+				/* right top back corner vz */
 				if (bnd.bac==4) {
 					iyo = mod.ieZy;
 					iye = mod.ieZy+bnd.ntap;
@@ -1282,63 +1372,74 @@ MID 	left 	mid 	mid
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
 								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
-											c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+										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*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
 				}
 			}
-			
+
 			ixo = mod.ioZx;
 			ixe = mod.ieZx;
 			iyo = mod.ioZy;
 			iye = mod.ieZy;
-			izo = mod.ieZz;
-			ize = mod.ieZz+bnd.ntap;
+			izo = mod.ioZz-bnd.ntap;
+			ize = mod.ioZz;
 
-			/* left bottom corner */
+			/* left top mid corner vz */
 			if (bnd.lef==4) {
 				ixo = mod.ioZx-bnd.ntap;
 				ixe = mod.ioZx;
-				ibz = (izo);
+				ibz = (bnd.ntap+izo-1);
 				ibx = (bnd.ntap+ixo-1);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
+#pragma ivdep
 					for (iy=iyo; iy<iye; iy++) {
-						#pragma ivdep
 						for (iz=izo; iz<ize; iz++) {
 							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]   - tzz[iy*n1*n2+ix*n1+iz-1]) +
-										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+										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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
-				/* left bottom front corner */
+				/* left top front corner vz */
 				if (bnd.fro==4) {
 					iyo = mod.ioZy-bnd.ntap;
 					iye = mod.ioZy;
 					iby = (bnd.ntap+iyo-1);
 #pragma omp for private(ix,iy,iz)
 					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
 						for (iy=iyo; iy<iye; iy++) {
-							#pragma ivdep
 							for (iz=izo; iz<ize; iz++) {
 								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
-											c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+										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*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
 				}
-				/* left bottom back corner */
+				/* left top back corner vz */
 				if (bnd.bac==4) {
 					iyo = mod.ieZy;
 					iye = mod.ieZy+bnd.ntap;
@@ -1349,22 +1450,26 @@ MID 	left 	mid 	mid
 						for (iy=iyo; iy<iye; iy++) {
 							for (iz=izo; iz<ize; iz++) {
 								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-											c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
-											c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+										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*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
 							}
 						}
 					}
 				}
 			}
-			/* front bottom corner */
+			/* mid top front corner vz */
 			if (bnd.fro==4) {
 				ixo = mod.ioZx;
 				ixe = mod.ieZx;
 				iyo = mod.ioZy-bnd.ntap;
 				iye = mod.ioZy;
-				ibz = (izo);
+				ibz = (bnd.ntap+izo-1);
 				iby = (bnd.ntap+iyo-1);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
@@ -1372,19 +1477,23 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
-										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+										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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
 			}
-			/* Back bottom corner */
+			/* mid top back corner vz */
 			if (bnd.bac==4) {
 				iyo = mod.ieZy;
 				iye = mod.ieZy+bnd.ntap;
-				ibz = (izo);
+				ibz = (bnd.ntap+izo-1);
 				iby = (iyo);
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
@@ -1392,37 +1501,1753 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
-										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+										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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibz-iz)];
 						}
 					}
 				}
 			}
-			
 		}
+		
 	}
 
+	/*********/
+	/* Bottom */
+	/*********/
+	if (bnd.bot==4) {
+		
+		if (mod.ischeme <= 2) { /* Acoustic scheme */
+			
+			/* Vx field */
+			/* mid bottom mid vx*/
+			ixo = mod.ioXx;
+			ixe = mod.ieXx;
+			iyo = mod.ioXy;
+			iye = mod.ieXy;
+			izo = mod.ieXz;
+			ize = mod.ieXz+bnd.ntap;
+			
+			ib = (ize-bnd.ntap);
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+									c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+									c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+
+						vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[iz-ib];
+					}
+				}
+			}
+			/* right bottom mid corner vx */
+			if (bnd.rig==4) {
+				ixo = mod.ieXx;
+				ixe = mod.ieXx+bnd.ntap;
+				ibz = (izo);
+				ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+		
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+				/* right bottom front corner vx */
+				if (bnd.fro==4) {
+					iyo = mod.ioXy-bnd.ntap;
+					iye = mod.ioXy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+											c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+			
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+				/* right bottom back corner vx */
+				if (bnd.bac==4) {
+					iyo = mod.ieXy;
+					iye = mod.ieXy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+											c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+			
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+
+			ixo = mod.ioXx;
+			ixe = mod.ieXx;
+			iyo = mod.ioXy;
+			iye = mod.ieXy;
+			izo = mod.ieXz;
+			ize = mod.ieXz+bnd.ntap;
+
+			/* left bottom mid corner vx*/
+			if (bnd.lef==4) {
+				ixo = mod.ioXx-bnd.ntap;
+				ixe = mod.ioXx;
+				ibz = (izo);
+				ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+							
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+				/* left bottom front corner vx */
+				if (bnd.fro==4) {
+					iyo = mod.ioXy-bnd.ntap;
+					iye = mod.ioXy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+											c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+			
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+				/* left bottom back corner vx*/
+				if (bnd.bac==4) {
+					iyo = mod.ieXy;
+					iye = mod.ieXy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+											c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+			
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+			/* mid bottom front corner vx */
+			if (bnd.fro==4) {
+				ixo = mod.ioXx;
+				ixe = mod.ieXx;
+				iyo = mod.ioXy-bnd.ntap;
+				iye = mod.ioXy;
+				ibz = (izo);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+		
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+			/*  mid bottom back corner vx */
+			if (bnd.bac==4) {
+				iyo = mod.ieXy;
+				iye = mod.ieXy+bnd.ntap;
+				ibz = (izo);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+		
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+
+
+			/* Vy field */
+			/* mid bottom mid vy */
+			ixo = mod.ioYx;
+			ixe = mod.ieYx;
+			iyo = mod.ioYy;
+			iye = mod.ieYy;
+			izo = mod.ieYz;
+			ize = mod.ieYz+bnd.ntap;
+	
+			ib = (ize-bnd.ntap);
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+
+						vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapy[iz-ib];
+					}
+				}
+			}
+			/* right bottom mid corner vy */
+			if (bnd.rig==4) {
+				ixo = mod.ieYx;
+				ixe = ixo+bnd.ntap;
+				ibz = (izo);
+				ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+		
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+				/* right bottom front corner vy */
+				if (bnd.fro==4) {
+					iyo = mod.ioYy-bnd.ntap;
+					iye = mod.ioYy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+											c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+			
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+				/* right bottom back corner vy */
+				if (bnd.bac==4) {
+					iyo = mod.ieYy;
+					iye = mod.ieYy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+											c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+			
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+			ixo = mod.ioYx;
+			ixe = mod.ieYx;
+			iyo = mod.ioYy;
+			iye = mod.ieYy;
+			izo = mod.ieYz;
+			ize = mod.ieYz+bnd.ntap;
+
+			/* left bottom mid corner vy */
+			if (bnd.lef==4) {
+				ixo = mod.ioYx-bnd.ntap;
+				ixe = mod.ioYx;
+				ibz = (izo);
+				ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+							
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+				/* left bottom front corner vy */
+				if (bnd.fro==4) {
+					iyo = mod.ioYy-bnd.ntap;
+					iye = mod.ioYy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+											c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+			
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+				/* left bottom back corner vy */
+				if (bnd.bac==4) {
+					iyo = mod.ieYy;
+					iye = mod.ieYy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+											c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+			
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+			/* mid bottom front corner vy */
+			if (bnd.fro==4) {
+				ixo = mod.ioYx;
+				ixe = mod.ieYx;
+				iyo = mod.ioYy-bnd.ntap;
+				iye = mod.ioYy;
+				ibz = (izo);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+		
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+			/* mid bottom back corner vy */
+			if (bnd.bac==4) {
+				iyo = mod.ieYy;
+				iye = mod.ieYy+bnd.ntap;
+				ibz = (izo);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+		
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+
+
+			/* Vz field */
+			ixo = mod.ioZx;
+			ixe = mod.ieZx;
+			iyo = mod.ioZy;
+			iye = mod.ieZy;
+			izo = mod.ieZz;
+			ize = mod.ieZz+bnd.ntap;
+			
+			ib = (ize-bnd.ntap);
+#pragma omp for private (ix, iy, iz) 
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+									c1*(tzz[iy*n1*n2+ix*n1+iz]   - tzz[iy*n1*n2+ix*n1+iz-1]) +
+									c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+						vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[iz-ib];
+					}
+				}
+			}
+
+
+			/* right bottom corner */
+			if (bnd.rig==4) {
+				ixo = mod.ieZx;
+				ixe = ixo+bnd.ntap;
+				ibz = (izo);
+				ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]   - tzz[iy*n1*n2+ix*n1+iz-1]) +
+										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+							
+							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+
+				/* right bottom front corner */
+				if (bnd.fro==4) {
+					iyo = mod.ioZy-bnd.ntap;
+					iye = mod.ioZy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
+											c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+			
+								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+
+
+
+				/* right bottom back corner */
+				if (bnd.bac==4) {
+					iyo = mod.ieZy;
+					iye = mod.ieZy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
+											c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+			
+								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+			
+			ixo = mod.ioZx;
+			ixe = mod.ieZx;
+			iyo = mod.ioZy;
+			iye = mod.ieZy;
+			izo = mod.ieZz;
+			ize = mod.ieZz+bnd.ntap;
+
+			/* left bottom corner */
+			if (bnd.lef==4) {
+				ixo = mod.ioZx-bnd.ntap;
+				ixe = mod.ioZx;
+				ibz = (izo);
+				ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+//#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						#pragma ivdep
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]   - tzz[iy*n1*n2+ix*n1+iz-1]) +
+										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+							
+							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+				/* left bottom front corner */
+				if (bnd.fro==4) {
+					iyo = mod.ioZy-bnd.ntap;
+					iye = mod.ioZy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+						for (iy=iyo; iy<iye; iy++) {
+							#pragma ivdep
+							for (iz=izo; iz<ize; iz++) {
+								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
+											c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+			
+								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+				/* left bottom back corner */
+				if (bnd.bac==4) {
+					iyo = mod.ieZy;
+					iye = mod.ieZy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+											c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
+											c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+			
+								vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+			/* front bottom corner */
+			if (bnd.fro==4) {
+				ixo = mod.ioZx;
+				ixe = mod.ieZx;
+				iyo = mod.ioZy-bnd.ntap;
+				iye = mod.ioZy;
+				ibz = (izo);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
+										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+		
+							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+			/* Back bottom corner */
+			if (bnd.bac==4) {
+				iyo = mod.ieZy;
+				iye = mod.ieZy+bnd.ntap;
+				ibz = (izo);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
+										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+		
+							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+			
+		}
+		else { /* Elastic scheme */
+			
+			/* Vx field */
+			/* mid bottom mid vx*/
+			ixo = mod.ioXx;
+			ixe = mod.ieXx;
+			iyo = mod.ioXy;
+			iye = mod.ieXy;
+			izo = mod.ieXz;
+			ize = mod.ieXz+bnd.ntap;
+			
+			ib = (ize-bnd.ntap);
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+
+						vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[iz-ib];
+					}
+				}
+			}
+			/* right bottom mid corner vx */
+			if (bnd.rig==4) {
+				ixo = mod.ieXx;
+				ixe = mod.ieXx+bnd.ntap;
+				ibz = (izo);
+				ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+		
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+				/* right bottom front corner vx */
+				if (bnd.fro==4) {
+					iyo = mod.ioXy-bnd.ntap;
+					iye = mod.ioXy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+			
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+				/* right bottom back corner vx */
+				if (bnd.bac==4) {
+					iyo = mod.ieXy;
+					iye = mod.ieXy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+			
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+
+			ixo = mod.ioXx;
+			ixe = mod.ieXx;
+			iyo = mod.ioXy;
+			iye = mod.ieXy;
+			izo = mod.ieXz;
+			ize = mod.ieXz+bnd.ntap;
+
+			/* left bottom mid corner vx*/
+			if (bnd.lef==4) {
+				ixo = mod.ioXx-bnd.ntap;
+				ixe = mod.ioXx;
+				ibz = (izo);
+				ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+							
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+				/* left bottom front corner vx */
+				if (bnd.fro==4) {
+					iyo = mod.ioXy-bnd.ntap;
+					iye = mod.ioXy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+			
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+				/* left bottom back corner vx*/
+				if (bnd.bac==4) {
+					iyo = mod.ieXy;
+					iye = mod.ieXy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+			
+								vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+			/* mid bottom front corner vx */
+			if (bnd.fro==4) {
+				ixo = mod.ioXx;
+				ixe = mod.ieXx;
+				iyo = mod.ioXy-bnd.ntap;
+				iye = mod.ioXy;
+				ibz = (izo);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+		
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+			/*  mid bottom back corner vx */
+			if (bnd.bac==4) {
+				iyo = mod.ieXy;
+				iye = mod.ieXy+bnd.ntap;
+				ibz = (izo);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+		
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+
+
+			/* Vy field */
+			/* mid bottom mid vy */
+			ixo = mod.ioYx;
+			ixe = mod.ieYx;
+			iyo = mod.ioYy;
+			iye = mod.ieYy;
+			izo = mod.ieYz;
+			ize = mod.ieYz+bnd.ntap;
+	
+			ib = (ize-bnd.ntap);
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+
+						vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapy[iz-ib];
+					}
+				}
+			}
+			/* right bottom mid corner vy */
+			if (bnd.rig==4) {
+				ixo = mod.ieYx;
+				ixe = ixo+bnd.ntap;
+				ibz = (izo);
+				ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+		
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+				/* right bottom front corner vy */
+				if (bnd.fro==4) {
+					iyo = mod.ioYy-bnd.ntap;
+					iye = mod.ioYy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+			
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+				/* right bottom back corner vy */
+				if (bnd.bac==4) {
+					iyo = mod.ieYy;
+					iye = mod.ieYy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+			
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+			ixo = mod.ioYx;
+			ixe = mod.ieYx;
+			iyo = mod.ioYy;
+			iye = mod.ieYy;
+			izo = mod.ieYz;
+			ize = mod.ieYz+bnd.ntap;
+
+			/* left bottom mid corner vy */
+			if (bnd.lef==4) {
+				ixo = mod.ioYx-bnd.ntap;
+				ixe = mod.ioYx;
+				ibz = (izo);
+				ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+							
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+				/* left bottom front corner vy */
+				if (bnd.fro==4) {
+					iyo = mod.ioYy-bnd.ntap;
+					iye = mod.ioYy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+			
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+				/* left bottom back corner vy */
+				if (bnd.bac==4) {
+					iyo = mod.ieYy;
+					iye = mod.ieYy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+			
+								vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+			/* mid bottom front corner vy */
+			if (bnd.fro==4) {
+				ixo = mod.ioYx;
+				ixe = mod.ieYx;
+				iyo = mod.ioYy-bnd.ntap;
+				iye = mod.ioYy;
+				ibz = (izo);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+		
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+			/* mid bottom back corner vy */
+			if (bnd.bac==4) {
+				iyo = mod.ieYy;
+				iye = mod.ieYy+bnd.ntap;
+				ibz = (izo);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+		
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+
+
+			/* Vz field */
+			ixo = mod.ioZx;
+			ixe = mod.ieZx;
+			iyo = mod.ioZy;
+			iye = mod.ieZy;
+			izo = mod.ieZz;
+			ize = mod.ieZz+bnd.ntap;
+			
+			ib = (ize-bnd.ntap);
+#pragma omp for private (ix, iy, iz) 
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz] *= bnd.tapz[iz-ib];
+					}
+				}
+			}
+
+
+			/* right bottom corner */
+			if (bnd.rig==4) {
+				ixo = mod.ieZx;
+				ixe = ixo+bnd.ntap;
+				ibz = (izo);
+				ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+
+				/* right bottom front corner */
+				if (bnd.fro==4) {
+					iyo = mod.ioZy-bnd.ntap;
+					iye = mod.ioZy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+
+
+
+				/* right bottom back corner */
+				if (bnd.bac==4) {
+					iyo = mod.ieZy;
+					iye = mod.ieZy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+			
+			ixo = mod.ioZx;
+			ixe = mod.ieZx;
+			iyo = mod.ioZy;
+			iye = mod.ieZy;
+			izo = mod.ieZz;
+			ize = mod.ieZz+bnd.ntap;
+
+			/* left bottom corner */
+			if (bnd.lef==4) {
+				ixo = mod.ioZx-bnd.ntap;
+				ixe = mod.ioZx;
+				ibz = (izo);
+				ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+//#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						#pragma ivdep
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+				/* left bottom front corner */
+				if (bnd.fro==4) {
+					iyo = mod.ioZy-bnd.ntap;
+					iye = mod.ioZy;
+					iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+						for (iy=iyo; iy<iye; iy++) {
+							#pragma ivdep
+							for (iz=izo; iz<ize; iz++) {
+								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+				/* left bottom back corner */
+				if (bnd.bac==4) {
+					iyo = mod.ieZy;
+					iye = mod.ieZy+bnd.ntap;
+					iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+					for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+						for (iy=iyo; iy<iye; iy++) {
+							for (iz=izo; iz<ize; iz++) {
+								vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz]   *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+							}
+						}
+					}
+				}
+			}
+			/* front bottom corner */
+			if (bnd.fro==4) {
+				ixo = mod.ioZx;
+				ixe = mod.ieZx;
+				iyo = mod.ioZy-bnd.ntap;
+				iye = mod.ioZy;
+				ibz = (izo);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+			/* Back bottom corner */
+			if (bnd.bac==4) {
+				iyo = mod.ieZy;
+				iye = mod.ieZy+bnd.ntap;
+				ibz = (izo);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+						}
+					}
+				}
+			}
+			
+		}
+	}
+
+
+	
+	/*********/
+	/* Left  */
+	/*********/
+	if (bnd.lef==4) {
+		
+		if (mod.ischeme <= 2) { /* Acoustic scheme */
+			
+			/* Vx field */
+			/* left mid mid vx */
+			ixo = mod.ioXx-bnd.ntap;
+			ixe = mod.ioXx;
+			iyo = mod.ioXy;
+			iye = mod.ieXy;
+			izo = mod.ioXz;
+			ize = mod.ieXz;
+			
+			ib = (bnd.ntap+ixo-1);
+
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+									c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+									c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+
+						vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[ib-ix]; 
+					}
+				}
+			}
+			/* left mid front corner vx */
+			if (bnd.fro==4) {
+				iyo = mod.ioXy-bnd.ntap;
+				iye = mod.ioXy;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+		
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)]; 
+						}
+					}
+				}
+			}
+			/* left mid back corner vx */
+			if (bnd.bac==4) {
+				iyo = mod.ieXy;
+				iye = mod.ieXy+bnd.ntap;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+		
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+						}
+					}
+				}
+			}
+
+			/* Vy field */
+			/* left mid mid vy */
+			ixo = mod.ioYx-bnd.ntap;
+			ixe = mod.ioYx;
+			iyo = mod.ioYy;
+			iye = mod.ieYy;
+			izo = mod.ioYz;
+			ize = mod.ieYz;
+			
+			ib = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+									c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+									c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+						
+						vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[ib-ix];
+					}
+				}
+			}
+			/* left mid front corner vy */
+			if (bnd.fro==4) {
+				iyo = mod.ioYy-bnd.ntap;
+				iye = mod.ioYy;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+		
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+						}
+					}
+				}
+			}
+			/* left mid back corner vy */
+			if (bnd.bac==4) {
+				iyo = mod.ieYy;
+				iye = mod.ieYy+bnd.ntap;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
+										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
+		
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+						}
+					}
+				}
+			}
+			
+			/* Vz field */
+			/* left mid mid vz */
+			ixo = mod.ioZx-bnd.ntap;
+			ixe = mod.ioZx;
+			iyo = mod.ioZy;
+			iye = mod.ieZy;
+			izo = mod.ioZz;
+			ize = mod.ieZz;
+
+			ib = (bnd.ntap+ixo-1);
+
+#pragma omp for private (ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+									c1*(tzz[iy*n1*n2+ix*n1+iz]   - tzz[iy*n1*n2+ix*n1+iz-1]) +
+									c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+						
+						vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[ib-ix];
+					}
+				}
+			}
+			/* left mid front corner vz*/
+			if (bnd.fro==4) {
+				iyo = mod.ioZy-bnd.ntap;
+				iye = mod.ioZy;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
+										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+		
+							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+						}
+					}
+				}
+			}
+			/* left mid back corner vz */
+			if (bnd.bac==4) {
+				iyo = mod.ieZy;
+				iye = mod.ieZy+bnd.ntap;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
+										c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
+										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+		
+							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+						}
+					}
+				}
+			}
+
+		}
+
+		else { /* Elastic scheme */
+			
+			/* Vx field */
+			/* left mid mid vx */
+			ixo = mod.ioXx-bnd.ntap;
+			ixe = mod.ioXx;
+			iyo = mod.ioXy;
+			iye = mod.ieXy;
+			izo = mod.ioXz;
+			ize = mod.ieXz;
+			
+			ib = (bnd.ntap+ixo-1);
+
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+
+						vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[ib-ix]; 
+					}
+				}
+			}
+			/* left mid front corner vx */
+			if (bnd.fro==4) {
+				iyo = mod.ioXy-bnd.ntap;
+				iye = mod.ioXy;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+		
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)]; 
+						}
+					}
+				}
+			}
+			/* left mid back corner vx */
+			if (bnd.bac==4) {
+				iyo = mod.ieXy;
+				iye = mod.ieXy+bnd.ntap;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+		
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+						}
+					}
+				}
+			}
+
+			/* Vy field */
+			/* left mid mid vy */
+			ixo = mod.ioYx-bnd.ntap;
+			ixe = mod.ioYx;
+			iyo = mod.ioYy;
+			iye = mod.ieYy;
+			izo = mod.ioYz;
+			ize = mod.ieYz;
+			
+			ib = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+						
+						vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[ib-ix];
+					}
+				}
+			}
+			/* left mid front corner vy */
+			if (bnd.fro==4) {
+				iyo = mod.ioYy-bnd.ntap;
+				iye = mod.ioYy;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+		
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+						}
+					}
+				}
+			}
+			/* left mid back corner vy */
+			if (bnd.bac==4) {
+				iyo = mod.ieYy;
+				iye = mod.ieYy+bnd.ntap;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+		
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+						}
+					}
+				}
+			}
+			
+			/* Vz field */
+			/* left mid mid vz */
+			ixo = mod.ioZx-bnd.ntap;
+			ixe = mod.ioZx;
+			iyo = mod.ioZy;
+			iye = mod.ieZy;
+			izo = mod.ioZz;
+			ize = mod.ieZz;
+
+			ib = (bnd.ntap+ixo-1);
+
+#pragma omp for private (ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz] *= bnd.tapz[ib-ix];
+					}
+				}
+			}
+			/* left mid front corner vz*/
+			if (bnd.fro==4) {
+				iyo = mod.ioZy-bnd.ntap;
+				iye = mod.ioZy;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+						}
+					}
+				}
+			}
+			/* left mid back corner vz */
+			if (bnd.bac==4) {
+				iyo = mod.ieZy;
+				iye = mod.ieZy+bnd.ntap;
+				ibx = (bnd.ntap+ixo-1);
+				iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+				for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+					for (iy=iyo; iy<iye; iy++) {
+						for (iz=izo; iz<ize; iz++) {
+							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+						}
+					}
+				}
+			}
+
+		}
+		
+	}
 
-	
 	/*********/
-	/* Left  */
+	/* Right */
 	/*********/
-	if (bnd.lef==4) {
+	if (bnd.rig==4) {
 		
 		if (mod.ischeme <= 2) { /* Acoustic scheme */
 			
 			/* Vx field */
-			/* left mid mid vx */
-			ixo = mod.ioXx-bnd.ntap;
-			ixe = mod.ioXx;
+			/* right mid mid vx */
+			ixo = mod.ieXx;
+			ixe = mod.ieXx+bnd.ntap;
 			iyo = mod.ioXy;
 			iye = mod.ieXy;
 			izo = mod.ioXz;
 			ize = mod.ieXz;
-			
-			ib = (bnd.ntap+ixo-1);
+		
+			ib = (ixo);
 
 #pragma omp for private(ix,iy,iz)
 			for (ix=ixo; ix<ixe; ix++) {
@@ -1433,16 +3258,17 @@ MID 	left 	mid 	mid
 									c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
 									c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
 
-						vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[ib-ix]; 
+						vx[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ix-ib]; 
 					}
 				}
 			}
-			/* left mid front corner vx */
+			/* right mid front corner vx */
 			if (bnd.fro==4) {
 				iyo = mod.ioXy-bnd.ntap;
 				iye = mod.ioXy;
-				ibx = (bnd.ntap+ixo-1);
+				ibx = (ixo);
 				iby = (bnd.ntap+iyo-1);
+
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
 #pragma ivdep
@@ -1452,17 +3278,18 @@ MID 	left 	mid 	mid
 										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
 										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
 		
-							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)]; 
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)]; 
 						}
 					}
 				}
 			}
-			/* left mid back corner vx */
+			/* right mid back corner vx*/
 			if (bnd.bac==4) {
 				iyo = mod.ieXy;
 				iye = mod.ieXy+bnd.ntap;
-				ibx = (bnd.ntap+ixo-1);
+				ibx = (ixo);
 				iby = (iyo);
+
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
 #pragma ivdep
@@ -1472,22 +3299,23 @@ MID 	left 	mid 	mid
 										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
 										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
 		
-							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
 						}
 					}
 				}
 			}
 
 			/* Vy field */
-			/* left mid mid vy */
-			ixo = mod.ioYx-bnd.ntap;
-			ixe = mod.ioYx;
+			/* right mid mid vy */
+			ixo = mod.ieYx;
+			ixe = mod.ieYx+bnd.ntap;
 			iyo = mod.ioYy;
 			iye = mod.ieYy;
 			izo = mod.ioYz;
 			ize = mod.ieYz;
-			
-			ib = (bnd.ntap+ixo-1);
+		
+			ib = (ixo);
+
 #pragma omp for private(ix,iy,iz)
 			for (ix=ixo; ix<ixe; ix++) {
 #pragma ivdep
@@ -1496,17 +3324,18 @@ MID 	left 	mid 	mid
 						vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
 									c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
 									c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
-						
-						vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[ib-ix];
+		
+						vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[ix-ib];
 					}
 				}
 			}
-			/* left mid front corner vy */
+			/* right mid front corner vy */
 			if (bnd.fro==4) {
 				iyo = mod.ioYy-bnd.ntap;
 				iye = mod.ioYy;
-				ibx = (bnd.ntap+ixo-1);
+				ibx = (ixo);
 				iby = (bnd.ntap+iyo-1);
+
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
 #pragma ivdep
@@ -1516,17 +3345,18 @@ MID 	left 	mid 	mid
 										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
 										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
 		
-							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)]; 
 						}
 					}
 				}
 			}
-			/* left mid back corner vy */
+			/* right mid back corner vy */
 			if (bnd.bac==4) {
 				iyo = mod.ieYy;
 				iye = mod.ieYy+bnd.ntap;
-				ibx = (bnd.ntap+ixo-1);
+				ibx = (ixo);
 				iby = (iyo);
+
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
 #pragma ivdep
@@ -1536,24 +3366,25 @@ MID 	left 	mid 	mid
 										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
 										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz]));
 		
-							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
 						}
 					}
 				}
 			}
-			
+
 			/* Vz field */
-			/* left mid mid vz */
-			ixo = mod.ioZx-bnd.ntap;
-			ixe = mod.ioZx;
+			/* right mid mid vz */
+
+			ixo = mod.ieZx;
+			ixe = mod.ieZx+bnd.ntap;
 			iyo = mod.ioZy;
 			iye = mod.ieZy;
 			izo = mod.ioZz;
 			ize = mod.ieZz;
+			
+			ib = (ixo);
 
-			ib = (bnd.ntap+ixo-1);
-
-#pragma omp for private (ix,iy,iz)
+#pragma omp for private (ix,iy,iz) 
 			for (ix=ixo; ix<ixe; ix++) {
 #pragma ivdep
 				for (iy=iyo; iy<iye; iy++) {
@@ -1561,17 +3392,18 @@ MID 	left 	mid 	mid
 						vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
 									c1*(tzz[iy*n1*n2+ix*n1+iz]   - tzz[iy*n1*n2+ix*n1+iz-1]) +
 									c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
-						
-						vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[ib-ix];
+		
+						vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[ix-ib]; 
 					}
 				}
 			}
-			/* left mid front corner vz*/
+			/* right mid front corner vz */
 			if (bnd.fro==4) {
 				iyo = mod.ioZy-bnd.ntap;
 				iye = mod.ioZy;
-				ibx = (bnd.ntap+ixo-1);
+				ibx = (ixo);
 				iby = (bnd.ntap+iyo-1);
+
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
 #pragma ivdep
@@ -1581,17 +3413,18 @@ MID 	left 	mid 	mid
 										c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
 										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
 		
-							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)];
 						}
 					}
 				}
 			}
-			/* left mid back corner vz */
+			/* right mid back corner vz */
 			if (bnd.bac==4) {
 				iyo = mod.ieZy;
 				iye = mod.ieZy+bnd.ntap;
-				ibx = (bnd.ntap+ixo-1);
+				ibx = (ixo);
 				iby = (iyo);
+
 #pragma omp for private(ix,iy,iz)
 				for (ix=ixo; ix<ixe; ix++) {
 #pragma ivdep
@@ -1601,22 +3434,15 @@ MID 	left 	mid 	mid
 										c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
 										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
 		
-							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+							vz[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
 						}
 					}
 				}
 			}
-
-		}
 		
-	}
+		}
 
-	/*********/
-	/* Right */
-	/*********/
-	if (bnd.rig==4) {
-		
-		if (mod.ischeme <= 2) { /* Acoustic scheme */
+		if (mod.ischeme <= 2) { /* Elastic scheme */
 			
 			/* Vx field */
 			/* right mid mid vx */
@@ -1635,8 +3461,12 @@ MID 	left 	mid 	mid
 				for (iy=iyo; iy<iye; iy++) {
 					for (iz=izo; iz<ize; iz++) {
 						vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-									c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-									c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 
 						vx[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ix-ib]; 
 					}
@@ -1655,8 +3485,12 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 		
 							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)]; 
 						}
@@ -1676,8 +3510,12 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
-										c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz]));
+										c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+											txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+										c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+											txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+											txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
 		
 							vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
 						}
@@ -1702,11 +3540,14 @@ MID 	left 	mid 	mid
 				for (iy=iyo; iy<iye; iy++) {
 					for (iz=izo; iz<ize; iz++) {
 						vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-									c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-									c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 						vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapx[ix-ib];
-						//vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapy[ix-ib]; 
 					}
 				}
 			}
@@ -1723,8 +3564,12 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)]; 
 						}
@@ -1744,8 +3589,12 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	   - tzz[(iy-1)*n1*n2+ix*n1+iz]) +
-										c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+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])  +
+										c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+											tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+											txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
 		
 							vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
 						}
@@ -1771,11 +3620,14 @@ MID 	left 	mid 	mid
 				for (iy=iyo; iy<iye; iy++) {
 					for (iz=izo; iz<ize; iz++) {
 						vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-									c1*(tzz[iy*n1*n2+ix*n1+iz]   - tzz[iy*n1*n2+ix*n1+iz-1]) +
-									c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+										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*n1*n2+ix*n1+iz] *= bnd.tapz[ix-ib]; 
-						//vz[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ix-ib]; 
 					}
 				}
 			}
@@ -1792,8 +3644,12 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
-										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+										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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)];
 						}
@@ -1813,8 +3669,12 @@ MID 	left 	mid 	mid
 					for (iy=iyo; iy<iye; iy++) {
 						for (iz=izo; iz<ize; iz++) {
 							vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*(
-										c1*(tzz[iy*n1*n2+ix*n1+iz]	 - tzz[iy*n1*n2+ix*n1+iz-1]) +
-										c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+										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*n1*n2+ix*n1+iz]   *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
 						}
@@ -1905,6 +3765,93 @@ MID 	left 	mid 	mid
 				}
 			}
 		}
+
+		else { /* Elastic scheme */
+			
+			/* Vx field */
+			/* mid mid front vx */
+			ixo = mod.ioXx;
+			ixe = mod.ieXx;
+			iyo = mod.ioXy-bnd.ntap;
+			iye = mod.ioXy;
+			izo = mod.ioXz;
+			ize = mod.ieXz;
+		
+			iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+									c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+										txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+									c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+										txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+		
+						vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapy[iby-iy]; 
+					}
+				}
+			}
+
+			/* Vy field */
+			/* mid mid front vy */
+			ixo = mod.ioYx;
+			ixe = mod.ieYx;
+			iyo = mod.ioYy-bnd.ntap;
+			iye = mod.ioYy;
+			izo = mod.ioYz;
+			ize = mod.ieYz;
+		
+			iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+									c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+										tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+										txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+		
+						vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapy[iby-iy]; 
+					}
+				}
+			}
+
+			/* Vz field */
+			/* mid mid front vz */
+			ixo = mod.ioZx;
+			ixe = mod.ieZx;
+			iyo = mod.ioZy-bnd.ntap;
+			iye = mod.ioZy;
+			izo = mod.ioZz;
+			ize = mod.ieZz;
+			
+			iby = (bnd.ntap+iyo-1);
+#pragma omp for private (ix,iy,iz) 
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz] *= bnd.tapy[iby-iy]; 
+					}
+				}
+			}
+		}
 		
 	}
 
@@ -1987,6 +3934,93 @@ MID 	left 	mid 	mid
 				}
 			}
 		}
+
+		else { /* Elastic scheme */
+			
+			/* Vx field */
+			/* mid mid back vx */
+			ixo = mod.ioXx;
+			ixe = mod.ieXx;
+			iyo = mod.ieXy;
+			iye = mod.ieXy+bnd.ntap;
+			izo = mod.ioXz;
+			ize = mod.ieXz;
+		
+			iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*(
+									c1*(txx[iy*n2*n1+ix*n1+iz]     - txx[iy*n2*n1+(ix-1)*n1+iz] +
+										txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+ix*n1+iz+1]   - txz[iy*n2*n1+ix*n1+iz])    +
+									c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+										txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+										txz[iy*n2*n1+ix*n1+iz+2]   - txz[iy*n2*n1+ix*n1+iz-1])  );
+		
+						vx[iy*n1*n2+ix*n1+iz]   *= bnd.tapy[iy-iby];
+					}
+				}
+			}
+
+			/* Vy field */
+			/* mid mid back vy */
+			ixo = mod.ioYx;
+			ixe = mod.ieYx;
+			iyo = mod.ieYy;
+			iye = mod.ieYy+bnd.ntap;
+			izo = mod.ioYz;
+			ize = mod.ieYz;
+		
+			iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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])  +
+									c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+										tyz[iy*n2*n1+ix*n1+iz+2]   - tyz[iy*n2*n1+ix*n1+iz-1] +
+										txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz])  );
+		
+						vy[iy*n1*n2+ix*n1+iz]   *= bnd.tapy[iy-iby];
+					}
+				}
+			}
+
+			/* Vz field */
+			/* mid mid back vz */
+			ixo = mod.ioZx;
+			ixe = mod.ieZx;
+			iyo = mod.ieZy;
+			iye = mod.ieZy+bnd.ntap;
+			izo = mod.ioZz;
+			ize = mod.ieZz;
+			
+			iby = (iyo);
+#pragma omp for private (ix,iy,iz) 
+			for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
+				for (iy=iyo; iy<iye; iy++) {
+					for (iz=izo; iz<ize; iz++) {
+						vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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*n1*n2+ix*n1+iz] *= bnd.tapy[iy-iby];
+					}
+				}
+			}
+		}
 		
 	}
 
@@ -1995,6 +4029,7 @@ MID 	left 	mid 	mid
 {
 		if (allocated) {
             free(Vxpml);
+			free(Vypml);
         	free(Vzpml);
         	free(sigmu);
         	free(RA);
@@ -2041,7 +4076,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa
     dy  = mod.dy;
     dt  = mod.dt;
     fac = dt/dx;
-    if ( (bnd.top==2) || (bnd.bot==2) || (bnd.lef==2) || (bnd.rig==2) ) pml=1;
+    if ( (bnd.top==2) || (bnd.bot==2) || (bnd.lef==2) || (bnd.rig==2) || (bnd.fro==2) || (bnd.bac==2) ) pml=1;
 	else pml=0;
 
 /************************************************************/
diff --git a/fdelmodc3D/fdelmodc3D b/fdelmodc3D/fdelmodc3D
index cdfa636fc5c87c6ed8961e34db637db53fb08adb..295f92c603f7e8aa6e19dda3809d39ad9687793b 100755
Binary files a/fdelmodc3D/fdelmodc3D and b/fdelmodc3D/fdelmodc3D differ