diff --git a/fdelmodc3D/fdelmodc3D.h b/fdelmodc3D/fdelmodc3D.h
index 76b05012b1037659747da57738ae0301135defee..158899df6305898453f54fc4344e21bd645b7fe8 100644
--- a/fdelmodc3D/fdelmodc3D.h
+++ b/fdelmodc3D/fdelmodc3D.h
@@ -201,6 +201,7 @@ typedef struct _boundPar { /* Boundary Parameters */
 	float *tapy;
 	float *tapx;
 	float *tapxz;
+	float *tapxyz;
 	long cfree;
 	long ntap;
 	long *surface;
diff --git a/fdelmodc3D/getParameters3D.c b/fdelmodc3D/getParameters3D.c
index 701269b268ef47d6b2301be49e010e4cecfb6f8e..df33bf7fd992374c828f96af4ad166ac0fd6c483 100644
--- a/fdelmodc3D/getParameters3D.c
+++ b/fdelmodc3D/getParameters3D.c
@@ -52,7 +52,7 @@ long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar
 	float *xsrca, *ysrca, *zsrca, rrcv;
 	float rsrc, oxsrc, oysrc, ozsrc, dphisrc, ncsrc;
 	size_t nsamp;
-	long i, j;
+	long i, j, l;
 	long cfree;
 	long tapleft,tapright,taptop,tapbottom,tapfront, tapback;
 	long nxsrc, nysrc, nzsrc;
@@ -347,10 +347,11 @@ long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar
 	bnd->npml=bnd->ntap;
 
 	if (bnd->ntap) {
-		bnd->tapx  = (float *)malloc(bnd->ntap*sizeof(float));
-        bnd->tapy  = (float *)malloc(bnd->ntap*sizeof(float));
-		bnd->tapz  = (float *)malloc(bnd->ntap*sizeof(float));
-		bnd->tapxz = (float *)malloc(bnd->ntap*bnd->ntap*sizeof(float));
+		bnd->tapx   = (float *)malloc(bnd->ntap*sizeof(float));
+        bnd->tapy   = (float *)malloc(bnd->ntap*sizeof(float));
+		bnd->tapz   = (float *)malloc(bnd->ntap*sizeof(float));
+		bnd->tapxz  = (float *)malloc(bnd->ntap*bnd->ntap*sizeof(float));
+		bnd->tapxyz = (float *)malloc(bnd->ntap*bnd->ntap*bnd->ntap*sizeof(float));
         if(!getparfloat("tapfact",&tapfact)) tapfact=0.30;
 		scl = tapfact/((float)bnd->ntap);
 		for (i=0; i<bnd->ntap; i++) {
@@ -368,6 +369,14 @@ long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar
 				bnd->tapxz[j*bnd->ntap+i] = exp(-(wfct*wfct));
 			}
 		}
+		for (j=0; j<bnd->ntap; j++) {
+			for (i=0; i<bnd->ntap; i++) {
+				for (l=0; l<bnd->ntap; i++) {
+					wfct = (scl*sqrt(i*i+j*j+l*l));
+					bnd->tapxyz[l*bnd->ntap*bnd->ntap+j*bnd->ntap+i] = exp(-(wfct*wfct));
+				}
+			}
+		}
 	}
     
     /* Vx: rox */