[cig-commits] r14356 - in mc/3D/CitcomS/trunk: CitcomS/Components lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Mon Mar 16 16:17:31 PDT 2009
Author: tan2
Date: 2009-03-16 16:17:30 -0700 (Mon, 16 Mar 2009)
New Revision: 14356
Added:
mc/3D/CitcomS/trunk/lib/Seismic_model.c
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Param.py
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Makefile.am
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Output_gzdir.c
mc/3D/CitcomS/trunk/lib/Output_h5.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/lib/output.h
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
Converting temperature/composition fields to seismic velocities, in a format that specfem3d portal and v4.1 can understand.
* solver.param.mineral_physics_mode: default to 3, the model of Tramper, Vacher and Vlaar's PEPI 2001.
* output_optional=seismic will write coordinates and seismic velocities in binary files for specfem comsumption.
* domain bounds file "datafile.domain": binary file with nproc * 10 doubles. The 10 doubles are (rmin, rmax) and four (theta, phi) pairs for the four bottom corner nodes. This file is written by rank-0 processor only.
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Param.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Param.py 2009-03-16 23:14:30 UTC (rev 14355)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Param.py 2009-03-16 23:17:30 UTC (rev 14356)
@@ -53,6 +53,8 @@
reference_state = pyre.inventory.int("reference_state", default=1)
refstate_file = pyre.inventory.str("refstate_file", default="refstate.dat")
+ mineral_physics_model = pyre.inventory.int("mineral_physics_model", default=3)
+
file_vbcs = pyre.inventory.bool("file_vbcs", default=False)
vel_bound_file = pyre.inventory.str("vel_bound_file", default="bvel.dat")
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2009-03-16 23:14:30 UTC (rev 14355)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2009-03-16 23:17:30 UTC (rev 14356)
@@ -525,6 +525,8 @@
input_string("refstate_file",E->refstate.filename,"refstate.dat",m);
}
+ input_int("mineral_physics_model",&(E->control.mineral_physics_model),"1",m);
+
input_int("mat_control",&(E->control.mat_control),"0",m);
input_string("mat_file",E->control.mat_file,"",m);
@@ -1563,6 +1565,16 @@
}
else if(strcmp(prev, "horiz_avg")==0)
E->output.horiz_avg = 1;
+ else if(strcmp(prev, "seismic")==0) {
+ E->output.seismic = E->output.coord_bin = 1;
+
+ /* Total temperature contrast is important when computing seismic velocity,
+ * but it is derived from several parameters. Output it clearly. */
+ if(E->parallel.me==0) {
+ fprintf(stderr, "Total temperature contrast = %f K\n", E->data.ref_temperature);
+ fprintf(E->fp, "Total temperature contrast = %f K\n", E->data.ref_temperature);
+ }
+ }
else if(strcmp(prev, "tracer")==0)
E->output.tracer = 1;
else if(strcmp(prev, "comp_el")==0)
Modified: mc/3D/CitcomS/trunk/lib/Makefile.am
===================================================================
--- mc/3D/CitcomS/trunk/lib/Makefile.am 2009-03-16 23:14:30 UTC (rev 14355)
+++ mc/3D/CitcomS/trunk/lib/Makefile.am 2009-03-16 23:17:30 UTC (rev 14356)
@@ -104,6 +104,7 @@
phase_change.h \
Problem_related.c \
Process_buoyancy.c \
+ Seismic_model.c \
Shape_functions.c \
Size_does_matter.c \
Solver_conj_grad.c \
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2009-03-16 23:14:30 UTC (rev 14355)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2009-03-16 23:17:30 UTC (rev 14356)
@@ -30,7 +30,9 @@
#include <stdlib.h>
+#include <stdint.h>
#include <math.h>
+#include <mpi.h>
#include "element_definitions.h"
#include "global_defs.h"
#include "parsing.h"
@@ -93,7 +95,11 @@
if (cycles == 0) {
output_coord(E);
+ output_domain(E);
/*output_mat(E);*/
+
+ if (E->output.coord_bin)
+ output_coord_bin(E);
}
@@ -120,6 +126,9 @@
if (E->output.horiz_avg)
output_horiz_avg(E, cycles);
+ if (E->output.seismic)
+ output_seismic(E, cycles);
+
if(E->output.tracer && E->control.tracer)
output_tracer(E, cycles);
@@ -177,6 +186,109 @@
}
+void output_domain(struct All_variables *E)
+{
+ /* This routine outputs the domain bounds in a single file. */
+ /* The file will be useful for external program to understand */
+ /* how the CitcomS mesh is domain decomposed. */
+
+ /* Note: rank-0 writes the domain bounds of all processors */
+
+ const int j = 1;
+ const int tag = 0;
+ const int receiver = 0;
+ const int nox = E->lmesh.nox;
+ const int noy = E->lmesh.noy;
+ const int noz = E->lmesh.noz;
+ const int corner_nodes[4] = {1,
+ 1 + noz*(nox-1),
+ nox*noy*noz - (noz -1),
+ 1 + noz*nox*(noy-1)};
+ /* Each line has so many columns:
+ * The columns are min(r) and max(r),
+ * then (theta, phi) of 4 bottom corners. */
+#define ncolumns 10
+
+ double buffer[ncolumns];
+
+ buffer[0] = E->sx[j][3][1];
+ buffer[1] = E->sx[j][3][noz];
+ buffer[2] = E->sx[j][1][corner_nodes[0]];
+ buffer[3] = E->sx[j][2][corner_nodes[0]];
+ buffer[4] = E->sx[j][1][corner_nodes[1]];
+ buffer[5] = E->sx[j][2][corner_nodes[1]];
+ buffer[6] = E->sx[j][1][corner_nodes[2]];
+ buffer[7] = E->sx[j][2][corner_nodes[2]];
+ buffer[8] = E->sx[j][1][corner_nodes[3]];
+ buffer[9] = E->sx[j][2][corner_nodes[3]];
+
+ if(E->parallel.me == 0) {
+ int i, rank;
+ char output_file[255];
+ FILE *fp1;
+ int32_t header[4];
+ MPI_Status status;
+
+ sprintf(output_file,"%s.domain",E->control.data_file);
+ fp1 = output_open(output_file, "wb");
+
+ /* header */
+ header[0] = E->parallel.nproc;
+ header[1] = ncolumns;
+ header[2] = 0x12345678; /* guard */
+ header[3] = sizeof(int32_t);
+ fwrite(header, sizeof(int32_t), 4, fp1);
+
+ /* bounds of self */
+ fwrite(buffer, sizeof(double), ncolumns, fp1);
+
+ /* bounds of other processors */
+ for(rank=1; rank<E->parallel.nproc; rank++) {
+ MPI_Recv(buffer, ncolumns, MPI_DOUBLE, rank, tag, E->parallel.world, &status);
+ fwrite(buffer, sizeof(double), ncolumns, fp1);
+ }
+
+ fclose(fp1);
+ }
+ else {
+ MPI_Send(buffer, ncolumns, MPI_DOUBLE, receiver, tag, E->parallel.world);
+ }
+
+#undef ncolumns
+
+ return;
+}
+
+
+/* write coordinates in binary double */
+void output_coord_bin(struct All_variables *E)
+{
+ int i, j;
+ char output_file[255];
+ FILE *fp1;
+
+ sprintf(output_file,"%s.coord_bin.%d",E->control.data_file,E->parallel.me);
+ fp1 = output_open(output_file, "wb");
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ int32_t header[4];
+ header[0] = E->lmesh.nox;
+ header[1] = E->lmesh.noy;
+ header[2] = E->lmesh.noz;
+ header[3] = 0x12345678; /* guard */
+ fwrite(header, sizeof(int32_t), 4, fp1);
+
+ fwrite(&(E->x[j][1][1]), sizeof(double), E->lmesh.nno, fp1);
+ fwrite(&(E->x[j][2][1]), sizeof(double), E->lmesh.nno, fp1);
+ fwrite(&(E->x[j][3][1]), sizeof(double), E->lmesh.nno, fp1);
+ }
+
+ fclose(fp1);
+
+ return;
+}
+
+
void output_visc(struct All_variables *E, int cycles)
{
int i, j;
@@ -216,7 +328,7 @@
for(j=1;j<=E->sphere.caps_per_proc;j++) {
fprintf(fp1,"%3d %7d\n",j,E->lmesh.nno);
for(i=1;i<=E->lmesh.nno;i++) {
- fprintf(fp1,"%.6e %.6e %.6e %.6e\n",E->sphere.cap[j].V[1][i],E->sphere.cap[j].V[2][i],E->sphere.cap[j].V[3][i],E->T[j][i]);
+ fprintf(fp1,"%.6e %.6e %.6e %.6e %.6e\n",E->sphere.cap[j].V[1][i],E->sphere.cap[j].V[2][i],E->sphere.cap[j].V[3][i],E->T[j][i],E->buoyancy[j][i]);
}
}
@@ -404,6 +516,46 @@
+void output_seismic(struct All_variables *E, int cycles)
+{
+ void compute_horiz_avg();
+ void compute_seismic_model(const struct All_variables*, double*, double*, double*);
+
+ char output_file[255];
+ FILE* fp;
+ int i, j;
+
+ double *rho, *vp, *vs;
+ const int len = 3 * E->lmesh.nno;
+
+ rho = malloc(len * sizeof(double));
+ vp = malloc(len * sizeof(double));
+ vs = malloc(len * sizeof(double));
+ if(rho==NULL || vp==NULL || vs==NULL) {
+ fprintf(stderr, "Error while allocating memory\n");
+ abort();
+ }
+
+ /* isotropic seismic velocity only */
+ /* XXX: update for anisotropy in the future */
+ compute_seismic_model(E, rho, vp, vs);
+
+ sprintf(output_file,"%s.seismic.%d.%d", E->control.data_file, E->parallel.me, cycles);
+ fp = output_open(output_file, "wb");
+
+ fwrite(rho, sizeof(double), E->lmesh.nno, fp);
+ fwrite(vp, sizeof(double), E->lmesh.nno, fp);
+ fwrite(vs, sizeof(double), E->lmesh.nno, fp);
+
+ fclose(fp);
+
+ free(rho);
+ free(vp);
+ free(vs);
+ return;
+}
+
+
void output_mat(struct All_variables *E)
{
int m, el;
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2009-03-16 23:14:30 UTC (rev 14355)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2009-03-16 23:17:30 UTC (rev 14356)
@@ -137,7 +137,11 @@
/* initial I/O */
gzdir_output_coord(E);
+ output_domain(E);
/*gzdir_output_mat(E);*/
+
+ if (E->output.coord_bin)
+ output_coord_bin(E);
}
/*
Modified: mc/3D/CitcomS/trunk/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_h5.c 2009-03-16 23:14:30 UTC (rev 14355)
+++ mc/3D/CitcomS/trunk/lib/Output_h5.c 2009-03-16 23:17:30 UTC (rev 14356)
@@ -228,6 +228,10 @@
#else
if (cycles == 0) {
h5output_const(E);
+ output_domain(E);
+
+ if (E->output.coord_bin)
+ output_coord_bin(E);
}
h5output_timedep(E, cycles);
#endif
Added: mc/3D/CitcomS/trunk/lib/Seismic_model.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Seismic_model.c (rev 0)
+++ mc/3D/CitcomS/trunk/lib/Seismic_model.c 2009-03-16 23:17:30 UTC (rev 14356)
@@ -0,0 +1,280 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<LicenseText>
+ *
+ * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
+ * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
+ * Copyright (C) 1994-2005, California Institute of Technology.
+ *
+ * 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+#include <stdio.h>
+#include "global_defs.h"
+
+
+void compute_horiz_avg(struct All_variables*);
+
+
+/*
+ * Given a radius (non-dimensional),
+ * returns Vp/Vs/rho (in km/s and g/cm^3) of PREM.
+ */
+static void get_prem(double r, double *vp, double *vs, double *rho)
+{
+#define NUM_PREM_LAYERS 11
+
+ /* radius of various layers */
+ const double r_cmb = 0.54622508;
+
+ const double prem_radius[NUM_PREM_LAYERS] =
+ {0.19164966, /* ICB */
+ 0.54622508, /* CMB */
+ 0.56976927, /* top of D'' */
+ 0.87898289, /* 771 */
+ 0.89483598, /* 670 */
+ 0.90582326, /* 600 */
+ 0.93721551, /* 400 */
+ 0.96546853, /* 220 */
+ 0.99617015, /* Moho */
+ 0.99764558, /* middle crust */
+ 1.00000000}; /* top surface */
+
+ /* polynomial coefficients of PREM */
+ const double prem_vs[NUM_PREM_LAYERS][4] =
+ {{ 3.6678, 0.0000, -4.4475, 0.0000},
+ { 0.0010, 0.0000, 0.0000, 0.0000},
+ { 6.9254, 1.4672, -2.0834, 0.9783},
+ {11.1671, -13.7818, 17.4575, -9.2777},
+ {22.3459, -17.2473, -2.0834, 0.9783},
+ { 9.9839, -4.9324, 0.0000, 0.0000},
+ {22.3512, -18.5856, 0.0000, 0.0000},
+ { 8.9496, -4.4597, 0.0000, 0.0000},
+ { 2.1519, 2.3481, 0.0000, 0.0000},
+ { 3.9000, 0.0000, 0.0000, 0.0000},
+ { 3.2000, 0.0000, 0.0000, 0.0000}};
+
+ const double prem_vp[NUM_PREM_LAYERS][4] =
+ {{11.2622, 0.0000, -6.3640, 0.0000},
+ {11.0487, -4.0362, 4.8023,-13.5732},
+ {15.3891, -5.3181, 5.5242, -2.5514},
+ {24.9520, -40.4673, 51.4832,-26.6419},
+ {29.2766, -23.6027, 5.5242, -2.5514},
+ {19.0957, -9.8672, 0.0000, 0.0000},
+ {39.7027, -32.6166, 0.0000, 0.0000},
+ {20.3926, -12.2569, 0.0000, 0.0000},
+ { 4.1875, 3.9382, 0.0000, 0.0000},
+ { 6.8000, 0.0000, 0.0000, 0.0000},
+ { 5.8000, 0.0000, 0.0000, 0.0000}};
+
+ const double prem_rho[NUM_PREM_LAYERS][4] =
+ {{13.0885, 0.0000, -8.8381, 0.0000},
+ {12.5815, -1.2638, -3.6426, -5.5281},
+ { 7.9565, -6.4761, 5.5283, -3.0807},
+ { 7.9565, -6.4761, 5.5283, -3.0807},
+ { 7.9565, -6.4761, 5.5283, -3.0807},
+ { 5.3197, -1.4836, 0.0000, 0.0000},
+ {11.2494, -8.0298, 0.0000, 0.0000},
+ { 7.1089, -3.8045, 0.0000, 0.0000},
+ { 2.6910, 0.6924, 0.0000, 0.0000},
+ { 2.9000, 0.0000, 0.0000, 0.0000},
+ { 2.6000, 0.0000, 0.0000, 0.0000}};
+
+
+ int j;
+ double r2, r3;
+
+ /* make sure r is above CMB */
+ r = (r < r_cmb) ? r_cmb : r;
+ r2 = r * r;
+ r3 = r2 * r;
+
+
+ /* find layer */
+ for (j = 0; j < NUM_PREM_LAYERS; ++j)
+ if (r < prem_radius[j]) break;
+
+ if (j < 0) j = 0;
+ if (j >= NUM_PREM_LAYERS) j = NUM_PREM_LAYERS - 1;
+
+ /* expand polynomials */
+ *vp = prem_vp[j][0]
+ + prem_vp[j][1] * r
+ + prem_vp[j][2] * r2
+ + prem_vp[j][3] * r3;
+ *vs = prem_vs[j][0]
+ + prem_vs[j][1] * r
+ + prem_vs[j][2] * r2
+ + prem_vs[j][3] * r3;
+ *rho = prem_rho[j][0]
+ + prem_rho[j][1] * r
+ + prem_rho[j][2] * r2
+ + prem_rho[j][3] * r3;
+
+ /** debug **
+ fprintf(stderr, "%e %d %f %f %f\n", r, j, *rho, *vp, *vs);
+ /**/
+
+#undef NUM_PREM_LAYERS
+
+ return;
+}
+
+
+
+static void modified_Tramper_Vacher_Vlaar_PEPI2001(struct All_variables *E,
+ double *rho, double *vp, double *vs)
+{
+
+ /* Table 2 in the paper, including quasi-harmonic and anelastic parts */
+ const double dlnvpdt[3] = {-5.71e-5, 2.44e-8, -3.84e-12};
+ const double dlnvsdt[3] = {-9.37e-5, 3.70e-8, -5.46e-12};
+ const double dlnvpdc[3] = {0.172, -0.098, 0.144};
+ const double dlnvsdc[3] = {0.150, -0.143, 0.192};
+
+ const int m = 1;
+ int i, j, nz;
+ double *rhor, *vpr, *vsr, *depthkm;
+ double d, d2, dT, dC, drho, dvp, dvs;
+
+ /* compute horizontal average, if not done yet */
+ if(!E->output.horiz_avg)
+ compute_horiz_avg(E);
+
+
+ /* reference model (PREM) */
+ rhor = malloc((E->lmesh.noz+1) * sizeof(double));
+ vpr = malloc((E->lmesh.noz+1) * sizeof(double));
+ vsr = malloc((E->lmesh.noz+1) * sizeof(double));
+ depthkm = malloc((E->lmesh.noz+1) * sizeof(double));
+
+ for(nz=1; nz<=E->lmesh.noz; nz++) {
+ get_prem(E->sx[m][3][nz], &vpr[nz], &vsr[nz], &rhor[nz]);
+ depthkm[nz] = (E->sphere.ro - E->sx[m][3][nz]) * E->data.radius_km;
+ }
+
+ /* deviation from reference */
+ dC = 0;
+ for(i=0; i<E->lmesh.nno; i++) {
+ nz = (i % E->lmesh.noz) + 1;
+
+ d = depthkm[nz];
+ d2 = d * d;
+ dT = (E->T[m][i+1] - E->Have.T[nz]) * E->data.ref_temperature;
+
+ drho = -dT * E->refstate.thermal_expansivity[nz] * E->data.therm_exp;
+
+ dvp = dT * (dlnvpdt[0] + dlnvpdt[1]*d + dlnvpdt[2]*d2);
+ dvs = dT * (dlnvsdt[0] + dlnvsdt[1]*d + dlnvsdt[2]*d2);
+
+ if(E->control.tracer && E->composition.on && E->composition.ichemical_buoyancy)
+ for(j=0; j<E->composition.ncomp; j++) {
+ dC = E->composition.comp_node[m][j][i+1] - E->Have.C[j][nz];
+
+ drho += dC * E->composition.buoyancy_ratio[j]
+ * E->data.ref_temperature * E->data.therm_exp;
+
+ dvp += dC * (dlnvpdc[0] + dlnvpdc[1]*d + dlnvpdc[2]*d2);
+ dvs += dC * (dlnvsdc[0] + dlnvsdc[1]*d + dlnvsdc[2]*d2);
+ }
+
+ rho[i] = rhor[nz] * (1 + drho);
+ vp[i] = vpr[nz] * (1 + dvp);
+ vs[i] = vsr[nz] * (1 + dvs);
+
+ /** debug **
+ fprintf(stderr, "node=%d dT=%f K, dC=%f, %e %e %e\n",
+ i, dT, dC, drho, dvp, dvs);
+ /**/
+ }
+
+
+ free(rhor);
+ free(vpr);
+ free(vsr);
+ free(depthkm);
+
+ return;
+}
+
+
+
+void compute_seismic_model(struct All_variables *E,
+ double *rho, double *vp, double *vs)
+{
+
+ switch(E->control.mineral_physics_model) {
+
+ case 0:
+ /* reserved for Stixrude and Lithgow-Bertelloni, GJI, 2005 model */
+ fprintf(stderr,"Invalid value: 'mineral_physics_model=%d'\n",
+ E->control.mineral_physics_model);
+ parallel_process_termination();
+
+ break;
+
+ case 1:
+ /* reserved for Karato, GRL, 1993 model */
+ fprintf(stderr,"Invalid value: 'mineral_physics_model=%d'\n",
+ E->control.mineral_physics_model);
+ parallel_process_termination();
+
+ break;
+
+ case 2:
+ /* reserved for Stacy, PEPI, 1998 model */
+ fprintf(stderr,"Invalid value: 'mineral_physics_model=%d'\n",
+ E->control.mineral_physics_model);
+ parallel_process_termination();
+
+ break;
+
+ case 3:
+ /* Based on the paper:
+ * Trampert, Vacher, and Vlaar, PEPI, 2001.
+ *
+ * Note that the paper has its own reference profile (which is not
+ * shown in the paper), and is only valid between
+ * 1000 km < depth < 2600 km.
+ *
+ * But here we use PREM as the reference model, and extend the model to the
+ * whole mantle.
+ */
+
+ modified_Tramper_Vacher_Vlaar_PEPI2001(E, rho, vp, vs);
+ break;
+
+ case 100:
+ /* user-defined initial temperature goes here */
+ fprintf(stderr,"Need user definition for mineral physics model: 'mineral_physics_model=%d'\n",
+ E->control.mineral_physics_model);
+ parallel_process_termination();
+ break;
+
+ default:
+ /* unknown option */
+ fprintf(stderr,"Invalid value: 'mineral_physics_model=%d'\n",
+ E->control.mineral_physics_model);
+ parallel_process_termination();
+ break;
+ }
+
+ return;
+}
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2009-03-16 23:14:30 UTC (rev 14355)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2009-03-16 23:17:30 UTC (rev 14356)
@@ -518,6 +518,7 @@
int vbcs_file;
int tbcs_file;
int mat_control;
+ int mineral_physics_model;
#ifdef USE_GGRD
struct ggrd_master ggrd;
float *surface_rayleigh;
@@ -635,6 +636,8 @@
int botm; /* whether to output bottom data */
int geoid; /* whether to output geoid/topo spherial harmonics */
int horiz_avg; /* whether to output horizontal averaged profile */
+ int seismic; /* whether to output seismic velocity model */
+ int coord_bin; /* whether to output coordinates in binary format */
int tracer; /* whether to output tracer coordinate */
int comp_el; /* whether to output composition at elements */
int comp_nd; /* whether to output composition at nodes */
Modified: mc/3D/CitcomS/trunk/lib/output.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/output.h 2009-03-16 23:14:30 UTC (rev 14355)
+++ mc/3D/CitcomS/trunk/lib/output.h 2009-03-16 23:17:30 UTC (rev 14356)
@@ -37,6 +37,9 @@
void output(struct All_variables *, int);
void output_time(struct All_variables *, int);
void output_checkpoint(struct All_variables *);
+void output_coord_bin(struct All_variables *);
+void output_domain(struct All_variables *);
+void output_seismic(struct All_variables *, int);
FILE* output_open(char *, char *);
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2009-03-16 23:14:30 UTC (rev 14355)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2009-03-16 23:17:30 UTC (rev 14356)
@@ -363,6 +363,8 @@
getStringProperty(properties, "refstate_file", E->refstate.filename, fp);
}
+ getIntProperty(properties, "mineral_physics_model", E->control.mineral_physics_model, fp);
+
getIntProperty(properties, "file_vbcs", E->control.vbcs_file, fp);
getStringProperty(properties, "vel_bound_file", E->control.velocity_boundary_file, fp);
More information about the CIG-COMMITS
mailing list