Merge branch 'develop' into pup

README.md

@@ -10,7 +10,7 @@ The [**manual**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/wikis/home) is in
 
 # Download and Newsletter
 
-You can download the [**latest Release**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tags/Release_2.0.1) or the current [**Beta-Version**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tree/develop).
+You can download the [**latest Release 2.0.2**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tags/Release_2.0.2) or the current [**Beta-Version**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tree/develop).
 
 To receive news and updates please [register](http://www.gpi.kit.edu/Software-FWI.php) on the email list [IFOS@lists.kit.edu](http://www.gpi.kit.edu/Software-FWI.php).  
 Please use this list also to ask questions or to report problems or bugs.
\ No newline at end of file

doc/latex/6_Parameter_Definition.tex

@@ -230,7 +230,7 @@ Default value is:
 
 With this option pure acoustic modelling and/or inversion can be performed (ACOUSTIC = 1). Only a P-wave and a density model need to be provided. Acoustic modelling and inversion can be a quick estimate, especially for marine environments.
 
-For acoustic modelling the option VELOCITY is not available and only PARAMETERIZATION = 1 is possible.
+For acoustic modelling only PARAMETERIZATION = 1 is possible.
 
 \section{PSV and SH modelling}
 {\color{blue}{\begin{verbatim}
@@ -322,7 +322,9 @@ The locations of the receivers may either be specified in a separate file REC\_F
 \end{verbatim}}}
 These receiver coordinates in REC\_FILE are shifted by REFREC[1], REFREC[2] into the  x- and y-direction, respectively. This allows for completely moving the receiver spread without modifying REC\_FILE. This may be useful for the simulation of moving profiles in reflection seismics.
 
-If READREC=0 the receiver locations must be specified in the parameter file. In this case, it is assumed that the receivers are located along a straight line. The first receiver position is defined by (XREC1, YREC1), and the last receiver position by (XREC1, YREC1). The spacing between receivers is NGEOPH grid points. 
+If READREC=0 the receiver locations must be specified in the parameter file. In this case, it is assumed that the receivers are located along a straight line. The first receiver position is defined by (XREC1, YREC1), and the last receiver position by (XREC1, YREC1). The spacing between receivers is NGEOPH grid points.
+
+READREC=2 is made for a moving streamer aquisition geometry. The receiver positions for shot 1 must be specified in the parameter file (as described for READREC=0). For all other shots the receiver positions are moved according to the movement of the following source. This implementation allows a memory saving usage of moving geometries. The receivers must be in the same depth and all receivers must be located in the model for all shots.
 
 Receivers are always located on full grid indices, i.e. a receiver that is located between two grid points will be shifted by the FD program to the closest next grid point. It is not yet possible
 to output seismograms for arbitrary receiver locations since this would require a certain wavefield interpolation.
@@ -504,12 +506,12 @@ Default values are:
 With the use of a workflow file, you can define different FWI stages. For instance, one FWI stage could refer to one corner frequency of a low-pass filter or to a specific time window. Every line in the workflow file corresponds to one FWI stage. To use a workflow file the switch USE\_WORKFLOW have to be set to 1.  The algorithm will automatically change to the next line of the workflow file if the abort criterium of the current line is reached or if no step length could be found which reduces the misfit. The structure of the variables inside the workflow file is as follow: \\
 
 \noindent
-\begin{tabular}{llllllllllll}
-\# & INV\_VS & INV\_VP & INV\_RHO & PRO & TIME\_FILT & FC & WAVETYPE & 0 & 0 & EPRECOND & EPSILON\_WE\\
+\begin{tabular}{lllllllllllll}
+\# & INV\_VS & INV\_VP & INV\_RHO & PRO & TIME\_FILT & F\_HIGH\_PASS & F\_LOW\_PASS & WAVETYPE & 0 & 0 & EPRECOND & EPSILON\_WE\\
 \end{tabular}
 
 \ \\
-The first column is the number of the line. With INV\_* etc. you can activate the inversion for VS, VP or RHO, respectively. The abort criterium in percent for this FWI stage will be the declared in the variable PRO. With TIME\_FILT you can activate the frequency filtering with the corner frequency FC. WAVETYPE can be used to switch between PSV and SH modeling, however it is recommended to use PSV modeling (WAVETYPE==1) due to the current development on SH modeling. The following to zeros are placeholders for an upcoming update. With EPRECOND and EPSILON\_WE you can control the approx. Hessian. Please note, that all features which are used eg. TIME\_FILT (see section \ref{sec:filtering}) within the workflow have to be activated in the .JSON file.
+The first column is the number of the line. With INV\_* etc. you can activate the inversion for VS, VP or RHO, respectively. The abort criterium in percent for this FWI stage will be the declared in the variable PRO. With TIME\_FILT you can activate the frequency filtering with the corner frequency F\_HIGH\_PASS of the high-pass and F\_LOW\_PASS of the low-pass filter. WAVETYPE can be used to switch between PSV and SH modeling, however it is recommended to use PSV modeling (WAVETYPE==1) due to the current development on SH modeling. The following to zeros are placeholders for an upcoming update. With EPRECOND and EPSILON\_WE you can control the approx. Hessian. Please note, that all features which are used eg. TIME\_FILT (see section \ref{sec:filtering}) within the workflow have to be activated in the .JSON file.
 
 For an example of a workflow file, have a look in the par/ folder.
 
@@ -857,7 +859,7 @@ For GRAD\_FILT\_WAVELENGTH = 1 (and TIME\_FILT=1) a new wavelength dependent fil
 \begin{equation}
  \mbox{FILT\_SIZE\_GRAD} = \frac{V_{s,\text{average}}\cdot \text{A}}{\text{FC}}
 \end{equation}
-where FC is the corner frequency of TIME\_FILT and A is a weighting factor.
+where F\_LOW\_PASS is the corner frequency of TIME\_FILT and A is a weighting factor.
 
 
 \section{Model manipulation}

genmod/toy_example_true_trench.c

+/*-----------------------------------------------------------------------------------------
+ * Copyright (C) 2005  <name of author>
+ *
+ * This file is part of DENISE.
+ * 
+ * DENISE 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.
+ * 
+ * DENISE 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 DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
+-----------------------------------------------------------------------------------------*/
+
+/* $Id: hh_el.c,v 2.3 2007/08/21 13:16:19 tbohlen Exp $*/
+/*
+ *   Model defined by flnode file. In this version a constant Q-model with Q_p=Q_s is used. It 
+ *   is defined via the relaxation frequencies and the TAU value given in the input file.  
+ */
+
+#include "fd.h"
+
+void model(float  **  rho, float **  pi, float **  u, float **  taus, float **  taup, float *  eta){
+
+	/*--------------------------------------------------------------------------*/
+	/* extern variables */
+
+	extern int NX, NY, NXG, NYG,  POS[3], MYID, L;
+	extern float DH, DT, *FL, TAU;
+	extern char  MFILE[STRING_SIZE];
+	extern char INV_MODELFILE[STRING_SIZE];	
+
+		/* local variables */
+	float muv, piv, vp, vs, rhov, ts, tp, *pts;
+	int i, j, ii, jj, l;
+	char modfile[STRING_SIZE];
+	
+	FILE *flfile;
+	int nodes;
+	char cline[256];
+	
+	float *fldepth, *flrho, *flvp, *flvs;
+	
+	/*---------------------------------------
+	 Trench "Ettlinger Linie" units in m  
+	----------------------------------------*/
+
+	/* trench=1 - elliptic shaped trench , trench=2 - V shaped trench */
+	int trench=2;
+
+	float width=10;
+	float depth=3;	
+	float pos=NXG/2*DH;
+
+	// Parameters of trench "Ettlinger Linie" velocities in km/s and density in kg/m^3
+	float evs=0.0917;
+	float evp=0.2569;
+	float erho=1.5;	
+	
+	/* fopen FL-node file*/
+	flfile=fopen("model_true/flnodes.toy_example","r");
+        if (flfile==NULL) declare_error(" FL-file could not be opened !");
+	/* nodes specified in FL-node File*/
+	nodes=7;
+
+	/*-----------No changes below that line needed ---------------------------*/
+
+	/* 2 linear lines as border for the v-shape */
+	float m1=2*depth/width;
+	float m2=-m1;
+	float c1=depth*(1-(2*pos/width));
+	float c2=depth*(1+(2*pos/width));
+
+	
+
+	
+	fldepth=vector(1,nodes);
+	flrho=vector(1,nodes);
+	flvp=vector(1,nodes);
+	flvs=vector(1,nodes);
+		
+	
+	/* vector for maxwellbodies */
+	pts=vector(1,L);
+	for (l=1;l<=L;l++) {
+		pts[l]=1.0/(2.0*PI*FL[l]);
+		eta[l]=DT/pts[l];
+	}
+	 
+	
+	
+	/*read FL nodes from File*/
+	
+	
+	for (l=1;l<=nodes;l++){
+		fgets(cline,255,flfile);
+		if (cline[0]!='#'){
+			sscanf(cline,"%f%f%f%f",&fldepth[l], &flrho[l], &flvp[l], &flvs[l]);
+		}
+		else l=l-1;
+	
+	}
+	
+	
+	if(MYID==0){
+	printf(" ------------------------------------------------------------------ \n\n");
+	printf(" Information of FL nodes: \n\n");
+	printf(" \t depth \t rho \t vp \t vs \n\n");
+	
+	for (l=1;l<=nodes;l++){
+	printf(" \t %f \t %f \t %f \t %f\n\n",fldepth[l],flrho[l],flvp[l],flvs[l]);
+	}
+	printf(" ------------------------------------------------------------------ \n\n");
+	}
+	/*-----------------------------------------------------------------------*/
+	
+	
+		
+	/*-----------------------------------------------------------------------*/
+
+	/* loop over global grid */
+		for (i=1;i<=NXG;i++){
+			for (l=1;l<nodes;l++){
+				if (fldepth[l]==fldepth[l+1]){
+					if ((i==1) && (MYID==0)){
+					printf("depth: %f m: double node\n",fldepth[l]);}}
+				else{
+					for (j=(int)(fldepth[l]/DH)+1;j<=(int)(fldepth[l+1]/DH);j++){
+			
+						vp=0.0; vs=0.0; rhov=0.0;
+				  
+						vp=(DH*(j-1)-fldepth[l])*(flvp[l+1]-flvp[l])/(fldepth[l+1]-fldepth[l])+flvp[l];
+						vp=vp*1000.0;
+						vs=(DH*(j-1)-fldepth[l])*(flvs[l+1]-flvs[l])/(fldepth[l+1]-fldepth[l])+flvs[l];
+						vs=vs*1000.0;
+						rhov=(DH*(j-1)-fldepth[l])*(flrho[l+1]-flrho[l])/(fldepth[l+1]-fldepth[l])+flrho[l];
+						rhov=rhov*1000.0;
+						
+						
+						muv=vs;
+						piv=vp;
+						ts=TAU;
+						tp=TAU;
+						
+						/* only the PE which belongs to the current global gridpoint 
+				  		is saving model parameters in his local arrays */
+						if ((POS[1]==((i-1)/NX)) && 
+				    		(POS[2]==((j-1)/NY))){
+						ii=i-POS[1]*NX;
+						jj=j-POS[2]*NY;
+
+						u[jj][ii]=muv;
+						rho[jj][ii]=rhov;
+						pi[jj][ii]=piv;
+						taus[jj][ii]=ts;
+						taup[jj][ii]=tp;
+						}
+					}
+				}
+			}
+			
+			for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
+			  
+				vp=0.0; vs=0.0; rhov=0.0;
+				vp=flvp[nodes]*1000.0; vs=flvs[nodes]*1000.0; rhov=flrho[nodes]*1000.0;
+				
+				muv=vs;
+				piv=vp;
+				ts=TAU;
+				tp=TAU;
+
+				/* only the PE which belongs to the current global gridpoint 
+				  		is saving model parameters in his local arrays */
+						if ((POS[1]==((i-1)/NX)) && 
+				    		(POS[2]==((j-1)/NY))){
+						ii=i-POS[1]*NX;
+						jj=j-POS[2]*NY;
+
+						u[jj][ii]=muv;
+						rho[jj][ii]=rhov;
+						pi[jj][ii]=piv;
+						taus[jj][ii]=ts;
+						taup[jj][ii]=tp;
+						}
+			}
+			/*Trench Ettlinger Linie*/
+			if (trench==1) {
+			 	for (j=1;j<=(int) depth/DH+1;j++){
+				
+					if ( ( ((pos-i*DH) * (pos-i*DH))/(width/2*width/2) + ((j*DH)*(j*DH))/(depth*depth) ) <=1){
+                                	vp=evp*1000.0; 
+					vs=evs*1000.0; 
+					rhov=erho*1000.0;
+				
+                                	muv=vs;
+                                	piv=vp;
+                                	ts=TAU;
+                                	tp=TAU;
+				
+
+                               	 	/* only the PE which belongs to the current global gridpoint 
+                                                is saving model parameters in his local arrays */
+                                                if ((POS[1]==((i-1)/NX)) &&
+                                                (POS[2]==((j-1)/NY))){
+                                                ii=i-POS[1]*NX;
+                                                jj=j-POS[2]*NY;
+
+                                                u[jj][ii]=muv;
+                                                rho[jj][ii]=rhov;
+                                                pi[jj][ii]=piv;
+                                                taus[jj][ii]=ts;
+                                                taup[jj][ii]=tp;
+                                                }
+                        		}
+				}
+
+			}
+			 /*Trench V-shape Ettlinger Linie*/
+                        if (trench==2) {
+                                for (j=1;j<=(int) depth/DH+1;j++){
+					
+                                        if ( ((i*DH)*m1+c1 >=0) && j*DH <= (i*DH)*m1+c1  && j*DH <= (i*DH)*m2+c2){
+                                        vp=evp*1000.0;
+                                        vs=evs*1000.0;
+                                        rhov=erho*1000.0;
+
+                                        muv=vs;
+                                        piv=vp;
+                                        ts=TAU;
+                                        tp=TAU;
+
+
+                                        /* only the PE which belongs to the current global gridpoint 
+                                                is saving model parameters in his local arrays */
+                                                if ((POS[1]==((i-1)/NX)) &&
+                                                (POS[2]==((j-1)/NY))){
+                                                ii=i-POS[1]*NX;
+                                                jj=j-POS[2]*NY;
+
+                                                u[jj][ii]=muv;
+                                                rho[jj][ii]=rhov;
+                                                pi[jj][ii]=piv;
+                                                taus[jj][ii]=ts;
+                                                taup[jj][ii]=tp;
+                                                }
+                                        }
+                                }
+
+                        }
+
+		}
+		
+		
+        sprintf(modfile,"%s_rho_it0.bin",INV_MODELFILE);
+	writemod(modfile,rho,3);
+	MPI_Barrier(MPI_COMM_WORLD);
+	if (MYID==0) mergemod(modfile,3);
+	MPI_Barrier(MPI_COMM_WORLD); 
+	sprintf(modfile,"%s_rho_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
+	remove(modfile);
+	
+	sprintf(modfile,"%s_vs_it0.bin",INV_MODELFILE);
+	writemod(modfile,u,3);
+	MPI_Barrier(MPI_COMM_WORLD);
+	if (MYID==0) mergemod(modfile,3);
+	MPI_Barrier(MPI_COMM_WORLD); 
+	sprintf(modfile,"%s_vs_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
+	remove(modfile);
+	
+	sprintf(modfile,"%s_vp_it0.bin",INV_MODELFILE);
+	writemod(modfile,pi,3);
+	MPI_Barrier(MPI_COMM_WORLD);
+	if (MYID==0) mergemod(modfile,3);
+	MPI_Barrier(MPI_COMM_WORLD); 
+	sprintf(modfile,"%s_vp_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
+	remove(modfile);
+	
+	sprintf(modfile,"%s_taus_it0.bin",INV_MODELFILE);
+	writemod(modfile,taus,3);
+	MPI_Barrier(MPI_COMM_WORLD);
+	if (MYID==0) mergemod(modfile,3);
+	MPI_Barrier(MPI_COMM_WORLD); 
+	sprintf(modfile,"%s_taus_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
+	remove(modfile);
+	
+	sprintf(modfile,"%s_taup_it0.bin",INV_MODELFILE);
+	writemod(modfile,taup,3);
+	MPI_Barrier(MPI_COMM_WORLD);
+	if (MYID==0) mergemod(modfile,3);
+	MPI_Barrier(MPI_COMM_WORLD); 
+	sprintf(modfile,"%s_taup_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
+	remove(modfile);
+	
+	
+	free_vector(fldepth,1,nodes);
+	free_vector(flrho,1,nodes);
+	free_vector(flvp,1,nodes);
+	free_vector(flvs,1,nodes);
+	
+	free_vector(pts,1,L);
+	 
+}
+
+
+

par/in_and_out/IFOS2D_INV_all_parameters.json

@@ -177,7 +177,6 @@
 			"F_LOW_PASS_END" : "75.0",
 			"F_LOW_PASS_INCR" : "10.0",
 			"ORDER" : "2",
-			"ZERO_PHASE" : "0",
 			"FREQ_FILE" : "frequencies.dat",  
 			"WRITE_FILTERED_DATA" : "0",                     
 

par/run_toy_example.sh

@@ -7,6 +7,14 @@
 # compiling all libraries and  IFOS
 make install MODEL=../genmod/toy_example_true.c
 
+#-------------------------------------------------------------------------------------
+# Insert the trench "Ettlinger Linie" into the toy_example model 
+# in the file ../genmod/toy_example_true_trench.c the size and position can be adjusted
+# Either an elliptic trench or a v-shaped trench can be choosen
+
+#make install MODEL=../genmod/toy_example_true_trench.c
+#---------------------------------------------------------------------------------------
+
 # starting IFOS for forward modeling
 # (depending on the MPI package you maybe have to adjust the following programme call, 
 #  e.g. when using openMPI you do not have to use the 'lamboot' command)

par/workflow.txt

-#	VS	VP	RHO	PRO	FIL	FC	WT	J	J	EPRE	EPSI
-1	1	0	0	0.01	0	0	2	0	0	0	0.005
-2	1	0	0	0.01	0	0	2	0	0	0	0.005	
-3	1	0	0	0.01	0	0	2	0	0	0	0.005
-4	1	0	0	0.01	0	0	2	0	0	0	0.005
+#	VS	VP	RHO	PRO	F_FILT	F_HIGH	F_LOW	WT	J	J	EPRE	EPSI
+1	1	0	0	0.01	0	0	0	2	0	0	0	0.005
+2	1	0	0	0.01	0	0	0	2	0	0	0	0.005	
+3	1	0	0	0.01	0	0	0	2	0	0	0	0.005
+4	1	0	0	0.01	0	0	0	2	0	0	0	0.005

src/Makefile

@@ -70,6 +70,7 @@ IFOS2D= \
 	IFOS2D.c \
 	stf.c \
 	window_cos.c \
+	alloc_sections.c \
 	calc_mat_change_test.c \
     	calc_res.c \
 	calc_misfit.c \
@@ -208,4 +209,4 @@ clean:
 
 install: clean all
 
--include $(IFOS_OBJ:.o=.d)
\ No newline at end of file
+-include $(IFOS_OBJ:.o=.d)

src/alloc_sections.c

+/*-----------------------------------------------------------------------------------------
+ * Copyright (C) 2016  For the list of authors, see file AUTHORS.
+ *
+ * This file is part of IFOS.
+ *
+ * IFOS 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.
+ *
+ * IFOS 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 IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
+ -----------------------------------------------------------------------------------------*/
+
+/*  Computation of local source coordinates
+ *
+ *
+ */
+#include "fd.h"
+ 
+void alloc_sections(int ntr,int ns,float ***sectionvx,float ***sectionvy,float ***sectionvz,float ***sectionp,float ***sectionpnp1,float ***sectionpn,float ***sectioncurl,float ***sectiondiv,
+                    float ***sectionpdata,float ***sectionpdiff,float ***sectionpdiffold,float ***sectionvxdata,float ***sectionvxdiff,float ***sectionvxdiffold,float ***sectionvydata,
+                    float ***sectionvydiff,float ***sectionvydiffold,float ***sectionvzdata,float ***sectionvzdiff,float ***sectionvzdiffold)
+{
+	extern int SEISMO, WAVETYPE;
+	extern FILE *FP;
+	
+	switch (SEISMO){
+            case 1 : /* particle velocities only */
+                switch (WAVETYPE) {
+                    case 1:
+                        *sectionvx=matrix(1,ntr,1,ns);
+                        *sectionvy=matrix(1,ntr,1,ns);
+                        break;
+                    case 2:
+                        *sectionvz=matrix(1,ntr,1,ns);
+                        break;
+                    case 3:
+                        *sectionvx=matrix(1,ntr,1,ns);
+                        *sectionvy=matrix(1,ntr,1,ns);
+                        *sectionvz=matrix(1,ntr,1,ns);
+                        break;
+                }
+                break;
+            case 2 : /* pressure only */
+                *sectionp=matrix(1,ntr,1,ns);
+                *sectionpnp1=matrix(1,ntr,1,ns);
+                *sectionpn=matrix(1,ntr,1,ns);
+                break;
+            case 3 : /* curl and div only */
+                *sectioncurl=matrix(1,ntr,1,ns);
+                *sectiondiv=matrix(1,ntr,1,ns);
+                break;
+            case 4 : /* everything */
+                switch (WAVETYPE) {
+                    case 1:
+                        *sectionvx=matrix(1,ntr,1,ns);
+                        *sectionvy=matrix(1,ntr,1,ns);
+                        break;
+                    case 2:
+                        *sectionvz=matrix(1,ntr,1,ns);
+                        break;
+                    case 3:
+                        *sectionvx=matrix(1,ntr,1,ns);
+                        *sectionvy=matrix(1,ntr,1,ns);
+                        *sectionvz=matrix(1,ntr,1,ns);
+                        break;
+                }
+                *sectioncurl=matrix(1,ntr,1,ns);
+                *sectiondiv=matrix(1,ntr,1,ns);
+                *sectionp=matrix(1,ntr,1,ns);
+                break;
+            case 5 : /* everything except curl and div*/
+                switch (WAVETYPE) {
+                    case 1:
+                        *sectionvx=matrix(1,ntr,1,ns);
+                        *sectionvy=matrix(1,ntr,1,ns);
+                        break;
+                    case 2:
+                        *sectionvz=matrix(1,ntr,1,ns);
+                        break;
+                    case 3:
+                        *sectionvx=matrix(1,ntr,1,ns);
+                        *sectionvy=matrix(1,ntr,1,ns);
+                        *sectionvz=matrix(1,ntr,1,ns);
+                        break;
+                }
+                *sectionp=matrix(1,ntr,1,ns);
+                break;
+        }
+        
+        *sectionpdata=matrix(1,ntr,1,ns);
+    *sectionpdiff=matrix(1,ntr,1,ns);
+    *sectionpdiffold=matrix(1,ntr,1,ns);
+    switch (WAVETYPE) {
+        case 1:
+            *sectionvxdata=matrix(1,ntr,1,ns);
+            *sectionvxdiff=matrix(1,ntr,1,ns);
+            *sectionvxdiffold=matrix(1,ntr,1,ns);
+            *sectionvydata=matrix(1,ntr,1,ns);
+            *sectionvydiff=matrix(1,ntr,1,ns);
+            *sectionvydiffold=matrix(1,ntr,1,ns);
+            break;
+            
+        case 2:
+            *sectionvzdata=matrix(1,ntr,1,ns);
+            *sectionvzdiff=matrix(1,ntr,1,ns);
+            *sectionvzdiffold=matrix(1,ntr,1,ns);
+            break;
+            
+        case 3:
+            *sectionvxdata=matrix(1,ntr,1,ns);
+            *sectionvxdiff=matrix(1,ntr,1,ns);
+            *sectionvxdiffold=matrix(1,ntr,1,ns);
+            *sectionvydata=matrix(1,ntr,1,ns);
+            *sectionvydiff=matrix(1,ntr,1,ns);
+            *sectionvydiffold=matrix(1,ntr,1,ns);
+            *sectionvzdata=matrix(1,ntr,1,ns);
+            *sectionvzdiff=matrix(1,ntr,1,ns);
+            *sectionvzdiffold=matrix(1,ntr,1,ns);
+            break;
+    }
+ }
+ 
+ 
+ 
+ void dealloc_sections(int ntr,int ns,int **recpos_loc,float **sectionvx,float **sectionvy,float **sectionvz,float **sectionp,float **sectionpnp1,float **sectionpn,float **sectioncurl,float **sectiondiv,
+                    float **sectionpdata,float **sectionpdiff,float **sectionpdiffold,float **sectionvxdata,float **sectionvxdiff,float **sectionvxdiffold,float **sectionvydata,
+                    float **sectionvydiff,float **sectionvydiffold,float **sectionvzdata,float **sectionvzdiff,float **sectionvzdiffold)
+{
+	extern int SEISMO, WAVETYPE;
+	extern FILE *FP;
+	
+	free_imatrix(recpos_loc,1,3,1,ntr);
+        switch (SEISMO){
+            case 1 : /* particle velocities only */
+                if (WAVETYPE==1 || WAVETYPE==3) {
+                    free_matrix(sectionvx,1,ntr,1,ns);
+                    free_matrix(sectionvy,1,ntr,1,ns);
+                }
+                if (WAVETYPE==2 || WAVETYPE==3) {
+                    free_matrix(sectionvz,1,ntr,1,ns);
+                }
+                break;
+            case 2 : /* pressure only */
+                if (WAVETYPE==1 || WAVETYPE==3) {
+                    free_matrix(sectionp,1,ntr,1,ns);
+                    free_matrix(sectionpn,1,ntr,1,ns);
+                    free_matrix(sectionpnp1,1,ntr,1,ns);
+                }
+                break;
+            case 3 : /* curl and div only */
+                if (WAVETYPE==1 || WAVETYPE==3) {
+                    free_matrix(sectioncurl,1,ntr,1,ns);
+                    free_matrix(sectiondiv,1,ntr,1,ns);
+                }
+                break;
+            case 4 : /* everything */
+                if (WAVETYPE==1 || WAVETYPE==3) {
+                    free_matrix(sectionvx,1,ntr,1,ns);
+                    free_matrix(sectionvy,1,ntr,1,ns);
+                    free_matrix(sectionp,1,ntr,1,ns);
+                    free_matrix(sectioncurl,1,ntr,1,ns);
+                    free_matrix(sectiondiv,1,ntr,1,ns);
+                }
+                if (WAVETYPE==2 || WAVETYPE==3) {
+                    free_matrix(sectionvz,1,ntr,1,ns);
+                }
+                break;
+            case 5 : /* everything except curl and div */
+                if (WAVETYPE==1 || WAVETYPE==3) {
+                    free_matrix(sectionvx,1,ntr,1,ns);
+                    free_matrix(sectionvy,1,ntr,1,ns);
+                    free_matrix(sectionp,1,ntr,1,ns);
+                }
+                if (WAVETYPE==2 || WAVETYPE==3) {
+                    free_matrix(sectionvz,1,ntr,1,ns);
+                }
+                break;
+        }
+    
+    
+    if (WAVETYPE==1 || WAVETYPE==3) {
+        free_matrix(sectionvxdata,1,ntr,1,ns);
+        free_matrix(sectionvxdiff,1,ntr,1,ns);
+        free_matrix(sectionvydata,1,ntr,1,ns);
+        free_matrix(sectionvydiff,1,ntr,1,ns);
+        free_matrix(sectionvydiffold,1,ntr,1,ns);
+        free_matrix(sectionvxdiffold,1,ntr,1,ns);
+        free_matrix(sectionpdata,1,ntr,1,ns);
+        free_matrix(sectionpdiff,1,ntr,1,ns);
+        free_matrix(sectionpdiffold,1,ntr,1,ns);
+    }
+    if (WAVETYPE==2 || WAVETYPE==3) {
+        free_matrix(sectionvzdata,1,ntr,1,ns);
+        free_matrix(sectionvzdiff,1,ntr,1,ns);