[cig-commits] r13332 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Tue Nov 18 14:27:18 PST 2008
Author: becker
Date: 2008-11-18 14:27:17 -0800 (Tue, 18 Nov 2008)
New Revision: 13332
Added:
mc/1D/hc/trunk/plot_geoid
mc/1D/hc/trunk/sh_extract_layer.c
mc/1D/hc/trunk/visc.A
mc/1D/hc/trunk/visc.C
mc/1D/hc/trunk/visc.D
mc/1D/hc/trunk/visc.sh08
Modified:
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/ggrd_grdtrack_util.c
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/hc_output.c
mc/1D/hc/trunk/hc_polsol.c
mc/1D/hc/trunk/main.c
Log:
Added tentative computation of geoid at all depths.
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2008-11-18 17:40:27 UTC (rev 13331)
+++ mc/1D/hc/trunk/Makefile 2008-11-18 22:27:17 UTC (rev 13332)
@@ -132,7 +132,7 @@
all: dirs libs hc hc_extract_sh_layer \
- sh_syn sh_corr sh_ana sh_power rotvec2vel
+ sh_syn sh_corr sh_ana sh_power sh_extract_layer rotvec2vel
libs: dirs hc_lib $(HEAL_LIBS) $(RICK_LIB)
@@ -170,6 +170,12 @@
$(CC) $(LIB_FLAGS) $(ODIR)/sh_ana.o \
-o $(BDIR)/sh_ana -lhc -lrick $(HEAL_LIBS_LINKLINE) $(GGRD_LIBS_LINKLINE) -lm $(LDFLAGS)
+sh_extract_layer: $(LIBS) $(INCS) $(ODIR)/sh_extract_layer.o
+ $(CC) $(LIB_FLAGS) $(ODIR)/sh_extract_layer.o \
+ -o $(BDIR)/sh_extract_layer \
+ -lhc -lrick $(HEAL_LIBS_LINKLINE) $(GGRD_LIBS_LINKLINE) \
+ -lm $(LDFLAGS)
+
rotvec2vel: rotvec2vel.c
$(CC) $(CFLAGS) rotvec2vel.c -o $(BDIR)/rotvec2vel -lm $(LDFLAGS)
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-11-18 17:40:27 UTC (rev 13331)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-11-18 22:27:17 UTC (rev 13332)
@@ -438,7 +438,7 @@
FILE *din;
float dz1,dz2;
struct GRD_HEADER ogrd;
- int i,j,one_or_zero,nx,ny,mx,my,nn;
+ int i,one_or_zero,nx,ny,mx,my,nn;
char filename[BUFSIZ*2],*cdummy;
static int gmt_init = FALSE;
/*
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2008-11-18 17:40:27 UTC (rev 13331)
+++ mc/1D/hc/trunk/hc.h 2008-11-18 22:27:17 UTC (rev 13332)
@@ -61,6 +61,10 @@
#define HC_PI M_PI
+
+#define HC_SCALAR 0
+#define HC_VECTOR 1
+
/*
PREM
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2008-11-18 17:40:27 UTC (rev 13331)
+++ mc/1D/hc/trunk/hc_init.c 2008-11-18 22:27:17 UTC (rev 13332)
@@ -24,10 +24,12 @@
p->free_slip = TRUE; /* free slip? */
p->no_slip = FALSE; /* no slip boundary condition? */
p->platebc = FALSE; /* plate velocities? */
- p->compute_geoid = TRUE; /* compute the geoid? */
+ p->compute_geoid = 1; /* compute the geoid? 1: surface 2: all layers */
p->compute_geoid_correlations = FALSE; /* compute the geoid
correlation with
- refernece only */
+ refernece only
+ (only works for
+ surface) */
p->dens_anom_scale = HC_D_LOG_V_D_LOG_D ; /* default density anomaly scaling to
go from PREM percent traveltime
anomalies to density anomalies */
@@ -368,9 +370,10 @@
fprintf(stderr,"\t\tThe file (e.g. %s) is based on a DT expansion of cm/yr velocity fields.\n\n",HC_PVEL_FILE);
fprintf(stderr,"solution procedure and I/O options:\n");
- fprintf(stderr,"-ng\t\tdo not compute and print the geoid (%s)\n",
- hc_name_boolean(!p->compute_geoid));
- fprintf(stderr,"-rg\tname\tcompute correlation of computed geoid with that in file \"name\",\n");
+ fprintf(stderr,"-ng\t\tdo not compute and print the geoid (%i)\n",
+ p->compute_geoid);
+ fprintf(stderr,"-ag\t\tcompute geoid at all layer depths, as opposed to the surface only\n");
+ fprintf(stderr,"-rg\tname\tcompute correlation of surface geoid with that in file \"name\",\n");
fprintf(stderr,"\t\tthis will not print out the geoid file, but only correlations (%s)\n",
hc_name_boolean(p->compute_geoid_correlations));
fprintf(stderr,"-pptsol\t\tprint pol[6] and tor[2] solution vectors (%s)\n",
@@ -400,7 +403,9 @@
parameters */
hc_toggle_boolean(&p->print_pt_sol);
}else if(strcmp(argv[i],"-ng")==0){ /* do not compute geoid */
- hc_toggle_boolean(&p->compute_geoid);
+ p->compute_geoid = 0;
+ }else if(strcmp(argv[i],"-ag")==0){ /* compute geoid at all layers */
+ p->compute_geoid = 2;
}else if(strcmp(argv[i],"-rg")==0){ /* compute geoid correlations */
hc_toggle_boolean(&p->compute_geoid_correlations);
p->compute_geoid = TRUE;
Modified: mc/1D/hc/trunk/hc_output.c
===================================================================
--- mc/1D/hc/trunk/hc_output.c 2008-11-18 17:40:27 UTC (rev 13331)
+++ mc/1D/hc/trunk/hc_output.c 2008-11-18 22:27:17 UTC (rev 13332)
@@ -15,10 +15,8 @@
*/
-void hc_print_spectral_solution(struct hcs *hc,struct sh_lms *sol,
- FILE *out,int sol_mode,
- hc_boolean binary,
- hc_boolean verbose)
+void hc_print_spectral_solution(struct hcs *hc,struct sh_lms *sol,FILE *out,int sol_mode,
+ hc_boolean binary, hc_boolean verbose)
{
int i,os;
static int ntype = 3; /* three sets of solutions, r/pol/tor */
Modified: mc/1D/hc/trunk/hc_polsol.c
===================================================================
--- mc/1D/hc/trunk/hc_polsol.c 2008-11-18 17:40:27 UTC (rev 13331)
+++ mc/1D/hc/trunk/hc_polsol.c 2008-11-18 22:27:17 UTC (rev 13332)
@@ -190,7 +190,7 @@
// SUBROUTINE EVPPOT (L,RATIO,PPOT): OBTAINS PROPAGATOR FOR
// NON-EQUILIBRIUM POTENTIAL AND DERIVATIVE (RATIO IS R(I)/
// R(I+1), FOR PROPAGATION FROM R(I) TO R(I+1) AT L),
- int i,i2,i3,i6,j,l,m,nih,nxtv,ivis,os,pos1,pos2,
+ int i,i2,i3,i6,j,l,m,nih,nxtv,ivis,os,pos1,pos2,gi,g1,g2,gic,
prop_s1,prop_s2,nvisp1,nzero,n6,ninho,nl=0,ip1;
int newprp,newpot,jpb,inho2,ibv,indx[3],a_or_b,ilayer,lmax,
nprops_max;
@@ -809,10 +809,10 @@
if(verbose)
fprintf(stderr,"hc_polsol: assigned nl: %i nprop: %i nrad: %i layers\n",
nl,hc->nprops,nrad);
- if(nl != nrad+2){
+ if(nl != hc->nradp2){
HC_ERROR("hc_polsol","nl not equal to nrad+2 at end of solution loop");
}
- if (compute_geoid){
+ if(compute_geoid){
//
// Calculating geoid coefficients. The factor gf comes from
// * u(5) is in units of r0 * pi/180 / Ma (about 111 km/Ma)
@@ -820,34 +820,54 @@
// * geoid is in units of meters
//
if(verbose > 1)
- fprintf(stderr,"hc_polsol: evaluating geoid\n");
+ fprintf(stderr,"hc_polsol: evaluating geoid%s\n",
+ (compute_geoid == 1)?(" at surface"):(", all layers"));
+ /*
+ select geoid solution
+ */
n6 = 4;
//n6 = -iformat-1;
/* first coefficients are zero */
clm[0] = clm[1] = 0.0;
- sh_write_coeff(geoid,0,0,0,FALSE,clm); /* 0,0 */
- sh_write_coeff(geoid,1,0,0,FALSE,clm); /* 1,0 */
- sh_write_coeff(geoid,1,1,2,FALSE,clm); /* 1,1 */
+ switch(compute_geoid){
+ case 1:
+ g1 = hc->nrad+1;g2=hc->nradp2; /* only surface */
+ break;
+ case 2:
+ g1 = 0;g2=hc->nradp2; /* all layers */
+ break;
+ default:
+ fprintf(stderr,"hc_polsol: error, geoid = %i undefined\n",compute_geoid);
+ exit(-1);
+ }
+ for(gic=0,gi=g1;gi < g2;gi++,gic++){ /* depth loop */
+ /*
+ first coefficients
+ */
+ sh_write_coeff((geoid+gic),0,0,0,FALSE,clm); /* 0,0 */
+ sh_write_coeff((geoid+gic),1,0,0,FALSE,clm); /* 1,0 */
+ sh_write_coeff((geoid+gic),1,1,2,FALSE,clm); /* 1,1 */
- os = (nrad+1) * 6 + n6; /* select component */
- for(l=2;l <= pol_sol[0].lmax;l++){
- for(m=0;m <= l;m++){
- if (m != 0){
- sh_get_coeff((pol_sol+os),l,m,2,FALSE,clm); /* internal convention */
- clm[0] *= hc->psp.geoid_factor;
- clm[1] *= hc->psp.geoid_factor;
- sh_write_coeff(geoid,l,m,2,FALSE,clm);
- }else{ /* m == 0 */
- sh_get_coeff((pol_sol+os),l,m,0,FALSE,clm);
- clm[0] *= hc->psp.geoid_factor;
- sh_write_coeff(geoid,l,m,0,FALSE,clm);
+ os = gi * 6 + n6; /* select component */
+ for(l=2;l <= pol_sol[0].lmax;l++){
+ for(m=0;m <= l;m++){
+ if (m != 0){
+ sh_get_coeff((pol_sol+os),l,m,2,FALSE,clm); /* internal convention */
+ clm[0] *= hc->psp.geoid_factor;
+ clm[1] *= hc->psp.geoid_factor;
+ sh_write_coeff((geoid+gic),l,m,2,FALSE,clm);
+ }else{ /* m == 0 */
+ sh_get_coeff((pol_sol+os),l,m,0,FALSE,clm);
+ clm[0] *= hc->psp.geoid_factor;
+ sh_write_coeff((geoid+gic),l,m,0,FALSE,clm);
+ }
}
}
}
if(verbose > 1)
fprintf(stderr,"hc_polsol: assigned geoid\n");
- }
+ } /* end geoid */
/*
free the local arrays
*/
Modified: mc/1D/hc/trunk/main.c
===================================================================
--- mc/1D/hc/trunk/main.c 2008-11-18 17:40:27 UTC (rev 13331)
+++ mc/1D/hc/trunk/main.c 2008-11-18 22:27:17 UTC (rev 13332)
@@ -52,7 +52,7 @@
struct hcs *model; /* main structure, make sure to initialize with
zeroes */
struct sh_lms *sol_spectral=NULL, *geoid = NULL; /* solution expansions */
- int nsol,lmax,solved;
+ int nsol,lmax,solved,i;
FILE *out;
struct hc_parameters p[1]; /* parameters */
char filename[HC_CHAR_LENGTH],file_prefix[10];
@@ -60,6 +60,8 @@
e.g. velocities */
HC_PREC corr[2]; /* correlations */
HC_PREC vl[4][3],v[4],dv; /* for viscosity scans */
+ static hc_boolean geoid_binary = FALSE; /* type of geoid output */
+ static HC_CPREC unitya[1] = {1.0};
/*
@@ -110,7 +112,17 @@
*/
hc_init_main(model,SH_RICK,p);
- nsol = (model->nrad+2) * 3; /* number of solution (r,pol,tor)*(nlayer+2) */
+ nsol = (model->nradp2) * 3; /*
+ number of solutions (r,pol,tor) * (nlayer+2)
+
+ total number of layers is nlayer +2,
+
+ because CMB and surface are added
+ to intermediate layers which are
+ determined by the spacing of the
+ density model
+
+ */
if(p->free_slip) /* maximum degree is determined by the
density expansion */
lmax = model->dens_anom[0].lmax;
@@ -131,12 +143,16 @@
/*
make room for the spectral solution on irregular grid
*/
- sh_allocate_and_init(&sol_spectral,nsol,lmax,model->sh_type,1,
+ sh_allocate_and_init(&sol_spectral,nsol,lmax,model->sh_type,HC_VECTOR,
p->verbose,FALSE);
- if(p->compute_geoid)
- /* make room for geoid solution */
+ if(p->compute_geoid == 1)
+ /* make room for geoid solution at surface */
sh_allocate_and_init(&geoid,1,model->dens_anom[0].lmax,
- model->sh_type,0,p->verbose,FALSE);
+ model->sh_type,HC_SCALAR,p->verbose,FALSE);
+ else if(p->compute_geoid == 2) /* all layers */
+ sh_allocate_and_init(&geoid,model->nradp2,model->dens_anom[0].lmax,
+ model->sh_type,HC_SCALAR,p->verbose,FALSE);
+
switch(p->solver_mode){
case HC_SOLVER_MODE_DEFAULT:
/*
@@ -177,6 +193,11 @@
fclose(out);
if(p->compute_geoid){
if(p->compute_geoid_correlations){
+ if(p->compute_geoid == 2){ /* check if all geoids were computed */
+ fprintf(stderr,"%s: can only compute correlation for surface geoid, geoid = %i\n",
+ argv[0],p->compute_geoid);
+ exit(-1);
+ }
if(p->verbose)
fprintf(stderr,"%s: correlation for geoid with %s\n",argv[0],
p->ref_geoid_file);
@@ -188,9 +209,18 @@
*/
sprintf(filename,"%s",HC_GEOID_FILE);
if(p->verbose)
- fprintf(stderr,"%s: writing geoid to %s\n",argv[0],filename);
- out = ggrd_open(filename,"w","main");
- hc_print_sh_scalar_field(geoid,out,FALSE,FALSE,p->verbose);
+ fprintf(stderr,"%s: writing geoid to %s, %s\n",
+ argv[0],filename,(p->compute_geoid == 1)?("at surface"):("all layers"));
+ out = ggrd_open(filename,"w","main");
+ if(p->compute_geoid == 1) /* surface layer */
+ hc_print_sh_scalar_field(geoid,out,FALSE,geoid_binary,p->verbose);
+ else{ /* all layers */
+ for(i=0;i < model->nradp2;i++){
+ sh_print_parameters_to_file((geoid+i),1,i,model->nradp2,
+ HC_Z_DEPTH(model->r[i]),out,FALSE,geoid_binary,p->verbose);
+ sh_print_coefficients_to_file((geoid+i),1,out,unitya,geoid_binary,p->verbose);
+ }
+ }
fclose(out);
}
}
@@ -276,9 +306,10 @@
*/
sh_free_expansion(sol_spectral,nsol);
- if(p->compute_geoid)
+ if(p->compute_geoid == 1) /* only one layer */
sh_free_expansion(geoid,1);
-
+ else if(p->compute_geoid == 2) /* all layers */
+ sh_free_expansion(geoid,model->nradp2);
free(sol_spatial);
if(p->verbose)
fprintf(stderr,"%s: done\n",argv[0]);
Added: mc/1D/hc/trunk/plot_geoid
===================================================================
--- mc/1D/hc/trunk/plot_geoid (rev 0)
+++ mc/1D/hc/trunk/plot_geoid 2008-11-18 22:27:17 UTC (rev 13332)
@@ -0,0 +1,33 @@
+#!/bin/bash
+#
+# simple geoid plotting tester
+#
+
+makecpt -T-200/200/10 -Chaxby > geoid.cpt
+inc=-I1 # default res of sh_sysn
+reg=-Rg
+proj=-JQ180/7 # projection
+ann=-Ba60f30WeSn
+#
+# surface geoid
+hc
+# expand on regular grid
+cat geoid.ab | sh_syn 0 0 0 | xyz2grd $reg $inc -Gtop1.grd
+
+#
+# geoid all layers
+#
+hc -ag
+
+# bottom
+sh_extract_layer geoid.ab 1 | sh_syn 0 0 0 | xyz2grd $reg $inc -Gbot.grd
+# top
+sh_extract_layer geoid.ab -1 | sh_syn 0 0 0 | xyz2grd $reg $inc -Gtop2.grd
+
+
+for g in top1 top2 bot ;do
+ grdimage $reg -Cgeoid.cpt $g.grd $proj $ann -P -K > $g.ps
+ pscoast -Dc -A50000 $reg $proj -O -K -W5 >> $g.ps
+ psscale -D3.5/-.5/3/.2h -Cgeoid.cpt -B50/:"[m]": -O >> $g.ps
+ echo $0: plotted to $g.ps
+done
Added: mc/1D/hc/trunk/sh_extract_layer.c
===================================================================
--- mc/1D/hc/trunk/sh_extract_layer.c (rev 0)
+++ mc/1D/hc/trunk/sh_extract_layer.c 2008-11-18 22:27:17 UTC (rev 13332)
@@ -0,0 +1,82 @@
+#include "hc.h"
+/*
+
+
+extract part of a spherical harmonics file that has several layers,
+such as the hc geoid.ab output for the case of all layers (hc -ag)
+
+*/
+
+int main(int argc, char **argv)
+{
+ int ilayer,i,shps,nset,ivec,lmax,type;
+ FILE *in;
+ struct sh_lms *exp=NULL;
+ HC_PREC unitya[3] = {1.0,1.0,1.0},*zdepth;
+ hc_boolean binary = FALSE, verbose = TRUE, short_format = FALSE;
+ /*
+ deal with parameters
+ */
+ hc_vecalloc(&zdepth,1,"");
+ ilayer = 0;
+ switch(argc){
+ case 3:
+ sscanf(argv[2],"%i",&ilayer);
+ break;
+ default:
+ fprintf(stderr,"%s: usage\n\n%s value.ab layer\n\n",argv[0],argv[0]);
+ fprintf(stderr,"extracts one SH layer (e.g. for use in sh_syn) from a (long format, hc) spherical harmonic file value.ab\n");
+ fprintf(stderr,"layer: 1...nset\n");
+ fprintf(stderr,"\tif ilayer= 1..nset, will print one layer\n");
+ fprintf(stderr,"\t -1, will select nset\n");
+ exit(-1);
+ break;
+ }
+ /*
+ read in spherical harmonics
+ */
+ in = ggrd_open(argv[1],"r","sh_extract_layer");
+ /* start loop */
+ i = 0;
+ sh_read_parameters_from_file(&type,&lmax, &shps,&i,&nset,(zdepth+i),
+ &ivec,in,short_format,binary,verbose);
+ if(verbose)
+ fprintf(stderr,"sh_extract_layer: detected %i layers, vec: %i, type %i SH file, shps: %i, layer %i at depth %g\n",
+ nset,ivec,type,shps,i+1,zdepth[i]);
+ sh_allocate_and_init(&exp,nset*shps,lmax,type,ivec,verbose,FALSE);
+ hc_vecrealloc(&zdepth,nset,"sh_extract_layer: zdepth mem error");
+ /* which layer to select */
+ if(ilayer == -1)
+ ilayer = nset - 1;
+ else
+ ilayer--;
+ /* check bounds */
+ if(ilayer < 0){
+ fprintf(stderr,"sh_extract_layer: ilayer should be given between 1 and %i for file %s\n",
+ nset,argv[1]);
+ exit(-1);
+ }
+ for(;i<nset;i++){
+ if(i != 0)
+ sh_read_parameters_from_file(&type,&lmax,&shps,&i,&nset,(zdepth+i),
+ &ivec,in,short_format,binary,verbose);
+ sh_read_coefficients_from_file((exp+i),shps,-1,in,binary,unitya,
+ verbose);
+ if(i == ilayer){ /* output */
+ if(verbose)
+ fprintf(stderr,"%s: printing SH from %s at layer %i out of %i to stdout (depth: %g)\n",
+ argv[0],argv[1],i+1,nset,zdepth[i]);
+ /*
+ output will remove the layer number information
+ */
+ sh_print_parameters_to_file((exp+i),shps,ilayer,1,zdepth[i],
+ stdout,short_format,FALSE,verbose);
+ sh_print_coefficients_to_file((exp+i),shps,stdout,unitya,FALSE,verbose);
+ }
+ }
+ fclose(in);
+
+ sh_free_expansion(exp,nset);
+
+ return 0;
+}
Added: mc/1D/hc/trunk/visc.A
===================================================================
--- mc/1D/hc/trunk/visc.A (rev 0)
+++ mc/1D/hc/trunk/visc.A 2008-11-18 22:27:17 UTC (rev 13332)
@@ -0,0 +1 @@
+0.546225 1e21
Added: mc/1D/hc/trunk/visc.C
===================================================================
--- mc/1D/hc/trunk/visc.C (rev 0)
+++ mc/1D/hc/trunk/visc.C 2008-11-18 22:27:17 UTC (rev 13332)
@@ -0,0 +1,3 @@
+0.546225 5e22
+0.8964 1e21
+0.984 5e22
Added: mc/1D/hc/trunk/visc.D
===================================================================
--- mc/1D/hc/trunk/visc.D (rev 0)
+++ mc/1D/hc/trunk/visc.D 2008-11-18 22:27:17 UTC (rev 13332)
@@ -0,0 +1,4 @@
+0.546225 5e22
+0.8964 1e21
+0.9356 1e20
+0.9843 5e22
Added: mc/1D/hc/trunk/visc.sh08
===================================================================
--- mc/1D/hc/trunk/visc.sh08 (rev 0)
+++ mc/1D/hc/trunk/visc.sh08 2008-11-18 22:27:17 UTC (rev 13332)
@@ -0,0 +1,22 @@
+0.5448 0.186E+21
+0.5668 0.354E+22
+0.5888 0.332E+23
+0.6107 0.885E+23
+0.6327 0.112E+24
+0.6547 0.111E+24
+0.6767 0.101E+24
+0.6986 0.875E+23
+0.7206 0.731E+23
+0.7426 0.586E+23
+0.7646 0.451E+23
+0.7865 0.332E+23
+0.8085 0.235E+23
+0.8305 0.159E+23
+0.8525 0.103E+23
+0.8744 0.624E+22
+0.8964 0.651E+21
+0.9184 0.431E+21
+0.9356 0.432E+21
+0.9466 0.290E+21
+0.9655 0.241E+21
+0.9843 0.264E+23
More information about the CIG-COMMITS
mailing list