fixed bug in checkpoint write/read

checkpoints were overwritten by zero functions.
Also PML boundaries and viscoelastic modelling now work with checkpointing.
I didnt't implement higher fd-orders in time because checkpoint files would become very large.

doc/guide_sofi3D/tex/guide_sofi3D.tex

@@ -1208,7 +1208,7 @@ CHECKPTREAD : read wavefield from checkpoint file (yes=1/no=0)\\
 CHECKPTWRITE : save wavefield to checkpoint file (yes=1/no=0)\\
 CHECKPTFILE : checkpoint file name\\
 
-On most supercomputers with a queuing system the run time of job is limited. Sometimes the allowed run time is not sufficient to finish a FD simulation. In such a case, check-pointing can be performed: the first job saves the complete elastic wavefield (CHECKPTWRITE=1) but does not! read the wavefield from a checkpoint file (CHECKPTREAD=0). The subsequent jobs read and write the wavefield to the CHECKPTFILE, i.e. CHECKPTREAD=1 and CHECKPTWRITE=1. In this manner, one simulation can be divided on different batch jobs. The resulting seismograms may be catenated using the SU-command sucat. But be aware that the checkpointing option saves the COMPLETE wavefield information of the last time step, which can produce a huge amount of data.
+On most supercomputers with a queuing system the run time of job is limited. Sometimes the allowed run time is not sufficient to finish a FD simulation. In such a case, check-pointing can be performed: the first job saves the complete elastic wavefield (CHECKPTWRITE=1) but does not! read the wavefield from a checkpoint file (CHECKPTREAD=0). The subsequent jobs read and write the wavefield to the CHECKPTFILE, i.e. CHECKPTREAD=1 and CHECKPTWRITE=1. In this manner, one simulation can be divided on different batch jobs. The resulting seismograms may be catenated using the SU-command suvcat. But be aware that the checkpointing option saves the COMPLETE wavefield information of the last time step, which can produce a huge amount of data.
 
 \subsection{''On the fly'' definition of material parameters}
 \label{model_def_func}

src/checkfd_ssg.c

@@ -301,7 +301,11 @@ void checkfd(FILE *fp, float *** prho, float *** ppi, float *** pu,
 		fprintf(fp," Check writability for checkpoint files %s... \n",xfile);
 		if (((fpcheck=fopen(xfile,xmod))==NULL) && (MYID==0)) err(" PE 0 cannot write checkpoints!");
 		else fclose(fpcheck); /* Is there any reason to remove it? */	
-	}	
+		if(MYID==0){
+			if (FDORDER_TIME > 2){warning(" Checkpoint writing for FDORDER_TIME > 2 is not implemented,\n wavefields calculated from the checkpoint will slightly differ.\n");}
+		}	
+		
+	}
 
 	fprintf(fp," Accessibility of output files from PE %d has been checked successfully.\n", MYID);
 

src/fd.h

@@ -132,8 +132,12 @@ double exchange_v(int nt, float *** vx, float *** vy, float *** vz,
 
 void read_checkpoint(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2,
 		float *** vx, float *** vy, float *** vz,
-		float *** sxx, float *** syy, float *** szz, float *** sxy,
-		float *** syz, float *** sxz);
+		float *** sxx, float *** syy, float *** szz, float *** sxy,float *** syz, float *** sxz,
+		float *** rxx, float *** ryy,float *** rzz, float *** rxy, float *** ryz, float *** rxz,
+		float *** psi_sxx_x, float *** psi_sxy_x, float *** psi_sxz_x, float *** psi_sxy_y,
+		float *** psi_syy_y, float *** psi_syz_y, float *** psi_sxz_z, float *** psi_syz_z, float *** psi_szz_z,
+		float *** psi_vxx, float *** psi_vyx, float *** psi_vzx, float *** psi_vxy, float *** psi_vyy, float *** psi_vzy,
+		float *** psi_vxz, float *** psi_vyz, float *** psi_vzz);
 
 void info(FILE *fp);
 
@@ -218,7 +222,11 @@ void saveseis_glob(FILE *fp, float **sectiondata,
 void save_checkpoint(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2,
 		float *** vx, float *** vy, float *** vz,
 		float *** sxx, float *** syy, float *** szz, float *** sxy,
-		float *** syz, float *** sxz);
+		float *** rxx, float *** ryy,float *** rzz, float *** rxy, float *** ryz, float *** rxz,
+		float *** syz, float *** sxz, float *** psi_sxx_x, float *** psi_sxy_x, float *** psi_sxz_x, 
+		float *** psi_sxy_y, float *** psi_syy_y, float *** psi_syz_y, float *** psi_sxz_z, float *** psi_syz_z, float *** psi_szz_z, 
+		float *** psi_vxx, float *** psi_vyx, float *** psi_vzx, float *** psi_vxy, float *** psi_vyy, float *** psi_vzy, 
+		float *** psi_vxz, float *** psi_vyz, float *** psi_vzz);
 
 void seismo_acoustic(int lsamp, int ntr, int **recpos, float **sectionvx, float **sectionvy, 
 		float **sectionvz, float **sectiondiv, float **sectioncurl, float **sectionp,

src/read_checkpoint.c

@@ -2,18 +2,18 @@
  * Copyright (C) 2011 For the list of authors, see file AUTHORS.
  *
  * This file is part of SOFI3D.
- * 
+ *
  * SOFI3D is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, version 2.0 of the License only.
- * 
+ *
  * SOFI3D is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  * GNU General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU General Public License
- * along with SOFI3D. See file COPYING and/or 
+ * along with SOFI3D. See file COPYING and/or
   * <http://www.gnu.org/licenses/gpl-2.0.html>.
 --------------------------------------------------------------------------*/
 /* ----------------------------------------------------------------------
@@ -21,18 +21,23 @@
  ----------------------------------------------------------------------*/
 #include "fd.h"
 
-void read_checkpoint(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2,
-float *** vx, float *** vy, float *** vz,
-float *** sxx, float *** syy, float *** szz, float *** sxy,
-float *** syz, float *** sxz){
+void read_checkpoint(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float *** vx, float *** vy, float *** vz,
+                     float *** sxx, float *** syy, float *** szz, float *** sxy, float *** syz, float *** sxz,
+		     float *** rxx, float *** ryy,float *** rzz, float *** rxy, float *** ryz, float *** rxz,
+		     float *** psi_sxx_x, float *** psi_sxy_x, float *** psi_sxz_x, float *** psi_sxy_y,
+		     float *** psi_syy_y, float *** psi_syz_y, float *** psi_sxz_z, float *** psi_syz_z, float *** psi_szz_z,
+                     float *** psi_vxx, float *** psi_vyx, float *** psi_vzx, float *** psi_vxy, float *** psi_vyy, float *** psi_vzy,
+                     float *** psi_vxz, float *** psi_vyz, float *** psi_vzz) {
 
 	int i,j,k;
 	char myid[5];
 	FILE *fp;
 	char checkptfile[STRING_SIZE];
 	extern int MYID;
+	extern int ABS_TYPE, NX,NY,NZ,FW,POS[4],NPROCX,NPROCY,NPROCZ,FREE_SURF;
+	extern int L;
 	extern char  CHECKPTFILE[STRING_SIZE];
-
+	
 
 
 	sprintf(checkptfile,"%s",CHECKPTFILE);
@@ -42,24 +47,136 @@ float *** syz, float *** sxz){
 
 
 	fp=fopen(checkptfile,"rb");
-	if (fp==NULL) err("CHECKPTFILE can't be opened !");
+	if (fp==NULL) {
+		err("CHECKPTFILE can't be opened !");
+	}
+
+
+	for (j=ny1; j<=ny2; j++) {
+		for (i=nx1; i<=nx2; i++) {
+			for (k=nz1; k<=nz2; k++) {
+
+				fread(&vx[j][i][k],sizeof(float),1,fp);
+				fread(&vy[j][i][k],sizeof(float),1,fp);
+				fread(&vz[j][i][k],sizeof(float),1,fp);
+				fread(&sxx[j][i][k],sizeof(float),1,fp);
+				fread(&syy[j][i][k],sizeof(float),1,fp);
+				fread(&szz[j][i][k],sizeof(float),1,fp);
+				fread(&sxy[j][i][k],sizeof(float),1,fp);
+				fread(&syz[j][i][k],sizeof(float),1,fp);
+				fread(&sxz[j][i][k],sizeof(float),1,fp);
+			}
+		}
+	}
+
+	
+	if (L) {
+		for (j=1; j<=NY; j++) {
+			for (i=1; i<=NX; i++) {
+				for (k=1; k<=NZ; k++) {
+					fread(&rxx[j][i][k],sizeof(float),1,fp);
+					fread(&ryy[j][i][k],sizeof(float),1,fp);
+					fread(&rzz[j][i][k],sizeof(float),1,fp);
+					fread(&rxy[j][i][k],sizeof(float),1,fp);
+					fread(&ryz[j][i][k],sizeof(float),1,fp);
+					fread(&rxz[j][i][k],sizeof(float),1,fp);
+				}
+			}
+		}
+	}
 	
+	
+	
+	if (ABS_TYPE == 1) {
+
+		if (POS[1]==0) {
+			for (j=1; j<=NY; j++) {
+				for (i=1; i<=FW; i++) {
+					for (k=1; k<=NZ; k++) {
+						fread(&psi_sxx_x[j][i][k],sizeof(float),1,fp);
+						fread(&psi_sxy_x[j][i][k],sizeof(float),1,fp);
+						fread(&psi_sxz_x[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vxx[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vyx[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vzx[j][i][k],sizeof(float),1,fp);
+					}
+				}
+			}
+		}
+
+		if (POS[1]==NPROCX-1) {
+			for (j=1; j<=NY; j++) {
+				for (i=FW+1; i<=2*FW; i++) {
+					for (k=1; k<=NZ; k++) {
+						fread(&psi_sxx_x[j][i][k],sizeof(float),1,fp);
+						fread(&psi_sxy_x[j][i][k],sizeof(float),1,fp);
+						fread(&psi_sxz_x[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vxx[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vyx[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vzx[j][i][k],sizeof(float),1,fp);
+					}
+				}
+			}
+		}
 
-	for (j=ny1;j<=ny2;j++){
-		for (i=nx1;i<=nx2;i++){
-			for (k=nz1;k<=nz2;k++){
-
-			fread( &vx[j][i][k],sizeof(float),1,fp);
-			fread( &vy[j][i][k],sizeof(float),1,fp);
-			fread( &vz[j][i][k],sizeof(float),1,fp);
-			fread(&sxx[j][i][k],sizeof(float),1,fp);
-			fread(&syy[j][i][k],sizeof(float),1,fp);
-			fread(&szz[j][i][k],sizeof(float),1,fp);
-			fread(&sxy[j][i][k],sizeof(float),1,fp);
-			fread(&syz[j][i][k],sizeof(float),1,fp);
-			fread(&sxz[j][i][k],sizeof(float),1,fp);
+		if (POS[2]==0 && FREE_SURF==0) {
+			for (j=1; j<=FW; j++) {
+				for (i=1; i<=NX; i++) {
+					for (k=1; k<=NZ; k++) {
+						fread(&psi_sxy_y[j][i][k],sizeof(float),1,fp);
+						fread(&psi_syy_y[j][i][k],sizeof(float),1,fp);
+						fread(&psi_syz_y[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vxy[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vyy[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vzy[j][i][k],sizeof(float),1,fp);
+					}
+				}
 			}
 		}
+		if (POS[2]==NPROCY-1) {
+			for (j=FW+1; j<=2*FW; j++) {
+				for (i=1; i<=NX; i++) {
+					for (k=1; k<=NZ; k++) {
+						fread(&psi_sxy_y[j][i][k],sizeof(float),1,fp);
+						fread(&psi_syy_y[j][i][k],sizeof(float),1,fp);
+						fread(&psi_syz_y[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vxy[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vyy[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vzy[j][i][k],sizeof(float),1,fp);
+					}
+				}
+			}
+		}
+
+		if (POS[3]==0) {
+			for (j=1; j<=NY; j++) {
+				for (i=1; i<=NX; i++) {
+					for (k=1; k<=FW; k++) {
+						fread(&psi_sxz_z[j][i][k],sizeof(float),1,fp);
+						fread(&psi_syz_z[j][i][k],sizeof(float),1,fp);
+						fread(&psi_szz_z[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vxz[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vyz[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vzz[j][i][k],sizeof(float),1,fp);
+					}
+				}
+			}
+		}
+		if (POS[3]==NPROCZ-1) {
+			for (j=1; j<=NY; j++) {
+				for (i=1; i<=NX; i++) {
+					for (k=FW+1; k<=2*FW; k++) {
+						fread(&psi_sxz_z[j][i][k],sizeof(float),1,fp);
+						fread(&psi_syz_z[j][i][k],sizeof(float),1,fp);
+						fread(&psi_szz_z[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vxz[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vyz[j][i][k],sizeof(float),1,fp);
+						fread(&psi_vzz[j][i][k],sizeof(float),1,fp);
+					}
+				}
+			}
+		}
+
 	}
 
 	fclose(fp);

src/save_checkpoint.c

@@ -2,18 +2,18 @@
  * Copyright (C) 2011 For the list of authors, see file AUTHORS.
  *
  * This file is part of SOFI3D.
- * 
+ *
  * SOFI3D is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, version 2.0 of the License only.
- * 
+ *
  * SOFI3D is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  * GNU General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU General Public License
- * along with SOFI3D. See file COPYING and/or 
+ * along with SOFI3D. See file COPYING and/or
   * <http://www.gnu.org/licenses/gpl-2.0.html>.
 --------------------------------------------------------------------------*/
 /*----------------------------------------------------------------------
@@ -22,16 +22,21 @@
  ---------------------------------------------------------------------- */
 #include "fd.h"
 
-void save_checkpoint(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2,
-float *** vx, float *** vy, float *** vz,
-float *** sxx, float *** syy, float *** szz, float *** sxy,
-float *** syz, float *** sxz){
+void save_checkpoint(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float *** vx, float *** vy, float *** vz,
+                     float *** sxx, float *** syy, float *** szz, float *** sxy,
+                     float *** syz, float *** sxz, float *** rxx, float *** ryy,float *** rzz, float *** rxy, float *** ryz, float *** rxz,
+		     float *** psi_sxx_x, float *** psi_sxy_x, float *** psi_sxz_x, float *** psi_sxy_y,
+		     float *** psi_syy_y, float *** psi_syz_y, float *** psi_sxz_z, float *** psi_syz_z, float *** psi_szz_z,
+                     float *** psi_vxx, float *** psi_vyx, float *** psi_vzx, float *** psi_vxy, float *** psi_vyy, float *** psi_vzy,
+                     float *** psi_vxz, float *** psi_vyz, float *** psi_vzz) {
 
 	int i,j, k;
 	char myid[5];
 	FILE *fp;
 	char checkptfile[STRING_SIZE];
 	extern int MYID;
+	extern int ABS_TYPE, NX,NY,NZ,FW,POS[4],NPROCX,NPROCY,NPROCZ,FREE_SURF;
+	extern int L;
 	extern char  CHECKPTFILE[STRING_SIZE];
 
 
@@ -43,26 +48,149 @@ float *** syz, float *** sxz){
 
 
 	fp=fopen(checkptfile,"wb");
-	if (fp==NULL) err("CHECKPTFILE can't be opened !");
-	
-
-	for (j=ny1;j<=ny2;j++){
-		for (i=nx1;i<=nx2;i++){
-			for (k=nz1;k<=nz2;k++){
-
-			fwrite( &vx[j][i][k],sizeof(float),1,fp);
-			fwrite( &vy[j][i][k],sizeof(float),1,fp);
-			fwrite( &vz[j][i][k],sizeof(float),1,fp);
-			fwrite(&sxx[j][i][k],sizeof(float),1,fp);
-			fwrite(&syy[j][i][k],sizeof(float),1,fp);
-			fwrite(&szz[j][i][k],sizeof(float),1,fp);
-			fwrite(&sxy[j][i][k],sizeof(float),1,fp);
-			fwrite(&syz[j][i][k],sizeof(float),1,fp);
-			fwrite(&sxz[j][i][k],sizeof(float),1,fp);
+	if (fp==NULL) {
+		err("CHECKPTFILE can't be opened !");
+	}
+
+
+	for (j=ny1; j<=ny2; j++) {
+		for (i=nx1; i<=nx2; i++) {
+			for (k=nz1; k<=nz2; k++) {
+				fwrite(&vx[j][i][k],sizeof(float),1,fp);
+				fwrite(&vy[j][i][k],sizeof(float),1,fp);
+				fwrite(&vz[j][i][k],sizeof(float),1,fp);
+				fwrite(&sxx[j][i][k],sizeof(float),1,fp);
+				fwrite(&syy[j][i][k],sizeof(float),1,fp);
+				fwrite(&szz[j][i][k],sizeof(float),1,fp);
+				fwrite(&sxy[j][i][k],sizeof(float),1,fp);
+				fwrite(&syz[j][i][k],sizeof(float),1,fp);
+				fwrite(&sxz[j][i][k],sizeof(float),1,fp);
+			}
+		}
+	}
+
+	if (L) {
+		for (j=1; j<=NY; j++) {
+			for (i=1; i<=NX; i++) {
+				for (k=1; k<=NZ; k++) {
+					fwrite(&rxx[j][i][k],sizeof(float),1,fp);
+					fwrite(&ryy[j][i][k],sizeof(float),1,fp);
+					fwrite(&rzz[j][i][k],sizeof(float),1,fp);
+					fwrite(&rxy[j][i][k],sizeof(float),1,fp);
+					fwrite(&ryz[j][i][k],sizeof(float),1,fp);
+					fwrite(&rxz[j][i][k],sizeof(float),1,fp);
+				}
+			}
+		}
+	}
+
+	if (ABS_TYPE == 1) {
+
+		if (POS[1]==0) {
+			for (j=1; j<=NY; j++) {
+				for (i=1; i<=FW; i++) {
+					for (k=1; k<=NZ; k++) {
+						fwrite(&psi_sxx_x[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_sxy_x[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_sxz_x[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vxx[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vyx[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vzx[j][i][k],sizeof(float),1,fp);
+					}
+				}
+			}
+		}
+
+		if (POS[1]==NPROCX-1) {
+			for (j=1; j<=NY; j++) {
+				for (i=FW+1; i<=2*FW; i++) {
+					for (k=1; k<=NZ; k++) {
+						fwrite(&psi_sxx_x[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_sxy_x[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_sxz_x[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vxx[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vyx[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vzx[j][i][k],sizeof(float),1,fp);
+					}
+				}
+			}
+		}
+
+		if (POS[2]==0 && FREE_SURF==0) {
+			for (j=1; j<=FW; j++) {
+				for (i=1; i<=NX; i++) {
+					for (k=1; k<=NZ; k++) {
+						fwrite(&psi_sxy_y[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_syy_y[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_syz_y[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vxy[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vyy[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vzy[j][i][k],sizeof(float),1,fp);
+					}
+				}
+			}
+		}
+		if (POS[2]==NPROCY-1) {
+			for (j=FW+1; j<=2*FW; j++) {
+				for (i=1; i<=NX; i++) {
+					for (k=1; k<=NZ; k++) {
+						fwrite(&psi_sxy_y[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_syy_y[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_syz_y[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vxy[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vyy[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vzy[j][i][k],sizeof(float),1,fp);
+					}
+				}
+			}
+		}
+
+		if (POS[3]==0) {
+			for (j=1; j<=NY; j++) {
+				for (i=1; i<=NX; i++) {
+					for (k=1; k<=FW; k++) {
+						fwrite(&psi_sxz_z[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_syz_z[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_szz_z[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vxz[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vyz[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vzz[j][i][k],sizeof(float),1,fp);
+					}
+				}
+			}
+		}
+		if (POS[3]==NPROCZ-1) {
+			for (j=1; j<=NY; j++) {
+				for (i=1; i<=NX; i++) {
+					for (k=FW+1; k<=2*FW; k++) {
+						fwrite(&psi_sxz_z[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_syz_z[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_szz_z[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vxz[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vyz[j][i][k],sizeof(float),1,fp);
+						fwrite(&psi_vzz[j][i][k],sizeof(float),1,fp);
+					}
+				}
 			}
 		}
+
 	}
 
 	fclose(fp);
 
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

src/sofi3D.c

@@ -740,7 +740,9 @@ int main(int argc, char **argv){
 			fprintf(FP," Reading wavefield from check-point file %s \n",CHECKPTFILE);
 		}
 
-		read_checkpoint(-1, NX+2, -1, NY+2, -1, NZ+2, vx,vy,vz,sxx,syy,szz,sxy,syz,sxz);
+		read_checkpoint(-1, NX+2, -1, NY+2, -1, NZ+2, vx,vy,vz,sxx,syy,szz,sxy,syz,sxz,rxx,ryy,rzz,rxy,ryz,rxz,
+			psi_sxx_x,psi_sxy_x,psi_sxz_x,psi_sxy_y,psi_syy_y,psi_syz_y,psi_sxz_z,psi_syz_z,psi_szz_z,
+			psi_vxx,psi_vyx,psi_vzx,psi_vxy,psi_vyy,psi_vzy,psi_vxz,psi_vyz,psi_vzz);
 		MPI_Barrier(MPI_COMM_WORLD);
 		if (MYID==0){
 			time4=MPI_Wtime();
@@ -808,20 +810,20 @@ int main(int argc, char **argv){
 		}
         
         /* initialize wavefield with zero */
-        if (L==1 && ABS_TYPE==2) {
+        if ((L==1) && (ABS_TYPE==2) && (CHECKPTREAD==0)){
             zero(1-FDORDER/2,NX+FDORDER/2,1-FDORDER/2,NY+FDORDER/2,1-FDORDER/2,NZ+FDORDER/2,vx,vy,vz,sxx,syy,szz,sxy,syz,sxz,
                  vxyyx,vyzzy,vxzzx,vxxyyzz,vyyzz,vxxzz,vxxyy,vxyyx_2,vyzzy_2,vxzzx_2,vxxyyzz_2,vyyzz_2,vxxzz_2,vxxyy_2,vxyyx_3,
                  vyzzy_3,vxzzx_3,vxxyyzz_3,vyyzz_3,vxxzz_3,vxxyy_3,vxyyx_4,vyzzy_4,vxzzx_4,
                  vxxyyzz_4,vyyzz_4,vxxzz_4,vxxyy_4,svx,svy,svz,svx_2,svy_2,svz_2,svx_3,svy_3,svz_3,svx_4,svy_4,svz_4,rxx,ryy,rzz,rxy,ryz,rxz,rxx_2,ryy_2,rzz_2,rxy_2,ryz_2,rxz_2,rxx_3,ryy_3,rzz_3,rxy_3,ryz_3,rxz_3,rxx_4,ryy_4,rzz_4,rxy_4,ryz_4,rxz_4);
         }
         
-		if(L==0 && ABS_TYPE==2){
+		if((L==0) && (ABS_TYPE==2) && (CHECKPTREAD==0)){
             zero_elastic(1-FDORDER/2,NX+FDORDER/2,1-FDORDER/2,NY+FDORDER/2,1-FDORDER/2,NZ+FDORDER/2,vx,vy,vz,sxx,syy,szz,sxy,syz,sxz,
                          vxyyx,vyzzy,vxzzx,vxxyyzz,vyyzz,vxxzz,vxxyy,vxyyx_2,vyzzy_2,vxzzx_2,vxxyyzz_2,vyyzz_2,vxxzz_2,vxxyy_2,vxyyx_3,
                          vyzzy_3,vxzzx_3,vxxyyzz_3,vyyzz_3,vxxzz_3,vxxyy_3,vxyyx_4,vyzzy_4,vxzzx_4,
                          vxxyyzz_4,vyyzz_4,vxxzz_4,vxxyy_4,svx,svy,svz,svx_2,svy_2,svz_2,svx_3,svy_3,svz_3,svx_4,svy_4,svz_4);
 		}
-		if(ABS_TYPE==1){
+		if((ABS_TYPE==1) && (CHECKPTREAD==0)){
 			zero_elastic_CPML(NX,NY,NZ,vx,vy,vz,sxx,syy,szz,sxy,syz,sxz,rxx,ryy,rzz,rxy,ryz,rxz,psi_sxx_x,psi_sxy_x,psi_sxz_x,psi_sxy_y,psi_syy_y,psi_syz_y,psi_sxz_z,psi_syz_z,psi_szz_z,psi_vxx,psi_vyx,psi_vzx,psi_vxy,psi_vyy,psi_vzy,psi_vxz,psi_vyz,psi_vzz,rxx_2,ryy_2,rzz_2,rxy_2,ryz_2,rxz_2,rxx_3,ryy_3,rzz_3,rxy_3,ryz_3,rxz_3,rxx_4,ryy_4,rzz_4,rxy_4,ryz_4,rxz_4);
 
 		}
@@ -1049,7 +1051,9 @@ int main(int argc, char **argv){
 			fprintf(FP," Saving wavefield to check-point file %s \n",CHECKPTFILE);
 		}
 
-		save_checkpoint(-1, NX+2, -1, NY+2, -1, NZ+2, vx,vy,vz,sxx,syy,szz,sxy,syz,sxz);
+		save_checkpoint(-1, NX+2, -1, NY+2, -1, NZ+2, vx,vy,vz,sxx,syy,szz,sxy,syz,sxz,rxx,ryy,rzz,rxy,ryz,rxz,
+				psi_sxx_x,psi_sxy_x,psi_sxz_x,psi_sxy_y,psi_syy_y,psi_syz_y,psi_sxz_z,psi_syz_z,psi_szz_z,
+				psi_vxx,psi_vyx,psi_vzx,psi_vxy,psi_vyy,psi_vzy,psi_vxz,psi_vyz,psi_vzz);
 		MPI_Barrier(MPI_COMM_WORLD);
 		if (MYID==0){
 			time4=MPI_Wtime();