Added tentative computation of geoid at all depths.

Makefile

@@ -132,7 +132,7 @@ LIBS = $(HC_LIBS) $(GGRD_LIBS) $(HEAL_LIBS) $(RICK_LIB)
 
 
 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 @@ sh_ana: $(LIBS) $(INCS) $(ODIR)/sh_ana.o
 	$(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)
 

ggrd_grdtrack_util.c

@@ -438,7 +438,7 @@ int ggrd_grdtrack_init(double *west, double *east,double *south, double *north,
   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;
   /* 

hc.h

@@ -61,6 +61,10 @@ struct hc_dcmplx{			/* double precision */
 
 #define HC_PI M_PI
 
+
+#define HC_SCALAR 0
+#define HC_VECTOR 1
+
 /* 
 
 PREM

hc_init.c

@@ -24,10 +24,12 @@ void hc_init_parameters(struct hc_parameters *p)
   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 @@ void hc_handle_command_line(int argc, char **argv,
       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 @@ void hc_handle_command_line(int argc, char **argv,
 						   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;

hc_output.c

@@ -15,10 +15,8 @@ print the spherical harmonics version of a solution set
 
 
 */
-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 */

hc_polsol.c

@@ -190,7 +190,7 @@ void hc_polsol(struct hcs *hc, 	/*
   //       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,45 +809,65 @@ void hc_polsol(struct hcs *hc, 	/*
   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) 
     //    * normalizing density is presumably 1 g/cm**3 = 1000 kg / m**3
     //    * 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 */
-
-    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);
+    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 = 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 
   */

main.c

@@ -52,14 +52,16 @@ int main(int argc, char **argv)
   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];
   float *sol_spatial = NULL;	/* spatial solution,
 				   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 @@ int main(int argc, char **argv)
 
   */
   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 @@ int main(int argc, char **argv)
   /* 
      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 @@ int main(int argc, char **argv)
     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 @@ int main(int argc, char **argv)
 	*/
 	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 @@ int main(int argc, char **argv)
 
   */
   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]);

plot_geoid

@@ -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