diff --git a/src/IFOS2D.c b/src/IFOS2D.c
index 50a03e8181fb7a75199b45b2ebfe54f8609a46be..85e008c6768939ba851c7369ca0c8f50461d95dc 100644
--- a/src/IFOS2D.c
+++ b/src/IFOS2D.c
@@ -99,6 +99,9 @@ int main(int argc, char **argv){
     int SOURCE_SHAPE_OLD=0;
     int SOURCE_SHAPE_OLD_SH=0;
     
+    /* Variables for conjungate gradient */
+    int PCG_iter_start=1;
+    
     /* Variables for L-BFGS */
     int LBFGS_NPAR=3;
     int LBFGS_iter_start=1;
@@ -222,9 +225,6 @@ int main(int argc, char **argv){
              at the receiver positions to calculate the adjoint sources and to do
              the backpropagation; look at function saveseis_glob.c to see that every
              NDT sample for the forward modeled wavefield is written to su files*/
-    lsnap=iround(TSNAP1/DT);      /* first snapshot at this timestep */
-    lsamp=NDT;
-    
     
     /* output of parameters to log-file or stdout */
     if (MYID==0) write_par(FP);
@@ -504,13 +504,7 @@ int main(int argc, char **argv){
         cjm =  vector(1,L);
     }
     
-    /*nf=4;
-     nfstart=4;*/
-    
     NTST=20;
-    NTSTI=NTST/DTINV;
-    
-    nxny=NX*NY;
     nxnyi=(NX/IDXI)*(NY/IDYI);
     
     /* Parameters for step length calculations */
@@ -2900,17 +2894,25 @@ int main(int argc, char **argv){
                 /* ----------------------------------------------------------------------*/
                 /* ----------- Preconditioned Conjugate Gradient Method (PCG)  ----------*/
                 /* ----------------------------------------------------------------------*/
-                if((GRAD_METHOD==1)){
+                if(GRAD_METHOD==1){
+                    
+                    if( (iter-PCG_iter_start) < 2 ) {
+                        fprintf(FP,"\n\n ----- Conjugate Gradient Method -----");
+                        fprintf(FP,"\n Will not use second last gradient for conjugate gradient");
+                        if( (iter-PCG_iter_start) < 1 ) {
+                            fprintf(FP,"\n Will not use last gradient for conjugate gradient");
+                        }
+                    }
                     
-                    if (WAVETYPE==1 || WAVETYPE==3) PCG(waveconv, taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, C_vp, gradp, nfstart_jac, waveconv_u, C_vs, gradp_u, waveconv_rho, C_rho, gradp_rho,Vs_avg,F_LOW_PASS);
-                    if (WAVETYPE==2 || WAVETYPE==3) PCG_SH(taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, nfstart_jac, waveconv_u_z, C_vs, gradp_u_z, waveconv_rho_z, C_rho, gradp_rho_z,Vs_avg,F_LOW_PASS);
+                    if (WAVETYPE==1 || WAVETYPE==3) PCG(waveconv, taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, C_vp, gradp, nfstart_jac, waveconv_u, C_vs, gradp_u, waveconv_rho, C_rho, gradp_rho,Vs_avg,F_LOW_PASS,PCG_iter_start);
+                    if (WAVETYPE==2 || WAVETYPE==3) PCG_SH(taper_coeff, nsrc, srcpos, recpos, ntr_glob, iter, nfstart_jac, waveconv_u_z, C_vs, gradp_u_z, waveconv_rho_z, C_rho, gradp_rho_z,Vs_avg,F_LOW_PASS,PCG_iter_start);
                     
                 }
                 
                 /* ---------------------------------------------------------*/
                 /* ----------- Beginn Joint Inversion PSV and SH  ----------*/
                 /* ---------------------------------------------------------*/
-                if((FORWARD_ONLY==0)){
+                if(FORWARD_ONLY==0){
                     switch (WAVETYPE) {
                         case 2:
                             /* If only SH is inverted, set these gradients to the actual ones */
@@ -2950,7 +2952,7 @@ int main(int argc, char **argv){
                 /* -------------------------------------------------------------------------*/
                 /* ----------- Limited Memory - Broyden-Fletcher-Goldfarb-Shanno  ----------*/
                 /* -------------------------------------------------------------------------*/
-                if((GRAD_METHOD==2)){
+                if(GRAD_METHOD==2){
                     
                     /*---------------------*/
                     /*         TAPER       */
@@ -3903,6 +3905,9 @@ int main(int argc, char **argv){
                 /* Restart L-BFGS at next iteration */
                 LBFGS_iter_start=iter+1;
                 
+                /* Restart conjugate gradient at next iteration */
+                PCG_iter_start=iter+1;
+                
                 wolfe_SLS_failed=0;
             }
             
@@ -3933,6 +3938,10 @@ int main(int argc, char **argv){
                 
                 /* Restart L-BFGS at next iteration */
                 LBFGS_iter_start=iter+1;
+                
+                /* Restart conjugate gradient at next iteration */
+                PCG_iter_start=iter+1;
+                
                 wolfe_SLS_failed=0;
                 alpha_SL_old=1;
             }
@@ -3964,6 +3973,10 @@ int main(int argc, char **argv){
                 
                 /* Restart L-BFGS at next iteration */
                 LBFGS_iter_start=iter+1;
+                
+                /* Restart conjungate gradient at next iteration */
+                PCG_iter_start=iter+1;
+                
                 wolfe_SLS_failed=0;
                 alpha_SL_old=1;
             }
diff --git a/src/PCG.c b/src/PCG.c
index 37814942a913bf11c7256239626e1196b110c5d2..dac01f04439149c36d0760941060fc55f3776aa3 100644
--- a/src/PCG.c
+++ b/src/PCG.c
@@ -24,7 +24,7 @@
 
 #include "fd.h"
 
-void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS){
+void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start){
 	
 	extern int NX, NY, IDX, IDY, SPATFILTER, GRAD_FILTER;
 	extern int FORWARD_ONLY, SWS_TAPER_GRAD_VERT, SWS_TAPER_GRAD_HOR, SWS_TAPER_GRAD_SOURCES, SWS_TAPER_FILE;
@@ -41,8 +41,14 @@ void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int
     int use_conjugate_2=1;
 	FILE *FP3, *FP4, *FP6, *FP5 = NULL;
 	
-
-
+    /* Check if conjugate gradient can be used */
+    if( (iter-PCG_iter_start) < 2 ) {
+        use_conjugate_2=0;
+        if( (iter-PCG_iter_start) < 1 ) {
+            use_conjugate_1=0;
+        }
+    }
+    
 	/* ===================================================================================================================================================== */
 	/* ===================================================== GRADIENT ZP ================================================================================== */
 	/* ===================================================================================================================================================== */
diff --git a/src/PCG_SH.c b/src/PCG_SH.c
index 6d56035f3e80de8b1edfb46c7125305f71ae184a..a1d96589039076dbf2584b8c424e09e72a57c131 100644
--- a/src/PCG_SH.c
+++ b/src/PCG_SH.c
@@ -24,7 +24,7 @@
 
 #include "fd.h"
 
-void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS){
+void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start){
     
     extern int NX, NY, IDX, IDY, SPATFILTER, GRAD_FILTER;
     extern int FORWARD_ONLY, SWS_TAPER_GRAD_VERT, SWS_TAPER_GRAD_HOR, SWS_TAPER_GRAD_SOURCES, SWS_TAPER_FILE;
@@ -40,6 +40,14 @@ void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int
     int use_conjugate_2=1;
 
     
+    /* Check if conjugate gradient can be used */
+    if( (iter-PCG_iter_start) < 2 ) {
+        use_conjugate_2=0;
+        if( (iter-PCG_iter_start) < 1 ) {
+            use_conjugate_1=0;
+        }
+    }
+    
     /* ===================================================================================================================================================== */
     /* ===================================================== GRADIENT Zs ================================================================================== */
     /* ===================================================================================================================================================== */
diff --git a/src/fd.h b/src/fd.h
index 1c55f3277a17e76db64fb871c28dc9b73aacf763..fe474d98401a1f0aa75cdebd01e6566c909568bd 100644
--- a/src/fd.h
+++ b/src/fd.h
@@ -186,9 +186,9 @@ void  taper(float *section, int ns, float fc);
 
 void  output_source_signal(FILE *fp, float **signals, int ns, int seis_form);
 
-void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS);
+void PCG(float ** waveconv, float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, float C_vp, float ** gradp, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start);
 
-void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS);
+void PCG_SH(float ** taper_coeff, int nsrc, float ** srcpos, int ** recpos, int ntr_glob, int iter, int nfstart_jac, float ** waveconv_u, float C_vs, float ** gradp_u, float ** waveconv_rho, float C_rho, float ** gradp_rho, float Vs_avg, float F_LOW_PASS, int PCG_iter_start);
 
 void PML_pro(float * d_x, float * K_x, float * alpha_prime_x, float * a_x, float * b_x,
              float * d_x_half, float * K_x_half, float * alpha_prime_x_half, float * a_x_half, float * b_x_half,