[cig-commits] r8974 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Sun Dec 30 18:25:01 PST 2007
Author: becker
Date: 2007-12-30 18:25:00 -0800 (Sun, 30 Dec 2007)
New Revision: 8974
Modified:
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_init.c
Log:
Added -dsf option for depth-dependent scaling factors.
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2007-12-29 07:02:49 UTC (rev 8973)
+++ mc/1D/hc/trunk/hc.h 2007-12-31 02:25:00 UTC (rev 8974)
@@ -74,7 +74,7 @@
#include "ggrd_grdtrack_util.h"
/*
-spherical haronics
+spherical harmonics
*/
#include "sh.h"
@@ -140,12 +140,14 @@
hc_boolean read_short_dens_sh; /* short SH format for density
files? */
+ hc_boolean read_dens_scale_from_file; /* read the density/velocity scaling from file? */
hc_boolean print_pt_sol; /* output of p[6] and t[2] vectors */
char visc_filename[HC_CHAR_LENGTH]; /* name of viscosity profile file */
char pvel_filename[HC_CHAR_LENGTH]; /* name of plate velocities file */
char dens_filename[HC_CHAR_LENGTH]; /* name of density model file */
char prem_model_filename[HC_CHAR_LENGTH]; /* PREM model filename */
+ char dens_scaling_filename[HC_CHAR_LENGTH]; /* */
};
/*
@@ -186,13 +188,13 @@
HC_PREC *rden; /* radii of density layers [normalized] */
struct sh_lms *dens_anom; /*
- expansions of density
- anomalies has to be [inho] (if
- those change, that's OK, no
- need to call with
- dens_fac_changed == TRUE)
- */
- HC_PREC dens_scale[1]; /* scale for density file */
+ expansions of density
+ anomalies has to be [inho] (if
+ those change, that's OK, no
+ need to call with
+ dens_fac_changed == TRUE)
+ */
+ HC_PREC dens_scale; /* scale for density file */
@@ -407,5 +409,14 @@
*/
#define HC_Z_DEPTH(x) ((HC_RE_KM * (1.0-(x))))
+/*
+
+default constant scaling for input density file, use 0.01 for %
+
+*/
+
+#define HC_DENSITY_SCALING 0.01
+
+
#define __HC_READ_HEADER_FILE__
#endif
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2007-12-29 07:02:49 UTC (rev 8973)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2007-12-31 02:25:00 UTC (rev 8974)
@@ -16,6 +16,8 @@
void ggrd_vecrealloc(double **, int, char *);
float ggrd_gt_rms(float *, int);
float ggrd_gt_mean(float *, int);
+void ggrd_global_bcr_assign(struct BCR *);
+void my_GMT_bcr_init(struct GRD_HEADER *, int *, int, struct BCR *);
/* ggrd_readgrds.c */
void ggrd_init_vstruc(struct ggrd_vel *);
int ggrd_read_vel_grids(struct ggrd_vel *, double, unsigned short, unsigned short, char *);
@@ -26,6 +28,7 @@
int ggrd_find_vel_and_der(double *, double, double, struct ggrd_vel *, int, unsigned short, unsigned short, double *, double *, double *);
void ggrd_get_velocities(double *, double *, double *, int, struct ggrd_vel *, double, double);
void ggrd_weights(double, double *, int, int, double [(5 +1)][(1 +1)]);
+/* grdinttester.c */
/* hc_extract_sh_layer.c */
/* hc_init.c */
void hc_init_parameters(struct hc_parameters *);
@@ -35,7 +38,8 @@
void hc_init_constants(struct hcs *, double, char *, unsigned short);
void hc_handle_command_line(int, char **, struct hc_parameters *);
void hc_assign_viscosity(struct hcs *, int, char [300], unsigned short);
-void hc_assign_density(struct hcs *, unsigned short, int, char *, int, unsigned short, unsigned short, unsigned short, unsigned short);
+void hc_assign_density(struct hcs *, unsigned short, int, char *, int, unsigned short, unsigned short, unsigned short, unsigned short, unsigned short, char *);
+double hc_find_dens_scale(double, double, unsigned short, double *, double *, int);
void hc_init_phase_boundaries(struct hcs *, int, unsigned short);
void hc_assign_plate_velocities(struct hcs *, int, char *, unsigned short, int, unsigned short, unsigned short);
void hc_init_l_factors(struct hcs *, int);
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2007-12-29 07:02:49 UTC (rev 8973)
+++ mc/1D/hc/trunk/hc_init.c 2007-12-31 02:25:00 UTC (rev 8974)
@@ -32,7 +32,8 @@
p->read_short_dens_sh = FALSE; /* read the density anomaly file in
the short format of Becker &
Boschi (2002)? */
-
+ p->read_dens_scale_from_file = FALSE; /* read a depth-dependent scaling from file? */
+
p->verbose = 0; /* debugging output? (0,1,2,3,4...) */
p->sol_binary_out = TRUE; /* binary or ASCII output of SH expansion */
p->solution_mode = HC_VEL; /* default: velocity */
@@ -125,14 +126,10 @@
HC_ERROR("hc_init","vel_bc_zero and free_slip doesn't make sense");
/* read in the densities */
hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,
- p->dens_filename,-1,FALSE,FALSE,
- p->verbose,
- p->read_short_dens_sh);
+ p->dens_filename,-1,FALSE,FALSE,p->verbose,p->read_short_dens_sh,
+ p->read_dens_scale_from_file,p->dens_scaling_filename);
/* assign all zeroes up to the lmax of the density expansion */
- hc_assign_plate_velocities(hc,HC_INIT_FROM_FILE,
- p->pvel_filename,
- p->vel_bc_zero,
- hc->dens_anom[0].lmax,
+ hc_assign_plate_velocities(hc,HC_INIT_FROM_FILE,p->pvel_filename,p->vel_bc_zero,hc->dens_anom[0].lmax,
FALSE,p->verbose);
}else{
/* presribed surface velocities */
@@ -142,17 +139,16 @@
p->pvel_filename,
p->vel_bc_zero,
dummy,FALSE,p->verbose);
- hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,
- p->dens_filename,hc->pvel[0].lmax,
- FALSE,FALSE,p->verbose,
- p->read_short_dens_sh);
+ hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,p->dens_filename,hc->pvel[0].lmax,
+ FALSE,FALSE,p->verbose,p->read_short_dens_sh,
+ p->read_dens_scale_from_file,p->dens_scaling_filename);
}else{
if(p->verbose)
fprintf(stderr,"hc_init: initializing for free-slip\n");
/* read in the density fields */
- hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,
- p->dens_filename,-1,FALSE,FALSE,
- p->verbose,p->read_short_dens_sh);
+ hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,p->dens_filename,-1,FALSE,FALSE,
+ p->verbose,p->read_short_dens_sh,
+ p->read_dens_scale_from_file,p->dens_scaling_filename);
}
}
/*
@@ -190,7 +186,7 @@
/*
density scale
*/
- hc->dens_scale[0] = dens_anom_scale;
+ hc->dens_scale = dens_anom_scale;
/*
constants
*/
@@ -265,20 +261,23 @@
argv[0]);
fprintf(stderr,"This code can compute velocities, tractions, and geoid for simple density distributions\n");
fprintf(stderr,"and plate velocities using the semi-analytical approach of Hager & O'Connell (1981).\n");
- fprintf(stderr,"This particular implementation illustrates one possible way to combine the HC solver routines.\n\n");
- fprintf(stderr,"Based on code by Brad Hager, Richard O'Connell, and Bernhard Steinberger.\n\n");
+ fprintf(stderr,"This particular implementation illustrates one possible way to combine the HC solver routines.\n");
+ fprintf(stderr,"Based on code by Brad Hager, Richard O'Connell, and Bernhard Steinberger.\n");
+ fprintf(stderr,"This version by Thorsten Becker and Craig O'Neill\n\n");
fprintf(stderr,"density anomaly options:\n");
fprintf(stderr,"-dens\tname\tuse name as a SH density anomaly model (%s)\n",
p->dens_filename);
- fprintf(stderr,"\t\tThis expects density anomalies in %% of PREM in Dahlen & Tromp format.\n");
+ fprintf(stderr,"\t\tAll density anomalies are in units of %g%% of PREM, all SH coefficients in Dahlen & Tromp convention.\n",
+ HC_DENSITY_SCALING*100);
fprintf(stderr,"-dshs\t\tuse the short, Becker & Boschi (2002) format for the SH density model (%s)\n",
hc_name_boolean(p->read_short_dens_sh));
- fprintf(stderr,"-ds\t\tadditional density anomaly scaling factor (%g)\n\n",
+ fprintf(stderr,"-ds\t\tadditional density anomaly scaling factor (%g)\n",
p->dens_anom_scale);
-
+ fprintf(stderr,"-dsf\tfile\tread depth dependent density scaling from file (overrides -ds, %s)\n\n",
+ hc_name_boolean(p->read_dens_scale_from_file));
fprintf(stderr,"Earth model options:\n");
fprintf(stderr,"-prem\tname\tset Earth model to name (%s)\n",
p->prem_model_filename);
@@ -338,6 +337,10 @@
}else if(strcmp(argv[i],"-dens")==0){ /* density filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->dens_filename,argv[i],HC_CHAR_LENGTH);
+ }else if(strcmp(argv[i],"-dsf")==0){ /* density scaling filename */
+ hc_toggle_boolean(&p->read_dens_scale_from_file);
+ hc_advance_argument(&i,argc,argv);
+ strncpy(p->dens_scaling_filename,argv[i],HC_CHAR_LENGTH);
}else if(strcmp(argv[i],"-visc")==0){ /* viscosity filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->visc_filename,argv[i],HC_CHAR_LENGTH);
@@ -452,7 +455,7 @@
mean += log(hc->visc[hc->nvis-1]) * mweight;
mean = exp(mean/mws);
- fprintf(stderr,"hc_assign_viscosity: read %i layers of non-dimensionalized viscosities from %s\n",
+ fprintf(stderr,"hc_assign_viscosity: read %i layered viscosity[Pas] from %s\n",
hc->nvis,filename);
fprintf(stderr,"hc_assign_viscosity: rough estimate of mean viscosity: %g x %g = %g Pas\n",
mean, hc->visnor, mean*hc->visnor);
@@ -496,21 +499,16 @@
hc_boolean layer_structure_changed,
hc_boolean density_in_binary,
hc_boolean verbose,
- hc_boolean use_short_format)
+ hc_boolean use_short_format,
+ hc_boolean read_dens_scale_from_file,
+ char *dens_scaling_filename)
{
FILE *in;
- int type,lmax,shps,ilayer,nset,ivec,i,j;
- HC_PREC *dtop,*dbot,zlabel,dens_scale[1],rho0;
- double dtmp;
+ int type,lmax,shps,ilayer,nset,ivec,i,j,ndf;
+ HC_PREC *dtop,*dbot,zlabel,local_scale,dens_scale[0],rho0,*rdf=NULL,*sdf=NULL,smean;
+ double dtmp[2];
hc_boolean reported = FALSE,read_on;
- static HC_PREC local_dens_fac = .01; /* this factor will be
- multiplied with the
- hc->dens_fac factor to
- arrive at fractional
- anomalies from input. set
- to 0.01 for percent input,
- for instance
- */
+
hc->compressible = compressible;
hc->inho = 0;
if(hc->dens_init) /* clear old expansions, if
@@ -523,7 +521,36 @@
case HC_INIT_FROM_FILE:
if(hc->dens_init)
HC_ERROR("hc_assign_density","really read dens anomalies again from file?");
-
+ if(read_dens_scale_from_file){
+ /* depth depending scaling, format same as for viscosity file */
+ if(verbose)
+ fprintf(stderr,"hc_assign_density: reading depth dependent dln\\rho/dln density scaling from %s\n",
+ dens_scaling_filename);
+ ndf=0;smean = 0.0;
+ in = hc_open(dens_scaling_filename,"r","hc_assign_density");
+ while(fscanf(in,"%lf %lf",dtmp,(dtmp+1)) == 2){
+ hc_vecrealloc(&rdf,(1+ndf),"hc_assign_density");
+ hc_vecrealloc(&sdf,(1+ndf),"hc_assign_density");
+ rdf[ndf] = dtmp[0];sdf[ndf] = dtmp[1];
+ smean+=sdf[ndf];
+ ndf++;
+ }
+ fclose(in);
+ if(!ndf){
+ fprintf(stderr,"hc_assign_density: error: did not read any density scaling factors from %s\n",
+ dens_scaling_filename);
+ exit(-1);
+ }
+ smean /= (HC_PREC)ndf;
+ if(verbose)
+ fprintf(stderr,"hc_assign_density: read scaling on %i layers, rough mean: %g\n",ndf,smean);
+ /* end init */
+ }else{
+ if(verbose)
+ fprintf(stderr,"hc_assign_density: using constant dln\\rho/dln input density scaling of %g\n",
+ hc->dens_scale);
+ ndf = 0;
+ }
/*
read in density anomalies in spherical harmonics format for
@@ -535,8 +562,8 @@
in = hc_open(filename,"r","hc_assign_density");
if(verbose)
- fprintf(stderr,"hc_assign_density: reading density anomalies in [%%] from %s, scaling by %g\n",
- filename,hc->dens_scale[0]);
+ fprintf(stderr,"hc_assign_density: reading density anomalies in [%g%%] from %s\n",
+ 100*HC_DENSITY_SCALING,filename);
hc->inho = 0; /* counter for density layers */
/* get one density expansion */
hc_get_blank_expansions(&hc->dens_anom,1,0,
@@ -555,20 +582,19 @@
if(verbose)
fprintf(stderr,"hc_assign_density: using default SH format for density\n");
}
+
read_on = TRUE;
while(read_on){
if(use_short_format){
/* short format I/O */
- i = fscanf(in,"%lf",&dtmp);zlabel = (HC_PREC)dtmp;
+ i = fscanf(in,"%lf",dtmp);zlabel = (HC_PREC)dtmp[0];
i += fscanf(in,"%i",&lmax);
read_on = (i == 2)?(TRUE):(FALSE);
ivec = 0;shps = 1;type = HC_DEFAULT_INTERNAL_FORMAT;
ilayer++;
}else{
- read_on = sh_read_parameters_from_file(&type,&lmax,&shps,
- &ilayer, &nset,
- &zlabel,&ivec,in,
- FALSE,density_in_binary,
+ read_on = sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer, &nset,
+ &zlabel,&ivec,in,FALSE,density_in_binary,
verbose);
}
if(read_on){
@@ -584,7 +610,10 @@
lmax);
}
reported = TRUE;
+ if(verbose >= 2)
+ fprintf(stderr,"hc_assign_density: non_dim radius d\\rho/dinput %% factor PREM \\rho layer # depth[km]\n");
}
+
/*
do tests
*/
@@ -607,14 +636,18 @@
prem_get_rho(&rho0,hc->rden[hc->inho],hc->prem);
rho0 /= 1000.0;
/*
- general density (add additional depth dependence here)
+ general density
*/
- dens_scale[0] = hc->dens_scale[0] * local_dens_fac * rho0;
+ /* depth dependent factor? */
+ local_scale = hc_find_dens_scale(hc->rden[hc->inho],hc->dens_scale,
+ read_dens_scale_from_file,rdf,sdf,ndf);
+ /* total scaling factor */
+ dens_scale[0] = local_scale * ((HC_PREC)HC_DENSITY_SCALING ) * rho0;
if(verbose >= 2){
fprintf(stderr,"hc_assign_density: r: %11g anom scales: %11g x %11g x %11g = %11g\t%5i out of %i, z: %11g\n",
- hc->rden[hc->inho],hc->dens_scale[0],
- local_dens_fac,rho0,dens_scale[0],hc->inho+1,nset,zlabel);
+ hc->rden[hc->inho],local_scale,
+ HC_DENSITY_SCALING,rho0,dens_scale[0],hc->inho+1,nset,zlabel);
}
if(hc->inho){
/*
@@ -648,8 +681,7 @@
will assume input is in physical convention
*/
- sh_read_coefficients_from_file((hc->dens_anom+hc->inho),1,
- lmax,in,density_in_binary,
+ sh_read_coefficients_from_file((hc->dens_anom+hc->inho),1,lmax,in,density_in_binary,
dens_scale,verbose);
hc->inho++;
} /* end actualy read on */
@@ -731,10 +763,37 @@
sqrt(sh_total_power((hc->dens_anom+i))));
free(dbot);free(dtop);
} /* end layer structure part */
+
+ /* free density scaling factors, if used */
+ if(read_dens_scale_from_file && ndf){
+ free(rdf);
+ free(sdf);
+ }
+
hc->dens_init = TRUE;
}
/*
+find depth dependent scaling
+
+*/
+HC_PREC hc_find_dens_scale(HC_PREC r, HC_PREC s0,hc_boolean depth_dependent, HC_PREC *rd,HC_PREC *sd,int n)
+{
+ int i;
+ if(depth_dependent){
+ i=0;
+ while((i<n) && (rd[i] < r)){
+ i++;
+ }
+ i--;
+ return sd[i];
+ }else{
+ return s0;
+ }
+}
+
+/*
+
assign phase boundary jumps
input:
More information about the cig-commits
mailing list