diff --git a/fdelmodc/applySource.c b/fdelmodc/applySource.c
index d29016e9b0e55cbe876f66c4fd246ad439614a49..ea1a6135a2799208f2dac312f8efe8a960de8b7b 100644
--- a/fdelmodc/applySource.c
+++ b/fdelmodc/applySource.c
@@ -97,7 +97,8 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
             ix = ixsrc + ibndx + is0 + isrc;
             iz = izsrc + ibndz;
 		}
-		time = itime*dt - src.tbeg[isrc];
+        /* mod.t0 negative time correction for dipping plane waves */
+		time = itime*dt - src.tbeg[isrc] + mod.t0;
 		id1 = floor(time/dt);
 		id2 = id1+1;
         
diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index 76229d47e7f019f30f47035cf2ca6ed8e79c0a56..694f71e2293e1d90247342af188aebe90726fd9d 100644
--- a/fdelmodc/fdelmodc.c
+++ b/fdelmodc/fdelmodc.c
@@ -610,7 +610,11 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
 				writeToFile = ! ( (((it-rec.delay)/rec.skipdt)+1)%rec.nt );
 				itwritten   = fileno*(rec.nt)*rec.skipdt;
                 /* Note that time step it=0 (t=0 for t**-fields t=-1/2 dt for v*-field) is not recorded */
-				isam        = (it-rec.delay-itwritten)/rec.skipdt+1;
+				//isam        = (it-rec.delay-itwritten)/rec.skipdt+1;
+				/* negative time correction with mod.t0 for dipping plane waves modeling */
+				isam        = (it-rec.delay-itwritten+NINT(mod.t0/mod.dt))/rec.skipdt+1;
+				if (isam < 0) isam = rec.nt+isam;
+
 				/* store time at receiver positions */
 				getRecTimes(mod, rec, bnd, it, isam, vx, vz, tzz, txx, txz, 
 					l2m, rox, roz, 
@@ -682,6 +686,8 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
 				}
 			}
 		}
+		/* negative time correction with mod.t0 for dipping plane waves modeling */
+        isam += NINT(-mod.t0/mod.dt)/rec.skipdt;
 		writeRec(rec, mod, bnd, wav, ixsrc, izsrc, isam+1, ishot, fileno,
 			rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz, 
 			rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
diff --git a/fdelmodc/fdelmodc.h b/fdelmodc/fdelmodc.h
index f68d00cfcb3fe8d3b0aff1ba4d2abb0c2969ae32..860a035a12535e5974b9a4ebf24f272455ab0fde 100644
--- a/fdelmodc/fdelmodc.h
+++ b/fdelmodc/fdelmodc.h
@@ -73,6 +73,7 @@ typedef struct _modelPar { /* Model Parameters */
 	int nt;
 	float z0;
 	float x0;
+	float t0;
 	/* medium max/min values */
 	float cp_min;
 	float cp_max;
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index 0894e7c7338f11b5a4fa779e5d7bda2cb3999583..988db0e788568924c2f2d1918fa3d88937eeacdc 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -174,6 +174,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
 	if (!getparfloat("rec_delay",&rdelay)) rdelay=0.0;
 	rec->delay=NINT(rdelay/mod->dt);
 	mod->nt = NINT(mod->tmod/mod->dt);
+    mod->t0 = 0.0;
 	dt = mod->dt;
 
 	if (!getparint("src_type",&src->type)) src->type=1;
@@ -878,6 +879,10 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
 			vwarn("number of gridpoints in X. Plane wave will be clipped to the edges of the model");
 			nsrc = mod->nx;
 		}
+	    src_ix0=MAX(0,NINT((xsrc-sub_x0)/dx));
+	    src_ix0=MIN(src_ix0,nx);
+	    src_iz0=MAX(0,NINT((zsrc-sub_z0)/dz));
+	    src_iz0=MIN(src_iz0,nz);
 
 	/* for a source defined on mutliple gridpoint calculate p delay factor */
 
@@ -887,6 +892,11 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
 		src->tend = (float *)malloc(nsrc*sizeof(float));
 		grad2rad = 17.453292e-3;
 		p = sin(src_angle*grad2rad)/src_velo;
+		for (is=0; is<nsrc; is++) {
+			src->tbeg[is] = (is-src_ix0)*dx*p;
+		}
+		mod->t0 = MIN(-src_ix0*dx*p,(nsrc-src_ix0)*dx*p);
+/*
 		if (p < 0.0) {
 			for (is=0; is<nsrc; is++) {
 				src->tbeg[is] = fabsf((nsrc-is-1)*dx*p);
@@ -899,11 +909,13 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
 				src->tbeg[is] = is*dx*p;
 			}
 		}
+*/
 		for (is=0; is<nsrc; is++) {
 			src->tend[is] = src->tbeg[is] + (wav->nt-1)*wav->dt;
 		}
 		
 		is0 = -1*floor((nsrc-1)/2);
+		is0 = -src_ix0;
 		for (is=0; is<nsrc; is++) {
 			src->x[is] = is0 + is;
 			src->z[is] = 0;
@@ -986,7 +998,10 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
 			vmess("Areal source array is defined with %d sources.",nsrc);
 /*			vmess("Memory requirement for sources = %.2f MB.",sizeof(float)*(wav->nx*(wav->nt/(1024.0*1024.0))));*/
 			vmess("Memory requirement for sources = %.2f MB.",sizeof(float)*(nsamp/(1024.0*1024.0)));
-			if (src->plane) vmess("Computed p-value = %f.",p);
+			if (src->plane) {
+				vmess("Computed p-value = %f.",p);
+        		vmess("Maximum negative time delay is %f\n", mod->t0);
+			}
 		}
 		if (src->random) {
 		vmess("Sources are placed at random locations in domain: ");