[cig-commits] r14034 - in seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D: . citcoms_prem_isotropic_no_crust
tan2 at geodynamics.org
tan2 at geodynamics.org
Tue Feb 10 15:21:26 PST 2009
Author: tan2
Date: 2009-02-10 15:21:26 -0800 (Tue, 10 Feb 2009)
New Revision: 14034
Added:
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/build.mk
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/citcoms_prem_isotropic_no_crust.c
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/citcoms_prem_isotropic_no_crust.f90
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/description.txt
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/prem
Log:
First step for citcoms earth model. This model is based on PREM, isotropic, and no crust.
Read citcoms domain bound. Need help from Leif to figure out specfem domain bound.
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/build.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/build.mk (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/build.mk 2009-02-10 23:21:26 UTC (rev 14034)
@@ -0,0 +1,9 @@
+
+include $S/MODELS/1D_ref/prem/build.mk
+
+mantle_model_OBJECTS = \
+ $O/mantle_model.o \
+ $(EMPTY_MACRO)
+
+$O/mantle_model.o: constants.h $S/MODELS/3D_mantle/citcoms_prem_isotropic_no_crust/citcoms_prem_isotropic_no_crust.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/mantle_model.o ${FCFLAGS_f90} $S/MODELS/3D_mantle/citcoms_prem_isotropic_no_crust/citcoms_prem_isotropic_no_crust.f90
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/citcoms_prem_isotropic_no_crust.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/citcoms_prem_isotropic_no_crust.c (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/citcoms_prem_isotropic_no_crust.c 2009-02-10 23:21:26 UTC (rev 14034)
@@ -0,0 +1,274 @@
+
+#include "config.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include "mpi.h"
+
+
+static int read_long_vector(FILE*, int, long*);
+static int read_double_vector(FILE*, int, double*);
+
+struct citcoms_domain {
+ int rank;
+
+ double corners[8];
+ double rmin, rmax;
+
+ struct citcoms_domain *next;
+};
+
+typedef struct citcoms_domain citcoms_domain_t;
+
+
+/* model static variable declaration */
+citcoms_domain_t *first;
+
+/* citcoms local mesh geometry */
+int nox, noy, noz;
+
+
+void FC_FUNC_(read_3d_mantle_model_in_c,READ_3D_MANTLE_MODEL_IN_C)
+ (double theta0, double phi0, double theta1, double phi1,
+ double theta2, double phi2, double theta3, double phi3,
+ char citcoms_bound_filename[], char citcoms_model_filename_base,
+ int len_citcoms_bound_filename, int len_citcoms_model_filename_base) {
+
+ const int ncolumns = 11;
+ int rank;
+ int i, buffer_size;
+ int total_citcoms_proc;
+ double corners[8];
+ double rmin, rmax;
+ double *buffer;
+ FILE *fp;
+ citcoms_domain_t *curr, *prev;
+
+ /***********************************************************************
+ * Read citcoms domain bound. Rank-0 processor reads the file and
+ * broadcast to other processors.
+ * See citcoms manual for file format
+ ***********************************************************************/
+
+ MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+ if(rank == 0) {
+ long ltmp[2];
+
+ fp = fopen(citcoms_bound_filename, "r");
+
+ /* read header line */
+ if(read_long_vector(fp, 2, ltmp) != 2) {
+ fprintf(stderr, "error while reading first line of '%s'\n",
+ citcoms_bound_filename);
+ abort();
+ }
+
+ if(ltmp[1] != ncolumns) {
+ fprintf(stderr, "'%s' contains wrong number of columns\n",
+ citcoms_bound_filename);
+ abort();
+ }
+
+ total_citcoms_proc = (int) ltmp[0];
+ }
+
+ MPI_Bcast(&total_citcoms_proc, 1, MPI_INT, 0, MPI_COMM_WORLD);
+
+ buffer_size = total_citcoms_proc * ncolumns;
+ buffer = malloc(buffer_size * sizeof(double));
+ if(buffer == NULL) {
+ fprintf(stderr, "error while allocating memory '%s' 1\n",
+ __FILE__);
+ abort();
+ }
+
+ if(rank == 0) {
+ /* read content of domain bound */
+ for(i=0; i<total_citcoms_proc; i++) {
+ if(read_double_vector(fp, ncolumns, &(buffer[i*ncolumns])) != ncolumns) {
+ fprintf(stderr, "error while reading %d-th line of '%s'\n",
+ i, citcoms_bound_filename);
+ abort();
+ }
+ }
+ fclose(fp);
+ }
+
+ MPI_Bcast(buffer, buffer_size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+
+ /***********************************************************************
+ * Find out which citcoms domain overlaps with my domain.
+ ***********************************************************************/
+
+ first = malloc(sizeof(citcoms_domain_t));
+ curr = first;
+ prev = NULL;
+
+ for(i=0; i<total_citcoms_proc; i++) {
+ double *tmp = &(buffer[i*ncolumns]);
+ //XXX: if buffer[i*ncolumns] overlapped with my domain
+ {
+ curr->rank = (int) tmp[0];
+ memcpy(curr->corners, &(tmp[1]), 8);
+ curr->rmin = tmp[9];
+ curr->rmax = tmp[10];
+
+ curr->next = malloc(sizeof(citcoms_domain_t));
+ prev = curr;
+ curr = curr->next;
+ }
+ }
+
+ if(prev == NULL) {
+ fprintf(stderr, "cannot find overllaping region '%s' 2\n",
+ __FILE__);
+ abort();
+ }
+
+ /* this pointer is over-allocated */
+ free(curr);
+ prev->next = NULL;
+
+ free(buffer);
+
+ /***********************************************************************
+ * Read the mantle model from those citcoms domain
+ ***********************************************************************/
+
+ curr = first;
+ do {
+ const char fmt[] = "%s.%d";
+ char filename[256];
+
+ snprintf(filename, 255, fmt, curr->rank);
+
+ fp = fopen(filename, "rb");
+
+ // XXX: read and store
+
+ fclose(fp);
+ } while (curr = curr->next)
+
+
+ return;
+}
+
+
+void FC_FUNC_(iso_mantle_model,ISO_MANTLE_MODEL)
+ (double *pRadius, double *pTheta, double *pPhi,
+ double *vpv, double *vph, double *vsv, double *vsh,
+ double *rho, double *eta_aniso)
+{
+ /* Define 3D isotropic models here. */
+
+ double radius, theta, phi;
+
+ /* inputs */
+ radius = *pRadius;
+ theta = *pTheta;
+ phi = *pPhi;
+
+ /* outputs */
+ *vpv = 0.0;
+ *vph = 0.0;
+ *vsv = 0.0;
+ *vsh = 0.0;
+ *rho = 0.0;
+ *eta_aniso = 0.0;
+}
+
+
+/*
+ * Scan input str to get a long int vector *values. The vector length is from
+ * input len. The input str contains white-space seperated numbers.
+ */
+static int scan_long_vector(const char *str, int len, long *values)
+{
+ char *nptr, *endptr;
+ int i;
+
+ /* cast to avoid compiler warning */
+ nptr = endptr = (char *) str;
+
+ for (i = 0; i < len; ++i) {
+ values[i] = strtol(nptr, &endptr, 10);
+ if (nptr == endptr) {
+ /* error: no conversion is performed */
+ return -1;
+ }
+ nptr = endptr;
+ }
+
+ /** debug **
+ for (i = 0; i < len; ++i) fprintf(stderr, "%ld, ", values[i]);
+ fprintf(stderr, "\n");
+ /**/
+ return len;
+}
+
+
+/*
+ * From input file, read a line, which contains white-space seperated numbers
+ * of lenght num_columns, store the numbers in a long int array.
+ * Return num_columns on success, 0 on file reading error, and -1 on conversion error.
+ */
+static int read_long_vector(FILE *in, int num_columns, long *fields)
+{
+ char buffer[256], *p;
+
+ p = fgets(buffer, 255, in);
+ if (!p) {
+ return 0;
+ }
+
+ return scan_long_vector(buffer, num_columns, fields);
+}
+
+
+/*
+ * Scan input str to get a double vector *values. The vector length is from
+ * input len. The input str contains white-space seperated numbers.
+ */
+static int scan_double_vector(const char *str, int len, double *values)
+{
+ char *nptr, *endptr;
+ int i;
+
+ /* cast to avoid compiler warning */
+ nptr = endptr = (char *) str;
+
+ for (i = 0; i < len; ++i) {
+ values[i] = strtod(nptr, &endptr);
+ if (nptr == endptr) {
+ /* error: no conversion is performed */
+ return -1;
+ }
+ nptr = endptr;
+ }
+
+ /** debug **
+ for (i = 0; i < len; ++i) fprintf(stderr, "%e, ", values[i]);
+ fprintf(stderr, "\n");
+ /**/
+ return len;
+}
+
+
+/*
+ * From input file, read a line, which contains white-space seperated numbers
+ * of lenght num_columns, store the numbers in a double array.
+ * Return num_columns on success, 0 on file reading error, and -1 on conversion error.
+ */
+static int read_double_vector(FILE *in, int num_columns, double *fields)
+{
+ char buffer[256], *p;
+
+ p = fgets(buffer, 255, in);
+ if (!p) {
+ return 0;
+ }
+
+ return scan_double_vector(buffer, num_columns, fields);
+}
+
+
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/citcoms_prem_isotropic_no_crust.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/citcoms_prem_isotropic_no_crust.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/citcoms_prem_isotropic_no_crust.f90 2009-02-10 23:21:26 UTC (rev 14034)
@@ -0,0 +1,172 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! This program 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; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program 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 this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+
+!---------------------------
+
+ subroutine get_model_properties(HONOR_1D_SPHERICAL_MOHO,ONE_CRUST, &
+ TRANSVERSE_ISOTROPY, &
+ ISOTROPIC_3D_MANTLE,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
+ CRUSTAL,CASE_3D, &
+ ATTENUATION_3D)
+
+ implicit none
+
+ logical HONOR_1D_SPHERICAL_MOHO,ONE_CRUST,&
+ TRANSVERSE_ISOTROPY,&
+ ISOTROPIC_3D_MANTLE,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,&
+ CRUSTAL,CASE_3D,&
+ ATTENUATION_3D
+
+ HONOR_1D_SPHERICAL_MOHO = .false.
+ ONE_CRUST = .true.
+
+ TRANSVERSE_ISOTROPY = .true.
+
+ ISOTROPIC_3D_MANTLE = .true.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+
+ CRUSTAL = .false.
+ CASE_3D = .false.
+
+ ATTENUATION_3D = .false.
+
+ end subroutine get_model_properties
+
+!---------------------------
+
+ subroutine read_3d_mantle_model()
+
+ implicit none
+ character(len=150) citcoms_bound_filename, citcoms_model_filename_base
+
+ ! coordinates of the four corner nodes in the current slice
+ double precision theta0, phi0
+ double precision theta1, phi1
+ double precision theta2, phi2
+ double precision theta3, phi3
+
+ read_value_string(citcoms_bound_filename, 'CITCOMS_BOUND_FILENAME')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+
+ if(len(trim(citcoms_bound_filename)) == 0) &
+ stop 'CITCOMS_BOUND_FILENAME undefined'
+
+ if(len(trim(citcoms_bound_filename)) >= 150-1) &
+ stop 'CITCOMS_BOUND_FILENAME is too long'
+
+ read_value_string(citcoms_model_filename_base, 'CITCOMS_MODEL_FILENAME_BASE')
+ if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
+
+ if(len(trim(citcoms_model_filename_base)) == 0) &
+ stop 'CITCOMS_MODEL_FILENAME_BASE undefined'
+
+ if(len(trim(citcoms_model_filename_base)) >= 150-1) &
+ stop 'CITCOMS_MODEL_FILENAME_BASE is too long'
+
+
+ ! XXX: fill in values of (theta?, phi?)...
+
+
+ call read_3d_mantle_model_in_c(theta0, phi0, theta1, phi1, &
+ theta2, phi2, theta3, phi3, &
+ trim(citcoms_bound_filename)//CHAR(0), &
+ trim(citcoms_model_filename_base)//CHAR(0))
+
+ end subroutine read_3d_mantle_model
+
+!---------------------------
+
+ subroutine iso_mantle_model(radius,theta,phi, &
+ vpv,vph,vsv,vsh,rho,eta_aniso)
+
+ implicit none
+ include "constants.h"
+
+ double precision radius,theta,phi
+ double precision vpv,vph,vsv,vsh,rho,eta_aniso
+
+ double precision dvp,dvs,drho
+
+ end subroutine iso_mantle_model
+
+!---------------------------
+
+ subroutine aniso_mantle_model(r,theta,phi,rho, &
+ c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
+
+ implicit none
+
+ double precision r,theta,phi
+ double precision rho
+ double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26, &
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+
+ ! no-op
+
+ end subroutine aniso_mantle_model
+
+!----------------------------------
+
+ subroutine add_moho_topography(myrank,xelm,yelm,zelm,RMOHO,R220)
+
+ implicit none
+
+ include "constants.h"
+
+ integer myrank
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+
+ double precision RMOHO,R220
+
+ ! no-np
+
+ end subroutine add_moho_topography
+
+!----------------------------------
+
+ subroutine add_topography_410_650(myrank,xelm,yelm,zelm,R220,R400,R670,R771)
+
+ implicit none
+
+ include "constants.h"
+
+ integer myrank
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+
+ double precision R220,R400,R670,R771
+
+ ! no-op
+
+ end subroutine add_topography_410_650
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/description.txt
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/description.txt (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/description.txt 2009-02-10 23:21:26 UTC (rev 14034)
@@ -0,0 +1 @@
+The code reads the (theta, phi, r, drho, dVp, dVs) output from CitcomS, relative to PREM, then interpolates to the GLL points.
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/prem
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/prem (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/prem 2009-02-10 23:21:26 UTC (rev 14034)
@@ -0,0 +1 @@
+link ../../1D_ref/prem/
\ No newline at end of file
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/prem
___________________________________________________________________
Name: svn:special
+ *
More information about the CIG-COMMITS
mailing list