readmod_viscac.c 4.3 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
 /*-----------------------------------------------------------------------------------------
 * Copyright (C) 2013  For the list of authors, see file AUTHORS.
 *
 * 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>.
-----------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
 *   Read viscoacoustic model properties (vp,density,Qp) from files  
 *
 *  ----------------------------------------------------------------------*/


/* This file contains function readmod, which has the purpose
   to read data from model-files for viscoelastic simulation */

#include "fd.h"

void readmod_viscac(float  **  rho, float **  pi, float ** taup, float * eta){
	
	extern float DT, *FL, TAU;
	extern int L;
34
	extern int NX, NY, NXG, NYG,  POS[3], MYID, PARAMETERIZATION;
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
	extern char  MFILE[STRING_SIZE];	
	extern FILE *FP;
	
	
	/* local variables */
	float rhov, piv, vp, qp, *pts;
	int i, j, ii, jj, l, sw_Qp=1;
	FILE *fp_vp, *fp_rho, *fp_qp;
	char filename[STRING_SIZE];
	
	
	/* 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];
	}
	
53
	fprintf(FP,"\n...reading model information from model-files...\n");
54 55 56
	
	/* read density and seismic velocities */
	/* ----------------------------------- */
57
	if(PARAMETERIZATION==1){
58 59 60
		fprintf(FP,"\t Vp:\n\t %s.vp\n\n",MFILE);
		sprintf(filename,"%s.vp",MFILE);
		fp_vp=fopen(filename,"r");
61
		if (fp_vp==NULL) declare_error(" Could not open model file for Vp ! ");
62 63 64 65
		
		fprintf(FP,"\t Density:\n\t %s.rho\n\n",MFILE);
		sprintf(filename,"%s.rho",MFILE);
		fp_rho=fopen(filename,"r");
66
		if (fp_rho==NULL) declare_error(" Could not open model file for densities ! ");
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
		
		fprintf(FP,"\t Qp:\n\t %s.qp\n\n",MFILE);
		sprintf(filename,"%s.qp",MFILE);
		fp_qp=fopen(filename,"r");
		if (fp_qp==NULL){
			if (MYID==0){
				printf(" Could not open model file for Qp-values ! \n");
				printf(" Uses input file value for tau! \n");
			}
			sw_Qp=0;
		}
	}
	
	
	/* loop over global grid */
	for (i=1;i<=NXG;i++){
	for (j=1;j<=NYG;j++){
84 85
        
        if(feof(fp_vp) && feof(fp_rho)){
86
            declare_error("Model file VP or RHO is to small. Check dimensions NX*NY of file.");
87 88 89
        }
		
        fread(&vp, sizeof(float), 1, fp_vp);
90
		fread(&rhov, sizeof(float), 1, fp_rho);
91 92 93 94
		
        if (sw_Qp){
            
            if(feof(fp_qp)){
95
                declare_error("Model file QP is to small. Check dimensions NX*NY of file.");
96 97
            }
            
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
			fread(&qp, sizeof(float), 1, fp_qp);
		}
		
		
		/* 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;
		
			rho[jj][ii]=rhov;
			pi[jj][ii]=vp;
			
			if (sw_Qp){
				taup[jj][ii]=2.0/qp;}
			else taup[jj][ii]=TAU;
		}
	}
	}
	
118 119
    fread(&vp, sizeof(float), 1, fp_vp);
    if(!feof(fp_vp)){
120
        declare_error("Model file VP is to big. Check dimensions NX*NY of file.");
121 122 123 124 125
    }
    fclose(fp_vp);
    
    fread(&rho, sizeof(float), 1, fp_rho);
    if(!feof(fp_rho)){
126
        declare_error("Model file RHO is to big. Check dimensions NX*NY of file.");
127 128
    }
    fclose(fp_rho);
129

130 131 132
    if (sw_Qp){
        fread(&qp, sizeof(float), 1, fp_qp);
        if(!feof(fp_qp)){
133
            declare_error("Model file QP is to big. Check dimensions NX*NY of file.");
134 135 136 137
        }
        fclose(fp_qp);
    }
    
138
	
139 140 141 142 143 144 145 146 147 148
	/* write model to disk */
	sprintf(filename,"%s.out.vp",MFILE);
    write_matrix_disk(pi,filename);
    
    sprintf(filename,"%s.out.rho",MFILE);
    write_matrix_disk(rho,filename);
    
    sprintf(filename,"%s.out.qp",MFILE);
    write_matrix_disk(taup,filename);
    
149 150
	free_vector(pts,1,L);
}