taper_grad.c 20 KB
Newer Older
Tilman Steinweg's avatar
Tilman Steinweg committed
1
/*-----------------------------------------------------------------------------------------
Florian Wittkamp's avatar
Florian Wittkamp committed
2
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
Tilman Steinweg's avatar
Tilman Steinweg committed
3
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
4 5 6
 * This file is part of IFOS.
 * 
 * IFOS is free software: you can redistribute it and/or modify
Tilman Steinweg's avatar
Tilman Steinweg committed
7 8
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
Florian Wittkamp's avatar
Florian Wittkamp committed
9 10
 * 
 * IFOS is distributed in the hope that it will be useful,
Tilman Steinweg's avatar
Tilman Steinweg committed
11 12 13
 * 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.
Florian Wittkamp's avatar
Florian Wittkamp committed
14
 * 
Tilman Steinweg's avatar
Tilman Steinweg committed
15
 * You should have received a copy of the GNU General Public License
Florian Wittkamp's avatar
Florian Wittkamp committed
16 17
 * along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
Tilman Steinweg's avatar
Tilman Steinweg committed
18 19 20 21 22 23 24 25 26 27 28 29 30

/*
 * taper gradient with a gaussian frame to damp inversion artefacts near the sources and receivers
 sws == 1 vertical taper (for tomography geometry)
 sws == 2 horizontal taper (for reflection geometry)
 sws == 3 local radial taper at the source and receiver positions
 sws == 4 custom taper
 */

#include "fd.h"

void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nshots, int **recpos, int ntr, int sws)
{
Florian Wittkamp's avatar
Florian Wittkamp committed
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

	/* extern variables */

        extern float DH;
	extern int FREE_SURF, NX, NY, NXG, NYG;
	extern int NPROCX, NPROCY, MYID, POS[3];
	extern FILE *FP;
	extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_U[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
	
	/* local variables */
	int i, j, h, ifw, ii, jj, n, xb, yb, xe, ye, taperlength,taperlength2, VTON, SRTON;
	int ijc, iy, ix, iii, jjj, xx, yy, srctaper_gridpt, i1, j1;

	extern int GRADT1, GRADT2, GRADT3, GRADT4;
	float amp, a, *window, grad_tap, **waveconvtmp;
	char modfile[STRING_SIZE];
	
	extern float SRTRADIUS;
	extern int SRTSHAPE, FILTSIZE;
        float **m, **edgemat, **mm, **msum, minm, maxm, x, y, rad, **taper_coeff_glob;
        float maxrad;
	char taper_file[STRING_SIZE];
53
    extern int VERBOSE;
Florian Wittkamp's avatar
Florian Wittkamp committed
54 55 56 57 58 59
	FILE *fp_taper;

	/*SRTSHAPE=2;
	SRTRADIUS=25.0;
	filtsize=2;*/

Tilman Steinweg's avatar
Tilman Steinweg committed
60 61 62
    /* =============== */
    /* Vertical taper  */
    /* =============== */
Florian Wittkamp's avatar
Florian Wittkamp committed
63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
	
	if(sws==1){
	
	/* define taper geometry */
	taperlength=GRADT2-GRADT1;
	taperlength2=GRADT4-GRADT3;
	ifw = GRADT4-GRADT1;
	
	/*printf("%d \t %d \t %d \t %d \n",GRADT1, GRADT2, GRADT3, GRADT4);
	printf("%d \t %d \t %d \n",taperlength, taperlength2,ifw);*/
	
	if (MYID==0&&VERBOSE)
	{
		fprintf(FP,"\n **Message from taper_grid (printed by PE %d):\n",MYID);
		fprintf(FP," Coefficients for gradient taper are now calculated.\n");
	}
	
	waveconvtmp = matrix(0,NY+1,0,NX+1);
	window=vector(1,ifw);
	 
	 /* Gaussian window */  
Tilman Steinweg's avatar
Tilman Steinweg committed
84 85
        a=3;    /* damping coefficient */
        for (i=1;i<=ifw;i++){
Florian Wittkamp's avatar
Florian Wittkamp committed
86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
	window[i] = 1.0;
	
	    if(i<=taperlength){
	       window[i] = exp(-(1.0/2.0)*(a*(i-taperlength)/(taperlength/2.0))*(a*(i-taperlength)/(taperlength/2.0)));
	    }
	    
	    if(i>=(ifw-taperlength2)){
	       window[i] = exp(-(1.0/2.0)*(a*(i-(ifw-taperlength2))/(taperlength2/2.0))*(a*(i-(ifw-taperlength2))/(taperlength2/2.0)));
	    }
	    
	}
	
	/* loop over global grid */
	for (j=1;j<=NYG;j++){
	
	   h=1;
	   for (i=1;i<=NXG;i++){

                        grad_tap=0.0;
			
			if((i>GRADT1)&&(i<GRADT4)){
			   grad_tap=window[h];
			   h++;
			}
			   			
			if ((POS[1]==((i-1)/NX)) && 
		   	 (POS[2]==((j-1)/NY))){
				ii=i-POS[1]*NX;
				jj=j-POS[2]*NY;

				taper_coeff[jj][ii]=grad_tap;

			}
		}
	}	
	
	/* apply taper on local gradient */
	for (j=1;j<=NY;j++){
	   for (i=1;i<=NX;i++){
           waveconv[j][i]*=taper_coeff[j][i];
	   waveconvtmp[j][i] = waveconv[j][i];
           }
	}
	
		
Tilman Steinweg's avatar
Tilman Steinweg committed
131
        /* apply filter at shot and receiver points */
Florian Wittkamp's avatar
Florian Wittkamp committed
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189
	for (n=1;n<=nshots;n++)
	{
		i = iround(srcpos[1][n]/DH);
		j = iround(srcpos[2][n]/DH);
		if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
			ii = i-POS[1]*NX;
			jj = j-POS[2]*NY;
			/*waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
						+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
						+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
			waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
			
			waveconvtmp[jj][ii] = 0.0;
		}
	}

	for (j=1;j<=NY;j++){
		for (i=1;i<=NX;i++){
        		waveconv[j][i] = waveconvtmp[j][i];
		}
	}	
	
	for (n=1;n<=ntr;n++)
	{
		i = recpos[1][n];
		j = recpos[2][n];
		if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
			ii = i-POS[1]*NX;
			jj = j-POS[2]*NY;
			/*waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
						+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
						+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
			waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
			
			waveconvtmp[jj][ii] = 0.0;
			
		}
	}

	for (j=1;j<=NY;j++){
		for (i=1;i<=NX;i++){
        		waveconv[j][i] = waveconvtmp[j][i];
		}
	}	

	sprintf(modfile,"taper_coeff_vert.bin");

	writemod(modfile,taper_coeff,3); 

	MPI_Barrier(MPI_COMM_WORLD);

	if (MYID==0) mergemod(modfile,3); 


	free_vector(window,1,ifw);
	free_matrix(waveconvtmp,0,NX+1,0,NY+1);
	}

190
        
Florian Wittkamp's avatar
Florian Wittkamp committed
191 192 193 194 195
	/* ======================== */
	/* Horizontal Taper         */
	/* ======================== */
	
	if(sws==2){
196 197 198 199
        /* define taper geometry */
        taperlength=GRADT2-GRADT1;
        taperlength2=GRADT4-GRADT3;
        ifw = GRADT2-GRADT1+1;
Florian Wittkamp's avatar
Florian Wittkamp committed
200

201 202 203 204 205 206
        if (MYID==0&&VERBOSE){
            printf("%d \t %d \t %d \t %d \n",GRADT1, GRADT2, GRADT3, GRADT4);
            printf("%d \t %d \t %d \n",taperlength, taperlength2,ifw);
            fprintf(FP,"\n **Message from taper_grid (printed by PE %d):\n",MYID);
            fprintf(FP," Coefficients for gradient taper are now calculated.\n");
        }
Florian Wittkamp's avatar
Florian Wittkamp committed
207

208 209
        waveconvtmp = matrix(0,NY+1,0,NX+1);
        window=vector(1,ifw);
Florian Wittkamp's avatar
Florian Wittkamp committed
210 211
         
        /* Gaussian window */  
212 213
        a=3;    /* damping coefficient */
        for (i=1;i<=ifw;i++){
Florian Wittkamp's avatar
Florian Wittkamp committed
214 215
        window[i] = 1.0;

216
            if(i<=taperlength){
Florian Wittkamp's avatar
Florian Wittkamp committed
217
               window[i] = exp(-(1.0/2.0)*(a*(i-taperlength)/(taperlength/2.0))*(a*(i-taperlength)/(taperlength/2.0)));
218
            }
Florian Wittkamp's avatar
Florian Wittkamp committed
219
          
220
        }
Florian Wittkamp's avatar
Florian Wittkamp committed
221

222
        /* loop over global grid */
Florian Wittkamp's avatar
Florian Wittkamp committed
223
	h=1;
224
        for (j=1;j<=NYG;j++){
Florian Wittkamp's avatar
Florian Wittkamp committed
225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
           for (i=1;i<=NXG;i++){

                        grad_tap=0.0;

                        if((j>=GRADT1)&&(j<=GRADT2)){
                           grad_tap=window[h];
                        }

                        if(j>GRADT2){
                          
                          /*grad_tap=((float)(j*DH))*((float)(j*DH))*((float)(j*DH));*/
			  /* grad_tap=((float)(j*DH))*((float)(j*DH));*/
			  
			  /*grad_tap=((float)(j*DH));*/    /*Version of Daniel */
			  grad_tap=1.0;

                        }
			
			  
                        if ((POS[1]==((i-1)/NX)) && 
                         (POS[2]==((j-1)/NY))){
                                ii=i-POS[1]*NX;
                                jj=j-POS[2]*NY;

                                taper_coeff[jj][ii]=grad_tap;

                        }

253
                }
Florian Wittkamp's avatar
Florian Wittkamp committed
254 255
		
	   if(j>=GRADT1){h++;}
256
        }
Florian Wittkamp's avatar
Florian Wittkamp committed
257

258 259
        /* apply taper on local gradient */
        for (j=1;j<=NY;j++){
Florian Wittkamp's avatar
Florian Wittkamp committed
260 261 262 263
           for (i=1;i<=NX;i++){
           waveconv[j][i]*=taper_coeff[j][i];
           waveconvtmp[j][i] = waveconv[j][i];
           }
264
        }
Florian Wittkamp's avatar
Florian Wittkamp committed
265 266 267



268 269
        /* apply filter at shot and receiver points */
        /*for (n=1;n<=nshots;n++)
Florian Wittkamp's avatar
Florian Wittkamp committed
270 271 272 273 274 275 276 277 278 279
        {
                i = iround(srcpos[1][n]/DH);
                j = iround(srcpos[2][n]/DH);
                if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
                        ii = i-POS[1]*NX;
                        jj = j-POS[2]*NY;
                        waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
                                                + 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
                                                + 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
                        waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;
Tilman Steinweg's avatar
Tilman Steinweg committed
280
                }
281
        }
Florian Wittkamp's avatar
Florian Wittkamp committed
282

Tilman Steinweg's avatar
Tilman Steinweg committed
283
        for (j=1;j<=NY;j++){
Florian Wittkamp's avatar
Florian Wittkamp committed
284 285 286
                for (i=1;i<=NX;i++){
                        waveconv[j][i] = waveconvtmp[j][i];
                }
287
        }
Florian Wittkamp's avatar
Florian Wittkamp committed
288 289

        for (n=1;n<=ntr;n++)
290
        {
Florian Wittkamp's avatar
Florian Wittkamp committed
291 292 293
                i = recpos[1][n];
                j = recpos[2][n];
                if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
294 295
                        ii = i-POS[1]*NX;
                        jj = j-POS[2]*NY;
Florian Wittkamp's avatar
Florian Wittkamp committed
296 297 298 299
                        waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
                                                + 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
                                                + 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
                        waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;
300 301
                }
        }
Florian Wittkamp's avatar
Florian Wittkamp committed
302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472

        for (j=1;j<=NY;j++){
                for (i=1;i<=NX;i++){
                        waveconv[j][i] = waveconvtmp[j][i];
                }
        }*/

        sprintf(modfile,"taper_coeff_hor.bin");
        writemod(modfile,taper_coeff,3); 
        MPI_Barrier(MPI_COMM_WORLD);
        if (MYID==0) mergemod(modfile,3); 


        free_vector(window,1,ifw);
        free_matrix(waveconvtmp,0,NX+1,0,NY+1);
        } /* end of sws==2 */

        /* =================================== */
        /* taper source and receiver positions */
	/* =================================== */
	
	if(sws==3) {
                /* Convert from meters to gridpoints -> minimum 5x5 gridpoints */
                srctaper_gridpt = (int)(ceil(2.0*SRTRADIUS/DH));
                if (srctaper_gridpt<5)  srctaper_gridpt = 5;

                m               = matrix(1,srctaper_gridpt,1,srctaper_gridpt);
                edgemat         = matrix(1,4,1,1);
                mm              = matrix(1,NYG,1,NXG);
                msum            = matrix(1,NYG,1,NXG);
                taper_coeff_glob= matrix(1,NYG,1,NXG);
		waveconvtmp     = matrix(0,NY+1,0,NX+1);

                for (iy=1;iy<=NYG;iy++)
                        for (ix=1;ix<=NXG;ix++)  msum[iy][ix] = 1.0;

                MPI_Barrier(MPI_COMM_WORLD);

                /*****************************/
                /* Taper at source positions */
                /*****************************/

                a = 1.0;
                maxrad = sqrt(2.0*SRTRADIUS*SRTRADIUS);
                for (j=1;j<=srctaper_gridpt;j++) {
                        for (i=1;i<=srctaper_gridpt;i++) {
                                x = ((float)i-((float)srctaper_gridpt)/2.0-0.5)*DH;
                                y = ((float)j-((float)srctaper_gridpt)/2.0-0.5)*DH;
                                rad = sqrt(x*x+y*y);

                                switch (SRTSHAPE) {
                                case 1:
                                        m[j][i] = erf(a*rad/maxrad);
                                        break;
                                case 2:
                                        if (rad>0)      m[j][i] = log(rad);
                                        else            m[j][i] = 0.0;
                                        break;
                                }
                        }
                }

                /* generate local taper matrix */
                minm = minimum_m(m,srctaper_gridpt,srctaper_gridpt);
                for (j=1;j<=srctaper_gridpt;j++)
                        for (i=1;i<=srctaper_gridpt;i++)  m[j][i] -= minm;

                /* normalize taper matrix to max of values at the centre of all 4 taper area edges,     */
                /* not the global maximum, which is located at the corners                              */
                edgemat[1][1] = m[1][srctaper_gridpt/2];
                edgemat[2][1] = m[srctaper_gridpt/2][1];
                edgemat[3][1] = m[srctaper_gridpt/2][srctaper_gridpt];
                edgemat[4][1] = m[srctaper_gridpt][srctaper_gridpt/2];
                maxm = maximum_m(edgemat,1,4);
                for (j=1;j<=srctaper_gridpt;j++)
                        for (i=1;i<=srctaper_gridpt;i++) {
                                m[j][i] /= maxm;
                                if (m[j][i]>1.0)  m[j][i] = 1.0;
                        }

                /* get central position within the taper */
                ijc = (int)(ceil((float)srctaper_gridpt/2));

                /*********************/
                /* loop over sources */
                for (n=1;n<=nshots;n++) {
                        for (iy=1;iy<=NYG;iy++)
                                for (ix=1;ix<=NXG;ix++)  mm[iy][ix] = 1.0;

                        i = iround(srcpos[1][n]/DH);
                        j = iround(srcpos[2][n]/DH);
                        for (iy=1;iy<=srctaper_gridpt;iy++) {
                                for (ix=1;ix<=srctaper_gridpt;ix++) {
                                        xx = i + ix - ijc;
                                        yy = j + iy - ijc;
                                        if ((xx<1) || (xx>NXG) || (yy<1) || (yy>NYG))  continue;
                                        mm[yy][xx] = m[iy][ix];
                                }
                        }

/*                      for (iy=1;iy<=NYG;iy++)
                                for (ix=1;ix<=NXG;ix++)  msum[iy][ix] += mm[iy][ix];
*/
                        for (iy=1;iy<=NYG;iy++)
                                for (ix=1;ix<=NXG;ix++)
                                        if (msum[iy][ix] > mm[iy][ix])
                                                msum[iy][ix] = mm[iy][ix];

                }

                /***********************/
                /* loop over receivers */
                /*for (n=1;n<=ntr;n++) {
                        for (iy=1;iy<=NYG;iy++)
                                for (ix=1;ix<=NXG;ix++)  mm[iy][ix] = 1.0;

                        i = recpos[1][n];
                        j = recpos[2][n];
                        for (iy=1;iy<=srctaper_gridpt;iy++) {
                                for (ix=1;ix<=srctaper_gridpt;ix++) {
                                        xx = i + ix - ijc;
                                        yy = j + iy - ijc;
                                        if ((xx<1) || (xx>NXG) || (yy<1) || (yy>NYG))  continue;
                                        mm[yy][xx] = m[iy][ix];
                                }
                        }*/

/*                      for (iy=1;iy<=NYG;iy++)    Die kommenden zwei Zeilen wurden von Daniel auskommentiert.
                                for (ix=1;ix<=NXG;ix++)  msum[iy][ix] += mm[iy][ix];
*/
                       /* for (iy=1;iy<=NYG;iy++)
                                for (ix=1;ix<=NXG;ix++)
                                        if (msum[iy][ix] > mm[iy][ix])
                                                msum[iy][ix] = mm[iy][ix];

                }*/

                minm = minimum_m(msum,NXG,NYG);
                for (iy=1;iy<=NYG;iy++)
                        for (ix=1;ix<=NXG;ix++)  msum[iy][ix] -= minm;

                maxm = maximum_m(msum,NXG,NYG);
                for (iy=1;iy<=NYG;iy++)
                        for (ix=1;ix<=NXG;ix++) {
                                taper_coeff_glob[iy][ix] = msum[iy][ix]/maxm;

                                if ((POS[1]==((ix-1)/NX)) && (POS[2]==((iy-1)/NY))){
                                        ii = ix-POS[1]*NX;
                                        jj = iy-POS[2]*NY;
                                        /*Diese Zeile wurde von Daniel auskommentiert: taper_coeff[jj][ii] = taper_coeff_glob[iy][ix] * ((float)(iy*DH));*/
					/*taper_coeff[jj][ii] = ((float)(iy*DH)) * ((float)(iy*DH)) * ((float)(iy*DH));*/
					taper_coeff[jj][ii]=msum[iy][ix]/maxm;
                                }
                        }
			
	/* apply taper on local gradient */
        for (j=1;j<=NY;j++){
           for (i=1;i<=NX;i++){
           waveconv[j][i]*=taper_coeff[j][i];
           waveconvtmp[j][i] = waveconv[j][i];
           }
        }	
		
	 /* apply filter at shot and receiver points */
	for (n=1;n<=nshots;n++)
	{
		i1 = iround(srcpos[1][n]/DH);
		j1 = iround(srcpos[2][n]/DH);
		
		for (i=i1-FILTSIZE;i<=i1+FILTSIZE;i++){
		   for (j=j1-FILTSIZE;j<=j1+FILTSIZE;j++){
Tilman Steinweg's avatar
Tilman Steinweg committed
473
		       if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
Florian Wittkamp's avatar
Florian Wittkamp committed
474 475 476 477 478
			     ii = i-POS[1]*NX;
			     jj = j-POS[2]*NY;
			     /*waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
				       		+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
						+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
Tilman Steinweg's avatar
Tilman Steinweg committed
479
			      waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
Florian Wittkamp's avatar
Florian Wittkamp committed
480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509
			      if (jj>0){
			         waveconvtmp[jj][ii] = 0.0;
			         taper_coeff[jj][ii] = 0.0;
			      }
		        }
		    }
		}
		
		
	}

	
	/*for (n=1;n<=ntr;n++)
	{
		i1 = recpos[1][n];
		j1 = recpos[2][n];
		
            for (i=i1-FILTSIZE;i<=i1+FILTSIZE;i++){
		   for (j=j1-FILTSIZE;j<=j1+FILTSIZE;j++){
		
		       if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
			    ii = i-POS[1]*NX;
			    jj = j-POS[2]*NY;   */
			    /* Die kommenden 4 Zeilen wurden von Daniel auskommentiert. waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
						+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
						+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
			      waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
			
		/*	waveconvtmp[jj][ii] = 0.0;
			
Tilman Steinweg's avatar
Tilman Steinweg committed
510
		       }
Florian Wittkamp's avatar
Florian Wittkamp committed
511 512 513 514
		   }
	     }	       
	}*/		
		
Tilman Steinweg's avatar
Tilman Steinweg committed
515 516
        /* apply taper on local gradient */
        for (j=1;j<=NY;j++){
Florian Wittkamp's avatar
Florian Wittkamp committed
517 518 519
           for (i=1;i<=NX;i++){
           waveconv[j][i] = waveconvtmp[j][i];
           }
Tilman Steinweg's avatar
Tilman Steinweg committed
520
        }
Florian Wittkamp's avatar
Florian Wittkamp committed
521 522 523 524 525 526 527 528 529 530 531 532 533 534 535
	
	
                free_matrix(m,1,srctaper_gridpt,1,srctaper_gridpt);
                free_matrix(edgemat,1,4,1,1);
                free_matrix(mm,1,NYG,1,NXG);
                free_matrix(msum,1,NYG,1,NXG);
                free_matrix(taper_coeff_glob,1,NYG,1,NXG);
		free_matrix(waveconvtmp,0,NX+1,0,NY+1);
        
	
	

 
		
	MPI_Barrier(MPI_COMM_WORLD);
Tilman Steinweg's avatar
Tilman Steinweg committed
536
        sprintf(modfile,"taper_coeff_sources.bin");
Florian Wittkamp's avatar
Florian Wittkamp committed
537
        writemod(modfile,taper_coeff,3); 
Tilman Steinweg's avatar
Tilman Steinweg committed
538
        MPI_Barrier(MPI_COMM_WORLD);
Florian Wittkamp's avatar
Florian Wittkamp committed
539 540 541 542
        if (MYID==0) mergemod(modfile,3); 
	
	}

Tilman Steinweg's avatar
Tilman Steinweg committed
543
    /* ======================== */
Florian Wittkamp's avatar
Florian Wittkamp committed
544
    /* Read Taper from file     */  
Tilman Steinweg's avatar
Tilman Steinweg committed
545 546
    /* ======================== */
        
Florian Wittkamp's avatar
Florian Wittkamp committed
547
    if((sws>=4)&&(sws<=6)){
548
        
Florian Wittkamp's avatar
Florian Wittkamp committed
549 550 551 552 553 554
          if (MYID==0&&VERBOSE)
          {
                  fprintf(FP,"\n **Message from taper_grid (printed by PE %d):\n",MYID);
                  fprintf(FP," Coefficients for gradient taper are now calculated.\n");
          }
          
Tilman Steinweg's avatar
Tilman Steinweg committed
555 556
        
        if(sws==4){
Florian Wittkamp's avatar
Florian Wittkamp committed
557 558 559
		sprintf(taper_file,"%s",TAPER_FILE_NAME);
		fp_taper=fopen(taper_file,"r");
	}
Tilman Steinweg's avatar
Tilman Steinweg committed
560
        if(sws==5){
Florian Wittkamp's avatar
Florian Wittkamp committed
561 562 563
		sprintf(taper_file,"%s",TAPER_FILE_NAME_U);
		fp_taper=fopen(taper_file,"r");
	}
Tilman Steinweg's avatar
Tilman Steinweg committed
564
        if(sws==6){
Florian Wittkamp's avatar
Florian Wittkamp committed
565 566
		sprintf(taper_file,"%s",TAPER_FILE_NAME_RHO);
		fp_taper=fopen(taper_file,"r");
567 568 569 570
        }
        
        if(fp_taper==NULL) {
            err("Taper file could not opened");
Florian Wittkamp's avatar
Florian Wittkamp committed
571 572
	}

Tilman Steinweg's avatar
Tilman Steinweg committed
573 574 575
        /* loop over global grid */
        for (i=1;i<=NXG;i++){
            for (j=1;j<=NYG;j++){
576 577 578 579
                
                if(feof(fp_taper)){
                    err("Taper file is to small. Check NX*NY of taper file.");
                }
Florian Wittkamp's avatar
Florian Wittkamp committed
580
           
Tilman Steinweg's avatar
Tilman Steinweg committed
581
                fread(&grad_tap, sizeof(float), 1, fp_taper);
Florian Wittkamp's avatar
Florian Wittkamp committed
582
                                                                                       
Tilman Steinweg's avatar
Tilman Steinweg committed
583
                if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
Florian Wittkamp's avatar
Florian Wittkamp committed
584 585 586 587 588
                   ii=i-POS[1]*NX;
                   jj=j-POS[2]*NY;
        
                   taper_coeff[jj][ii]=grad_tap;
                             
Tilman Steinweg's avatar
Tilman Steinweg committed
589 590 591
                }
            }
        }
Florian Wittkamp's avatar
Florian Wittkamp committed
592 593 594 595 596
 
        for (j=1;j<=NY;j++){   
           for (i=1;i<=NX;i++){     
              waveconv[j][i]*=taper_coeff[j][i];
           }                            
Tilman Steinweg's avatar
Tilman Steinweg committed
597
        }   
Florian Wittkamp's avatar
Florian Wittkamp committed
598
                                         
Tilman Steinweg's avatar
Tilman Steinweg committed
599 600 601 602
        fclose(fp_taper);
        
        sprintf(modfile,"taper_coeff_file.bin");
        writemod(modfile,taper_coeff,3);
Florian Wittkamp's avatar
Florian Wittkamp committed
603
                        
Tilman Steinweg's avatar
Tilman Steinweg committed
604 605
        MPI_Barrier(MPI_COMM_WORLD);
        if (MYID==0) mergemod(modfile,3);       
Florian Wittkamp's avatar
Florian Wittkamp committed
606
        } 
Tilman Steinweg's avatar
Tilman Steinweg committed
607 608 609 610
}