[cig-commits] r11265 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Tue Feb 26 21:54:15 PST 2008
Author: becker
Date: 2008-02-26 21:54:14 -0800 (Tue, 26 Feb 2008)
New Revision: 11265
Modified:
mc/1D/hc/trunk/ggrd.h
mc/1D/hc/trunk/ggrd_grdtrack_util.c
mc/1D/hc/trunk/ggrd_grdtrack_util.h
mc/1D/hc/trunk/ggrd_readgrds.c
mc/1D/hc/trunk/ggrd_struc.h
mc/1D/hc/trunk/ggrd_test.c
mc/1D/hc/trunk/ggrd_velinterpol.c
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_extract_sh_layer.c
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/hc_misc.c
mc/1D/hc/trunk/hc_output.c
mc/1D/hc/trunk/main.c
mc/1D/hc/trunk/sh_test.c
Log:
More internal changes to interact better with CitcomS CIG version.
Modified: mc/1D/hc/trunk/ggrd.h
===================================================================
--- mc/1D/hc/trunk/ggrd.h 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/ggrd.h 2008-02-27 05:54:14 UTC (rev 11265)
@@ -26,7 +26,7 @@
/* errors */
#define GGRD_PE(x) {fprintf(stderr,"ggrd: %s\n",x);}
-
+#define GGRD_MEMERROR(x) {fprintf(stderr,"%s: memory allocation error, exiting\n",x);exit(-1);}
/* radius of CMB */
#define GGRD_RCMB_ND 0.546225
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-02-27 05:54:14 UTC (rev 11265)
@@ -24,10 +24,15 @@
void ggrd_init_master(struct ggrd_master *ggrd)
{
ggrd->mat_control = ggrd->mat_control_init = 0;
- ggrd->vel_control = ggrd->vel_control_init = 0;
+ ggrd->vtop_control = ggrd->vtop_control_init = 0;
+ ggrd->age_control = ggrd->age_control_init = 0;
+ ggrd->nage = 0;ggrd->age_bandlim = 200.;
+ ggrd->sf_init = 0;
ggrd->time_hist.init = 0;
ggrd->temp_init.init = 0;
ggrd->time_hist.vstage_transition = 0.1; /* in Ma, transition */
+ /* 3-D velocity settings */
+ ggrd_init_vstruc(ggrd);
}
@@ -1088,78 +1093,78 @@
*/
int interpolate_seafloor_ages(GGRD_CPREC xt, GGRD_CPREC xp,
- GGRD_CPREC age,struct ggrd_vel *v,
+ GGRD_CPREC age,struct ggrd_master *ggrd,
GGRD_CPREC *seafloor_age)
{
int left, right,i;
GGRD_CPREC f1,f2;
double a1,a2;
- if(!v->sf_init){
+ if(!ggrd->sf_init){
/*
init the constants
*/
- v->sf_ntlim = v->nage - 1;
- if(!v->nage){
+ ggrd->sf_ntlim = ggrd->nage - 1;
+ if(!ggrd->nage){
fprintf(stderr,"interpolate_seafloor_ages: not initialized?!\n");
return -2;
}
- for(i=1;i < v->nage;i++)
- if(v->age_time[i] < v->age_time[i-1]){
+ for(i=1;i < ggrd->nage;i++)
+ if(ggrd->age_time[i] < ggrd->age_time[i-1]){
fprintf(stderr,"interpolate_seafloor_ages: error: times need to be sorted monotnically increasing\n");
return -3;
}
- v->sf_old_age = v->age_time[0] - 1000; /* so that we get interpolation
+ ggrd->sf_old_age = ggrd->age_time[0] - 1000; /* so that we get interpolation
factors */
- v->sf_init = TRUE;
+ ggrd->sf_init = TRUE;
}
/* check range */
- if((age < v->age_time[0]) || (age > v->age_time[v->sf_ntlim])){
+ if((age < ggrd->age_time[0]) || (age > ggrd->age_time[ggrd->sf_ntlim])){
fprintf(stderr,"interpolate_seafloor_ages: age: %g out of bounds [%g;%g]\n",
- age,v->age_time[0],v->age_time[v->sf_ntlim]);
+ age,ggrd->age_time[0],ggrd->age_time[ggrd->sf_ntlim]);
return -3;
}
- if(fabs(age- v->sf_old_age) > 1e-8){
+ if(fabs(age- ggrd->sf_old_age) > 1e-8){
/*
time interpolation
*/
right = 0;
- while((right < v->sf_ntlim) && (v->age_time[right] < age))
+ while((right < ggrd->sf_ntlim) && (ggrd->age_time[right] < age))
right++;
if(right == 0)
right++;
left = right - 1;
- f2 = (age - v->age_time[left])/
- (v->age_time[right]-v->age_time[left]);
+ f2 = (age - ggrd->age_time[left])/
+ (ggrd->age_time[right]-ggrd->age_time[left]);
f1 = 1.0-f2;
- //fprintf(stderr,"sai: %g %g %g\t%g %g\n",v->age_time[left],age,v->age_time[right],f1,f2);
- v->sf_old_age = age;
- v->sf_old_left = left;v->sf_old_right = right;
- v->sf_old_f1 = f1;v->sf_old_f2 = f2;
+ //fprintf(stderr,"sai: %g %g %g\t%g %g\n",ggrd->age_time[left],age,ggrd->age_time[right],f1,f2);
+ ggrd->sf_old_age = age;
+ ggrd->sf_old_left = left;ggrd->sf_old_right = right;
+ ggrd->sf_old_f1 = f1;ggrd->sf_old_f2 = f2;
}else{ /* reuse time interpolation */
- left = v->sf_old_left;right = v->sf_old_right;
- f1 = v->sf_old_f1;f2 = v->sf_old_f2;
+ left = ggrd->sf_old_left;right = ggrd->sf_old_right;
+ f1 = ggrd->sf_old_f1;f2 = ggrd->sf_old_f2;
}
if(!ggrd_grdtrack_interpolate_tp((double)xt,(double)xp,
- (v->ages+left),&a1,FALSE)){
+ (ggrd->ages+left),&a1,FALSE)){
fprintf(stderr,"interpolate_seafloor_ages: interpolation error left\n");
return -6;
}
- if(a1 > v->ages[left].fmaxlim[0]) /* limit to bandlim */
- a1 = v->ages[left].bandlim;
+ if(a1 > ggrd->ages[left].fmaxlim[0]) /* limit to bandlim */
+ a1 = ggrd->ages[left].bandlim;
if(!ggrd_grdtrack_interpolate_tp((double)xt,(double)xp,
- (v->ages+right),&a2,FALSE)){
+ (ggrd->ages+right),&a2,FALSE)){
fprintf(stderr,"interpolate_seafloor_ages: interpolation error right\n");
return -6;
}
- if(a2 > v->ages[right].fmaxlim[0])
- a2 = v->ages[right].bandlim;
+ if(a2 > ggrd->ages[right].fmaxlim[0])
+ a2 = ggrd->ages[right].bandlim;
*seafloor_age = (GGRD_CPREC)(f1 * a1 + f2 * a2);
if(*seafloor_age < 0)
*seafloor_age = 0.0; /* and >= 0 */
@@ -1167,6 +1172,20 @@
}
+/*
+ open a file safely and give an error message if there was
+ a problem
+*/
+FILE *ggrd_open(char *name, char *mode, char *program)
+{
+ FILE *in;
+ if((in=fopen(name,mode)) == NULL){
+ fprintf(stderr,"%s: error: can not open file %s for mode %s access\n",
+ program,name,mode);
+ exit(-1);
+ }
+ return in;
+}
/* general floating point vector allocation */
void ggrd_vecalloc(double **x,int n,char *message)
@@ -1188,7 +1207,154 @@
}
}
+/*
+ calculate mean and standard deviation of x_i
+ if hypoth is set, calc mean and stddev of sqrt(x_i^2 + y_i^2)
+ if weighted is set, uses weights[n] for weighting the mean and so on
+
+ this way of calculating the stddev is inaccurate but fast
+
+*/
+void ggrd_calc_mean_and_stddev(GGRD_CPREC *x, GGRD_CPREC *y,int n,GGRD_CPREC *mean,
+ GGRD_CPREC *stddev,GGRD_CPREC *rms,
+ ggrd_boolean hypoth, ggrd_boolean weighted,GGRD_CPREC *weight)
+{
+ GGRD_CPREC sum1=0.0,sum2=0.0,tmp,ws;
+ int i;
+ if(n <= 1){
+ fprintf(stderr,"ggrd_calc_mean_and_stddev: error: n: %i\n",n);
+ exit(-1);
+ }
+ ws=0.0;
+ if(hypoth){// sqrt(x^2+y+2)
+ for(i=0;i<n;i++){
+ if(weighted){
+ tmp = hypot(x[i],y[i]) * weight[i];
+ ws += weight[i];
+ }else{
+ tmp = hypot(x[i],y[i]);
+ ws += 1.0;
+ }
+ sum1 += tmp;sum2 += tmp * tmp;
+ }
+ }else{
+ for(i=0;i<n;i++){
+ if(weighted){
+ tmp = x[i] * weight[i];
+ ws += weight[i];
+ }else{
+ tmp = x[i];
+ ws += 1.0;
+ }
+ sum1 += tmp;sum2 += tmp*tmp;
+ }
+ }
+ // standard deviation
+ tmp = (ws * sum2 - sum1 * sum1) / (ws*(ws-1.0));
+ if(tmp > 0)
+ *stddev = sqrt(tmp);
+ else
+ *stddev = 0.0;
+ *rms = sqrt(sum2 / ws);// RMS
+ *mean = sum1 / ws; // mean
+}
+/*
+
+ sort array by order, returns index
+
+ GIVE INPUT INDEX SHIFTED
+ (ie. x-1)
+ output will be 0...n-1
+
+ based on numerical recipes
+
+
+*/
+
+
+#define GGRD_INDEXX_SWAP(a,b) itemp=(a);(a)=(b);(b)=itemp;
+#define GGRD_INDEXX_M 7
+#define GGRD_INDEXX_NSTACK 5000
+
+void ggrd_indexx(int n,GGRD_CPREC *arr, int *indx)
+{
+ int i,indxt,ir=n,itemp,j,k,l=1;
+ int jstack=0,*istack;
+ GGRD_CPREC a;
+ //
+ istack=(int *)malloc(sizeof(int)*(GGRD_INDEXX_NSTACK+1));
+ if(!istack)
+ GGRD_MEMERROR("indexx");
+ //
+ for (j=1;j<=n;j++)
+ indx[j]=j;
+ for (;;) {
+ if (ir-l < GGRD_INDEXX_M) {
+ for (j=l+1;j<=ir;j++) {
+ indxt=indx[j];
+ a=arr[indxt];
+ for (i=j-1;i>=1;i--) {
+ if (arr[indx[i]] <= a) break;
+ indx[i+1]=indx[i];
+ }
+ indx[i+1]=indxt;
+ }
+ if (jstack == 0) break;
+ ir=istack[jstack--];
+ l=istack[jstack--];
+ } else {
+ k=(l+ir) >> 1;
+ GGRD_INDEXX_SWAP(indx[k],indx[l+1]);
+ if (arr[indx[l+1]] > arr[indx[ir]]) {
+ GGRD_INDEXX_SWAP(indx[l+1],indx[ir]);
+ }
+ if (arr[indx[l]] > arr[indx[ir]]) {
+ GGRD_INDEXX_SWAP(indx[l],indx[ir]);
+ }
+ if (arr[indx[l+1]] > arr[indx[l]]) {
+ GGRD_INDEXX_SWAP(indx[l+1],indx[l]);
+ }
+ i=l+1;
+ j=ir;
+ indxt=indx[l];
+ a=arr[indxt];
+ for (;;) {
+ do i++; while (arr[indx[i]] < a);
+ do j--; while (arr[indx[j]] > a);
+ if (j < i) break;
+ GGRD_INDEXX_SWAP(indx[i],indx[j]);
+ }
+ indx[l]=indx[j];
+ indx[j]=indxt;
+ jstack += 2;
+ if (jstack > GGRD_INDEXX_NSTACK) {
+ fprintf(stderr,
+ "indexx: GGRD_INDEXX_NSTACK too small");
+ exit(-1);
+ }
+ if (ir-i+1 >= j-l) {
+ istack[jstack]=ir;
+ istack[jstack-1]=i;
+ ir=j-1;
+ } else {
+ istack[jstack]=j-1;
+ istack[jstack-1]=l;
+ l=i;
+ }
+ }
+ }
+ free(istack);
+ // go back to normal numbering
+ for(i=1;i<=n;i++)
+ indx[i]--;
+}
+#undef GGRD_INDEXX_M
+#undef GGRD_INDEXX_NSTACK
+#undef GGRD_INDEXX_SWAP
+
+
+
/*
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.h
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.h 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.h 2008-02-27 05:54:14 UTC (rev 11265)
@@ -18,6 +18,12 @@
#include <string.h>
#include <math.h>
+#ifdef USE_GMT4
+#define GGRD_GMT_GLOBAL_STRING "-fg"
+#else
+#define GGRD_GMT_GLOBAL_STRING "-Lg"
+#endif
+
/*
wrappers
@@ -55,7 +61,7 @@
void ggrd_init_master(struct ggrd_master *);
int ggrd_init_thist_from_file(struct ggrd_t *,char *,ggrd_boolean ,ggrd_boolean);
-int ggrd_read_vel_grids(struct ggrd_vel *, double, unsigned short, unsigned short, char *);
+int ggrd_read_vel_grids(struct ggrd_master *, double, unsigned short, unsigned short, char *);
#ifdef USE_GMT4
/* GMT4.1.2 */
@@ -108,6 +114,7 @@
ggrd_boolean *); /* */
float ggrd_gt_rms(float *,int );
float ggrd_gt_mean(float *,int );
+FILE *ggrd_open(char *, char *, char *);
void ggrd_print_layer_avg(float *,float *,int , int ,FILE *);
@@ -121,3 +128,8 @@
void ggrd_vecalloc(double **,int,char *);
void ggrd_vecrealloc(double **,int,char *);
+void ggrd_indexx(int ,GGRD_CPREC *, int *);
+void ggrd_calc_mean_and_stddev(GGRD_CPREC *, GGRD_CPREC *,int ,GGRD_CPREC *,
+ GGRD_CPREC *,GGRD_CPREC *,
+ ggrd_boolean , ggrd_boolean,GGRD_CPREC *);
+
Modified: mc/1D/hc/trunk/ggrd_readgrds.c
===================================================================
--- mc/1D/hc/trunk/ggrd_readgrds.c 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/ggrd_readgrds.c 2008-02-27 05:54:14 UTC (rev 11265)
@@ -44,34 +44,25 @@
int finite(double ); /* why? */
/* init a v structure */
-void ggrd_init_vstruc(struct ggrd_vel *v)
+void ggrd_init_vstruc(struct ggrd_master *ggrd)
{
- v->read_gmt = TRUE; /* read GMT by default */
- v->amode = GGRD_NORMAL;
- v->init = FALSE;
- v->history = FALSE;
- v->use_age = FALSE;
- v->age_bandlim = 200;
- v->nage = 0;
- v->vr = v->vt = v->vp = NULL;
- v->thist.init = FALSE;
- v->velscale = 1.0;
- v->rcmb = GGRD_RCMB_ND;
+ /* directly velocity related */
+ ggrd->v.read_gmt = TRUE; /* read GMT by default */
+ ggrd->v.init = FALSE;
+ ggrd->v.history = FALSE;
+
+ ggrd->v.vr = ggrd->v.vt = ggrd->v.vp = NULL;
+ ggrd->v.velscale = 1.0;
+ ggrd->v.rcmb = GGRD_RCMB_ND;
/* rlevels */
- v->rlevels = NULL;
- v->rl_warned = FALSE;
- /* seafloor */
- v->sf_init = FALSE;
- /* thist structure */
- v->thist.init = FALSE;
- v->thist.called = FALSE;
-
+ ggrd->v.rlevels = NULL;
+ ggrd->v.rl_warned = FALSE;
/* interpolation structure */
- v->vd.init = FALSE;
- v->vd.reduce_r_stencil = FALSE;
- v->vd.z_warned = FALSE;
- v->vd.w_warned = FALSE;
-
+ ggrd->v.vd.init = FALSE;
+ ggrd->v.vd.reduce_r_stencil = FALSE;
+ ggrd->v.vd.z_warned = FALSE;
+ ggrd->v.vd.w_warned = FALSE;
+
}
/*
@@ -87,10 +78,9 @@
*/
-int ggrd_read_vel_grids(struct ggrd_vel *v, /* velocity structure,
- should be initialized first
-
- */
+int ggrd_read_vel_grids(struct ggrd_master *ggrd, /* ggrd master structure
+ should be initialized first
+ */
GGRD_CPREC scale, /* divide all velocities by this
factor */
hc_boolean verbose, /* verbosity level */
@@ -126,21 +116,21 @@
maxtheta = HC_FLT_MIN;
weights=NULL;
- v->velscale = scale;
- if(fabs(v->velscale) < HC_EPS_PREC){
+ ggrd->v.velscale = scale;
+ if(fabs(ggrd->v.velscale) < HC_EPS_PREC){
fprintf(stderr,"ggrd_read_vel_grids: error: velocity scale is zero\n");
return(-1);
}
- if(!v->init){
+ if(!ggrd->v.init){
//
// read time intervals for velocities from file
sprintf(tfilename,"%s%s",prefix,GGRD_THFILE);
/*
- if v->history is set, will look for different time intervals
+ if ggrd->v.history is set, will look for different time intervals
*/
- ggrd_init_thist_from_file(&v->thist,tfilename,v->history,verbose);
- if(v->use_age){
+ ggrd_init_thist_from_file(&ggrd->time_hist,tfilename,ggrd->v.history,verbose);
+ if(ggrd->age_control){
/*
initialize seafloor ages specified at the beginning of each
@@ -149,43 +139,43 @@
*/
- if(!v->history){
+ if(!ggrd->v.history){
fprintf(stderr,"ggrd_read_vel_grids: error: for ages, need history input\n");
return(-1);
}
if(verbose)
fprintf(stderr,"ggrd_read_vel_grids: expecting %i (nt) + 1 age grids\n",
- v->thist.nvtimes);
- v->nage = v->thist.nvtimes + 1;
+ ggrd->time_hist.nvtimes);
+ ggrd->nage = ggrd->time_hist.nvtimes + 1;
/* important to use calloc so that some flags are set to zero */
- v->ages = (struct ggrd_gt *)calloc(v->nage,
+ ggrd->ages = (struct ggrd_gt *)calloc(ggrd->nage,
sizeof(struct ggrd_gt));
- v->age_time = (GGRD_CPREC *)malloc(v->nage*sizeof(GGRD_CPREC));
- if(!v->ages || ! v->age_time){
+ ggrd->age_time = (GGRD_CPREC *)malloc(ggrd->nage*sizeof(GGRD_CPREC));
+ if(!ggrd->ages || ! ggrd->age_time){
fprintf(stderr,"ggrd_read_vel_grids: memory error\n");
return -5;
}
/*
read in the age grids
*/
- for(ivt=0;ivt < v->nage;ivt++){
- v->ages[ivt].bandlim = v->age_bandlim;
+ for(ivt=0;ivt < ggrd->nage;ivt++){
+ ggrd->ages[ivt].bandlim = ggrd->age_bandlim;
sprintf(tfilename,"%s%i/age.grd",prefix,ivt+1);
if(ggrd_grdtrack_init_general(FALSE,tfilename,char_dummy, /* load file */
- "-Lx",(v->ages+ivt),verbose,
+ "-Lx",(ggrd->ages+ivt),verbose,
FALSE)){
fprintf(stderr,"ggrd_read_vel_grids: file error\n");
return -10;
}
- if(ivt < v->nage-1) /* assign beginning of stage as time
+ if(ivt < ggrd->nage-1) /* assign beginning of stage as time
for seafloor age */
- v->age_time[ivt] = v->thist.vtimes[ivt*3];
+ ggrd->age_time[ivt] = ggrd->time_hist.vtimes[ivt*3];
else /* end of last stage */
- v->age_time[ivt] = v->thist.vtimes[(ivt-1)*3+2];
+ ggrd->age_time[ivt] = ggrd->time_hist.vtimes[(ivt-1)*3+2];
if(verbose)
fprintf(stderr,"ggrd_read_vel_grids: read %s for seafloor age at time %g\n",
- tfilename,v->age_time[ivt]);
+ tfilename,ggrd->age_time[ivt]);
}
/* end age init */
}
@@ -194,15 +184,15 @@
// this also creates a sorting array
//
sprintf(vsfile_loc,"%s%s",prefix,GGRD_DFILE);
- ggrd_read_depth_levels(v,&index,vsfile_loc,verbose);
+ ggrd_read_depth_levels(ggrd,&index,vsfile_loc,verbose);
/*
read the velocities in binary format, either GMT grd or double bin
*/
- if(v->n[HC_R] < 1)
+ if(ggrd->v.n[HC_R] < 1)
GGRD_PE("ggrd_read_vel_grids: error: should have more than one layer, check depth file");
- if(v->read_gmt){
+ if(ggrd->v.read_gmt){
/* prepare filenames */
if(verbose)
fprintf(stderr,"ggrd_read_vel_grids: reading grd files\n");
@@ -212,24 +202,24 @@
fprintf(stderr,"ggrd_read_vel_grids: reading bin files\n");
strcpy(suffix,"bin");
}
- if(v->amode == GGRD_ONLY_VEL_STATS){
+ if(ggrd->amode == GGRD_ONLY_VEL_STATS){
sprintf(vsfile_loc,"%s.%s",prefix,GGRD_VSFILE);
fprintf(stderr,"ggrd_read_vel_grids: writing z rms_vr rms_vt rms_vp rms_vh to %s\n",
vsfile_loc);
- out = hc_open(vsfile_loc,"w","ggrd_read_vel_grids");
+ out = ggrd_open(vsfile_loc,"w","ggrd_read_vel_grids");
}
- for(ivt=0;ivt < v->thist.nvtimes;ivt++){
- if((v->history)&&(verbose))
+ for(ivt=0;ivt < ggrd->time_hist.nvtimes;ivt++){
+ if((ggrd->v.history)&&(verbose))
fprintf(stderr,"ggrd_read_vel_grids: reading velocities for time [%12g, %12g] from %3i/\n",
- v->thist.vtimes[ivt*3],
- v->thist.vtimes[ivt*3+2],ivt+1);
- for(i=0;i < v->n[HC_R];i++){
+ ggrd->time_hist.vtimes[ivt*3],
+ ggrd->time_hist.vtimes[ivt*3+2],ivt+1);
+ for(i=0;i < ggrd->v.n[HC_R];i++){
//
// determine number of grd file based on resorted arrays
//
level = index[i]+1;// level numbers should go from 1 .. N
for(j=0;j<3;j++){
- if(v->history)
+ if(ggrd->v.history)
sprintf(loc_prefix,"%i/",ivt+1);
else
sprintf(loc_prefix,"./");
@@ -243,7 +233,7 @@
else
sprintf(sname,"%s%svp.%i.%s",
prefix,loc_prefix,level,suffix);
- if(v->read_gmt){
+ if(ggrd->v.read_gmt){
#ifndef USE_GMT4 /* old */
if(GMT_cdf_read_grd_info (sname,header) == -1){
fprintf(stderr,"ggrd_read_vel_grids: error opening GMT grd file %s\n",sname);
@@ -257,7 +247,7 @@
}
#endif
}else{
- in = hc_open(sname,"r","ggrd_read_vel_grids");
+ in = ggrd_open(sname,"r","ggrd_read_vel_grids");
// read header type of information
header->node_offset=FALSE;
fread(&header->x_min, sizeof(double), 1, in);
@@ -280,16 +270,16 @@
omaxphi= LON2PHI(header->x_max-(pixelreg?header->x_inc/2.0:0.0));
maxtheta=LAT2THETA(header->y_min+(pixelreg?header->y_inc/2.0:0.0));
mintheta=LAT2THETA(header->y_max-(pixelreg?header->y_inc/2.0:0.0));
- v->dphi= DEG2RAD( header->x_inc);
- v->dtheta=DEG2RAD( header->y_inc);
+ ggrd->v.dphi= DEG2RAD( header->x_inc);
+ ggrd->v.dtheta=DEG2RAD( header->y_inc);
if(HC_DIFFERENT(minphi,0.0) ||
- HC_DIFFERENT(mintheta,v->dtheta*0.5) ||
- HC_DIFFERENT(maxtheta,GGRD_PI - v->dtheta*0.5) ||
+ HC_DIFFERENT(mintheta,ggrd->v.dtheta*0.5) ||
+ HC_DIFFERENT(maxtheta,GGRD_PI - ggrd->v.dtheta*0.5) ||
(HC_DIFFERENT(omaxphi,GGRD_TWOPI) &&
- HC_DIFFERENT(omaxphi,GGRD_TWOPI - v->dphi))){
+ HC_DIFFERENT(omaxphi,GGRD_TWOPI - ggrd->v.dphi))){
fprintf(stderr,"ggrd_read_vel_grids: expecting 0/360(or %g)/%g/%g range, problem with %s\n",
- 360-RAD2DEG(v->dphi),-90+RAD2DEG(v->dtheta*0.5),
- 90-RAD2DEG(v->dtheta*0.5),sname);
+ 360-RAD2DEG(ggrd->v.dphi),-90+RAD2DEG(ggrd->v.dtheta*0.5),
+ 90-RAD2DEG(ggrd->v.dtheta*0.5),sname);
fprintf(stderr,"ggrd_read_vel_grids: expected range in radians: t: %g/%g p: %g/%g\n",
mintheta,maxtheta,minphi,omaxphi);
fprintf(stderr,"ggrd_read_vel_grids: expected range in degrees: %g/%g/%g/%g\n",
@@ -303,75 +293,75 @@
//
// check if we should throw away double entries at 0 and 360
if(!HC_DIFFERENT(omaxphi,GGRD_TWOPI)){
- v->n[HC_PHI] = header->nx - 1;
+ ggrd->v.n[HC_PHI] = header->nx - 1;
wraparound = TRUE;
}else{
- v->n[HC_PHI] = header->nx;
+ ggrd->v.n[HC_PHI] = header->nx;
wraparound = FALSE;
}
- v->n[HC_THETA] = header->ny;
- if(HC_DIFFERENT(v->dtheta,GGRD_PI /
- ((GGRD_CPREC)(v->n[HC_THETA])))||
- HC_DIFFERENT(v->dphi,GGRD_TWOPI/
- ((GGRD_CPREC)(v->n[HC_PHI])))){
+ ggrd->v.n[HC_THETA] = header->ny;
+ if(HC_DIFFERENT(ggrd->v.dtheta,GGRD_PI /
+ ((GGRD_CPREC)(ggrd->v.n[HC_THETA])))||
+ HC_DIFFERENT(ggrd->v.dphi,GGRD_TWOPI/
+ ((GGRD_CPREC)(ggrd->v.n[HC_PHI])))){
fprintf(stderr,"ggrd_read_vel_grids: spacing error: ndx/dx phi: %g/%g theta: %g/%g\n",
- GGRD_TWOPI/v->n[HC_PHI],v->dphi,
- GGRD_PI/v->n[HC_THETA],v->dtheta);
+ GGRD_TWOPI/ggrd->v.n[HC_PHI],ggrd->v.dphi,
+ GGRD_PI/ggrd->v.n[HC_THETA],ggrd->v.dtheta);
return(-3);
}
//
// set auxiliary grid dimensions
//
- v->n[HC_TPPROD] = v->n[HC_THETA] * v->n[HC_PHI];// ny * nx
- v->n[HC_NRNTNP] = v->n[HC_TPPROD] * v->n[HC_R]; // ny * nx * nr
- os = v->n[HC_NRNTNP] * v->thist.nvtimes;// ny * nx * nr *nt
+ ggrd->v.n[HC_TPPROD] = ggrd->v.n[HC_THETA] * ggrd->v.n[HC_PHI];// ny * nx
+ ggrd->v.n[HC_NRNTNP] = ggrd->v.n[HC_TPPROD] * ggrd->v.n[HC_R]; // ny * nx * nr
+ os = ggrd->v.n[HC_NRNTNP] * ggrd->time_hist.nvtimes;// ny * nx * nr *nt
//
// allocate space
- hc_vecalloc(&v->vr,os,"ggrd_readgrds: vr");
- hc_vecalloc(&v->vt,os,"ggrd_readgrds: vt");
- hc_vecalloc(&v->vp,os,"ggrd_readgrds: vp");
- if(v->read_gmt){
+ ggrd_vecalloc(&ggrd->v.vr,os,"ggrd_readgrds: vr");
+ ggrd_vecalloc(&ggrd->v.vt,os,"ggrd_readgrds: vt");
+ ggrd_vecalloc(&ggrd->v.vp,os,"ggrd_readgrds: vp");
+ if(ggrd->v.read_gmt){
// this has to be of the original GRD file size
// NOT the new grid dimensions
fgrd = (float *)malloc(sizeof(float) * header->nx * header->ny);
}else{
dgrd = (double *)malloc(sizeof(double) * header->nx * header->ny);
}
- if((v->read_gmt && !fgrd) ||(!v->read_gmt && !dgrd))
+ if((ggrd->v.read_gmt && !fgrd) ||(!ggrd->v.read_gmt && !dgrd))
HC_MEMERROR("ggrd_read_vel_grids: velocity fields:");
if(weighted){
//
// need to construct 2-D array with area weights
//
- hc_vecalloc(&weights,v->n[HC_TPPROD],"readgrds");
+ ggrd_vecalloc(&weights,ggrd->v.n[HC_TPPROD],"readgrds");
for(theta=mintheta,
- k=0;k < v->n[HC_THETA];k++,theta += v->dtheta){
+ k=0;k < ggrd->v.n[HC_THETA];k++,theta += ggrd->v.dtheta){
tmp = sin(theta);
- for(l=0;l < v->n[HC_PHI];l++)
- weights[k*v->n[HC_PHI]+l] = tmp;
+ for(l=0;l < ggrd->v.n[HC_PHI];l++)
+ weights[k*ggrd->v.n[HC_PHI]+l] = tmp;
}
}
if(verbose)
fprintf(stderr,"ggrd_read_vel_grids: x: %g/%g/%g nx: %i y: %g/%g/%g ny: %i wrap: %i v_c: %g\n",
PHI2LON(minphi),PHI2LON(omaxphi),
- RAD2DEG(v->dphi),v->n[HC_PHI],
+ RAD2DEG(ggrd->v.dphi),ggrd->v.n[HC_PHI],
THETA2LAT(maxtheta),THETA2LAT(mintheta),
- RAD2DEG(v->dtheta),v->n[HC_THETA],wraparound,
- v->velscale);
+ RAD2DEG(ggrd->v.dtheta),ggrd->v.n[HC_THETA],wraparound,
+ ggrd->v.velscale);
init = TRUE;
}else{
if(HC_DIFFERENT(minphi,LON2PHI(header->x_min+(pixelreg?header->x_inc/2.0:0.0)))||
HC_DIFFERENT(omaxphi,LON2PHI(header->x_max-(pixelreg?header->x_inc/2.0:0.0)))||
HC_DIFFERENT(maxtheta,LAT2THETA(header->y_min+(pixelreg?header->y_inc/2.0:0.0)))||
HC_DIFFERENT(mintheta,LAT2THETA(header->y_max-(pixelreg?header->y_inc/2.0:0.0)))||
- HC_DIFFERENT(v->dphi,DEG2RAD(header->x_inc))||
- HC_DIFFERENT(v->dtheta,DEG2RAD( header->y_inc))){
+ HC_DIFFERENT(ggrd->v.dphi,DEG2RAD(header->x_inc))||
+ HC_DIFFERENT(ggrd->v.dtheta,DEG2RAD( header->y_inc))){
fprintf(stderr,"ggrd_read_vel_grids: grd files have different size, grd: %s\n",
sname);
exit(-1);
}
}
- if(v->read_gmt){
+ if(ggrd->v.read_gmt){
#ifdef USE_GMT4
// read the netcdf GRD file
sprintf(header->name,"%s",sname);
@@ -391,42 +381,42 @@
// AND: leave those pointer calculations here, since we
// do not initially have the size of the arrays
//
- os1 = v->n[HC_NRNTNP] * ivt;
- os1 += v->n[HC_TPPROD] * i;
+ os1 = ggrd->v.n[HC_NRNTNP] * ivt;
+ os1 += ggrd->v.n[HC_TPPROD] * i;
// these should theoretically be == zero
if(j == HC_R){
//
// vr
//
if(zero_boundary_vr &&
- (1.0 - v->rlevels[i] < HC_EPS_PREC)){
+ (1.0 - ggrd->v.rlevels[i] < HC_EPS_PREC)){
if(verbose)
fprintf(stderr,"ggrd_read_vel_grids: WARNING: assuming level %3i is at surface and setting vr to zero\n",
level);
- ggrd_resort_and_check((v->vr+os1),fgrd,dgrd,v->n[HC_PHI],
- v->n[HC_THETA],wraparound,1.0/v->velscale,
- v->read_gmt,TRUE,0.0,&v->vd.w_warned);
- }else if((zero_boundary_vr)&&(v->rlevels[i] < v->rcmb)){
+ ggrd_resort_and_check((ggrd->v.vr+os1),fgrd,dgrd,ggrd->v.n[HC_PHI],
+ ggrd->v.n[HC_THETA],wraparound,1.0/ggrd->v.velscale,
+ ggrd->v.read_gmt,TRUE,0.0,&ggrd->v.vd.w_warned);
+ }else if((zero_boundary_vr)&&(ggrd->v.rlevels[i] < ggrd->v.rcmb)){
if(verbose)
fprintf(stderr,"ggrd_read_vel_grids: WARNING: assuming level %3i is at CMB and setting vr to zero\n",
level);
- ggrd_resort_and_check((v->vr+os1),fgrd,dgrd,v->n[HC_PHI],
- v->n[HC_THETA],wraparound,1.0/v->velscale,
- v->read_gmt,TRUE,0.0,&v->vd.w_warned);
+ ggrd_resort_and_check((ggrd->v.vr+os1),fgrd,dgrd,ggrd->v.n[HC_PHI],
+ ggrd->v.n[HC_THETA],wraparound,1.0/ggrd->v.velscale,
+ ggrd->v.read_gmt,TRUE,0.0,&ggrd->v.vd.w_warned);
}else
- ggrd_resort_and_check((v->vr+os1),fgrd,dgrd,v->n[HC_PHI],
- v->n[HC_THETA],wraparound,1.0/v->velscale,
- v->read_gmt,FALSE,ddummy,&v->vd.w_warned);
- hc_calc_mean_and_stddev((v->vr+os1),&ddummy,v->n[HC_TPPROD],
+ ggrd_resort_and_check((ggrd->v.vr+os1),fgrd,dgrd,ggrd->v.n[HC_PHI],
+ ggrd->v.n[HC_THETA],wraparound,1.0/ggrd->v.velscale,
+ ggrd->v.read_gmt,FALSE,ddummy,&ggrd->v.vd.w_warned);
+ ggrd_calc_mean_and_stddev((ggrd->v.vr+os1),&ddummy,ggrd->v.n[HC_TPPROD],
(mean+j),(std+j),(rms+j),FALSE,weighted,weights);
}else if(j == HC_THETA){
//
// vtheta
//
- ggrd_resort_and_check((v->vt+os1),fgrd,dgrd,v->n[HC_PHI],
- v->n[HC_THETA],wraparound,1.0/v->velscale,
- v->read_gmt,FALSE,ddummy,&v->vd.w_warned);
- hc_calc_mean_and_stddev((v->vt+os1),&ddummy,v->n[HC_TPPROD],
+ ggrd_resort_and_check((ggrd->v.vt+os1),fgrd,dgrd,ggrd->v.n[HC_PHI],
+ ggrd->v.n[HC_THETA],wraparound,1.0/ggrd->v.velscale,
+ ggrd->v.read_gmt,FALSE,ddummy,&ggrd->v.vd.w_warned);
+ ggrd_calc_mean_and_stddev((ggrd->v.vt+os1),&ddummy,ggrd->v.n[HC_TPPROD],
(mean+j),(std+j),(rms+j),FALSE,weighted,weights);
}else{
//
@@ -434,49 +424,49 @@
//
if(j != HC_PHI)
GGRD_PE("ggrd_read_vel_grds: index error");
- ggrd_resort_and_check((v->vp+os1),fgrd,dgrd,v->n[HC_PHI],v->n[HC_THETA],
- wraparound,1.0/v->velscale,
- v->read_gmt,FALSE,ddummy,&v->vd.w_warned);
- hc_calc_mean_and_stddev((v->vp+os1),&ddummy,v->n[HC_TPPROD],
+ ggrd_resort_and_check((ggrd->v.vp+os1),fgrd,dgrd,ggrd->v.n[HC_PHI],ggrd->v.n[HC_THETA],
+ wraparound,1.0/ggrd->v.velscale,
+ ggrd->v.read_gmt,FALSE,ddummy,&ggrd->v.vd.w_warned);
+ ggrd_calc_mean_and_stddev((ggrd->v.vp+os1),&ddummy,ggrd->v.n[HC_TPPROD],
(mean+j),(std+j),(rms+j),FALSE,weighted,weights);
//
// and horizontal stats, put those in the 4th element
// of mean
//
- hc_calc_mean_and_stddev((v->vp+os1),(v->vt+os1),v->n[HC_TPPROD],
+ ggrd_calc_mean_and_stddev((ggrd->v.vp+os1),(ggrd->v.vt+os1),ggrd->v.n[HC_TPPROD],
(mean+3),(std+3),(rms+3),TRUE,weighted,weights);
}
}
if(verbose)
fprintf(stderr,"ggrd_read_depth_levels: %13s: l: %3i i: %3i r: %9.7f z: %9.2f %s mean/RMS: vr: %9.2e/%9.2e vt: %9.2e/%9.2e vp: %9.2e/%9.2e\n",
- sname,level,i,v->rlevels[i],
- HC_Z_DEPTH(v->rlevels[i]),
- (weighted?"weighted":"unweighted"),mean[HC_R]*v->velscale,
- rms[HC_R]*v->velscale,mean[HC_THETA]*v->velscale,
- rms[HC_THETA]*v->velscale,mean[HC_PHI]*v->velscale,
- rms[HC_PHI]*v->velscale);
- if(v->amode == GGRD_ONLY_VEL_STATS)// velocity statistics output
+ sname,level,i,ggrd->v.rlevels[i],
+ HC_Z_DEPTH(ggrd->v.rlevels[i]),
+ (weighted?"weighted":"unweighted"),mean[HC_R]*ggrd->v.velscale,
+ rms[HC_R]*ggrd->v.velscale,mean[HC_THETA]*ggrd->v.velscale,
+ rms[HC_THETA]*ggrd->v.velscale,mean[HC_PHI]*ggrd->v.velscale,
+ rms[HC_PHI]*ggrd->v.velscale);
+ if(ggrd->amode == GGRD_ONLY_VEL_STATS)// velocity statistics output
fprintf(out,"%14.5e %14.5e %14.5e %14.5e %14.5e %5i %13.5f\n",
- HC_Z_DEPTH(v->rlevels[i]),rms[HC_R]*v->velscale,
- rms[HC_THETA]*v->velscale,rms[HC_PHI]*v->velscale,
- rms[3]*v->velscale,ivt+1,
- ((v->history)?(v->thist.vtimes[ivt*3+1]):(0.0)));
+ HC_Z_DEPTH(ggrd->v.rlevels[i]),rms[HC_R]*ggrd->v.velscale,
+ rms[HC_THETA]*ggrd->v.velscale,rms[HC_PHI]*ggrd->v.velscale,
+ rms[3]*ggrd->v.velscale,ivt+1,
+ ((ggrd->v.history)?(ggrd->time_hist.vtimes[ivt*3+1]):(0.0)));
}
}
/* free sorting array */
free(index);
- if(v->read_gmt)
+ if(ggrd->v.read_gmt)
free(fgrd);
else
free(dgrd);
if(weighted)
free(weights);
- if(v->amode == GGRD_ONLY_VEL_STATS){
+ if(ggrd->amode == GGRD_ONLY_VEL_STATS){
fclose(out);
fprintf(stderr,"ggrd_read_vel_grids: exiting after printing vel stats\n");
return(0);
}
- v->init = TRUE;
+ ggrd->v.init = TRUE;
}else{
GGRD_PE("ggrd_read_vel_grds: error, already initialized");
}
@@ -563,7 +553,7 @@
and create sorting index
*/
-void ggrd_read_depth_levels(struct ggrd_vel *v,
+void ggrd_read_depth_levels(struct ggrd_master *ggrd,
int **index,char *filename,
hc_boolean verbose)
{
@@ -571,53 +561,53 @@
int i;
GGRD_CPREC *rnew;
- in = hc_open(filename,"r","ggrd_read_depth_levels");
+ in = ggrd_open(filename,"r","ggrd_read_depth_levels");
/* set counters */
- v->n[HC_R]=0;
- v->rlevels=(GGRD_CPREC *)realloc(v->rlevels,sizeof(GGRD_CPREC));
- if(!v->rlevels)
+ ggrd->v.n[HC_R]=0;
+ ggrd->v.rlevels=(GGRD_CPREC *)realloc(ggrd->v.rlevels,sizeof(GGRD_CPREC));
+ if(!ggrd->v.rlevels)
HC_MEMERROR("ggrd_read_depth_levels");
- while(fscanf(in,HC_FLT_FORMAT,(v->rlevels + v->n[HC_R]))==1){
- if(v->n[HC_R] > 1) /* test, if sorted */
- if(fabs(v->rlevels[v->n[HC_R]] - v->rlevels[v->n[HC_R]-1]) < 1e-7)
+ while(fscanf(in,HC_FLT_FORMAT,(ggrd->v.rlevels + ggrd->v.n[HC_R]))==1){
+ if(ggrd->v.n[HC_R] > 1) /* test, if sorted */
+ if(fabs(ggrd->v.rlevels[ggrd->v.n[HC_R]] - ggrd->v.rlevels[ggrd->v.n[HC_R]-1]) < 1e-7)
GGRD_PE("ggrd_read_depth_levels: error: two radii are at same level");
- if(v->rlevels[v->n[HC_R]] < 0){
+ if(ggrd->v.rlevels[ggrd->v.n[HC_R]] < 0){
/* flip sign */
- v->rlevels[v->n[HC_R]] = -v->rlevels[v->n[HC_R]];
- if((!v->rl_warned) && (verbose)){
+ ggrd->v.rlevels[ggrd->v.n[HC_R]] = -ggrd->v.rlevels[ggrd->v.n[HC_R]];
+ if((!ggrd->v.rl_warned) && (verbose)){
fprintf(stderr,"ggrd_read_depth_levels: WARNING: flipping sign of depth levels in %s\n",
GGRD_DFILE);
- v->rl_warned = TRUE;
+ ggrd->v.rl_warned = TRUE;
}
}
/* radius of levels */
- v->rlevels[v->n[HC_R]] = HC_ND_RADIUS(v->rlevels[v->n[HC_R]]);
- if((v->rlevels[v->n[HC_R]] > 1)||
- (v->rlevels[v->n[HC_R]] < GGRD_RCMB_ND)){
+ ggrd->v.rlevels[ggrd->v.n[HC_R]] = HC_ND_RADIUS(ggrd->v.rlevels[ggrd->v.n[HC_R]]);
+ if((ggrd->v.rlevels[ggrd->v.n[HC_R]] > 1)||
+ (ggrd->v.rlevels[ggrd->v.n[HC_R]] < GGRD_RCMB_ND)){
// check for above surface or below CMB
- fprintf(stderr,"ggrd_read_depth_levels: radius %g out of range\n",v->rlevels[v->n[HC_R]]);
+ fprintf(stderr,"ggrd_read_depth_levels: radius %g out of range\n",ggrd->v.rlevels[ggrd->v.n[HC_R]]);
exit(-1);
}
- v->n[HC_R]++;
- v->rlevels=(GGRD_CPREC *)realloc(v->rlevels,sizeof(GGRD_CPREC)*
- (v->n[HC_R]+1));
- if(!v->rlevels)
+ ggrd->v.n[HC_R]++;
+ ggrd->v.rlevels=(GGRD_CPREC *)realloc(ggrd->v.rlevels,sizeof(GGRD_CPREC)*
+ (ggrd->v.n[HC_R]+1));
+ if(!ggrd->v.rlevels)
HC_MEMERROR("ggrd_read_depth_levels");
}
fclose(in);
// sort and create index
- *index=(int *)malloc(sizeof(int)*v->n[HC_R]);
+ *index=(int *)malloc(sizeof(int)*ggrd->v.n[HC_R]);
if(! *index)
HC_MEMERROR("ggrd_read_depth_levels");
- hc_indexx(v->n[HC_R],(v->rlevels-1),(*index-1));
+ ggrd_indexx(ggrd->v.n[HC_R],(ggrd->v.rlevels-1),(*index-1));
// reassign
- rnew=(GGRD_CPREC *)malloc(sizeof(GGRD_CPREC)*v->n[HC_R]);
- for(i=0;i < v->n[HC_R];i++)
- rnew[i] = v->rlevels[(*index)[i]];
- for(i=0;i < v->n[HC_R];i++)
- v->rlevels[i] = rnew[i];
+ rnew=(GGRD_CPREC *)malloc(sizeof(GGRD_CPREC)*ggrd->v.n[HC_R]);
+ for(i=0;i < ggrd->v.n[HC_R];i++)
+ rnew[i] = ggrd->v.rlevels[(*index)[i]];
+ for(i=0;i < ggrd->v.n[HC_R];i++)
+ ggrd->v.rlevels[i] = rnew[i];
free(rnew);
if(verbose)
fprintf(stderr,"ggrd_read_depth_levels: read %i levels from %s, r_min: %g r_max: %g \n",
- v->n[HC_R],GGRD_DFILE,v->rlevels[0],v->rlevels[v->n[HC_R]-1]);
+ ggrd->v.n[HC_R],GGRD_DFILE,ggrd->v.rlevels[0],ggrd->v.rlevels[ggrd->v.n[HC_R]-1]);
}
Modified: mc/1D/hc/trunk/ggrd_struc.h
===================================================================
--- mc/1D/hc/trunk/ggrd_struc.h 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/ggrd_struc.h 2008-02-27 05:54:14 UTC (rev 11265)
@@ -93,37 +93,24 @@
/*
-velocity structure for interpolation
+structure for 3-D velocity interpolation
*/
struct ggrd_vel{
-
GGRD_CPREC *vr,*vt,*vp; /* velocity field */
int n[5]; /* dimensions in r, theta, and
phi directions */
int ntnp,nrntnp; /* */
GGRD_CPREC *rlevels; /* levels where velocities are
specified */
- GGRD_CPREC dtheta,dphi; /* spacing in theta and phi */
+ GGRD_CPREC dtheta,dphi; /* spacing in theta and phi */
GGRD_CPREC velscale,rcmb;
- struct ggrd_t thist;
unsigned char init, /* initialized? */
history, /* time-dependent? */
- use_age, /* use an additional age file */
read_gmt; /* read GMT grd files or binary format?*/
unsigned char rl_warned,vd_init,vd_reduce_r_stencil; /* */
- int amode;
- struct ggrd_gt *ages; /* for seafloor ages */
- int nage; /* ntime for velo + 1 */
- GGRD_CPREC *age_time; /* times */
- float age_bandlim; /* bandlim for age to decide on
- continent
- */
+
struct ggrd_vip vd; /* velocity interpolation structure */
- /* seafloor stuff */
- unsigned short sf_init;
- GGRD_CPREC sf_old_age,sf_old_f1,sf_old_f2;
- int sf_old_left,sf_old_right,sf_ntlim;
};
@@ -144,18 +131,32 @@
struct ggrd_master{ /* master structure */
int mat_control,mat_control_init;
- int vel_control,vel_control_init;
+ int vtop_control,vtop_control_init;
+ int age_control,age_control_init;
char mat_file[1000];
- char vel_file[1000];
+ char vtop_file[1000];
+ char age_file[1000];
/* grid structures */
struct ggrd_gt *mat; /* material grids */
+ /* surface velocity */
+ struct ggrd_gt *svp,*svt; /* phi/theta surface velocities */
+ /* age stuff */
+ struct ggrd_gt *ages;
+ int nage,amode;
+ GGRD_CPREC *age_time; /* times */
+ float age_bandlim; /* bandlim for age to decide on
+ continent */
- /* different for velocities */
- struct ggrd_vel *ggrd_v; /* velocity grids */
+ unsigned short sf_init; /* seafloor stuff */
+ GGRD_CPREC sf_old_age,sf_old_f1,sf_old_f2;
+ int sf_old_left,sf_old_right,sf_ntlim;
+
+ struct ggrd_vel v; /* 3D velocity grid structure */
+
+ /* time histtory */
struct ggrd_t time_hist; /* time history structure */
-
/* temperature init */
struct ggrd_temp_init temp_init;
};
Modified: mc/1D/hc/trunk/ggrd_test.c
===================================================================
--- mc/1D/hc/trunk/ggrd_test.c 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/ggrd_test.c 2008-02-27 05:54:14 UTC (rev 11265)
@@ -2,7 +2,7 @@
int main(int argc, char **argv)
{
- struct ggrd_vel *v;
+ struct ggrd_master *ggrd;
HC_PREC xloc[3],time,vr[4],vphi[4],vtheta[4],dtrange;
static int order = 3;
hc_boolean calc_derivatives ;
@@ -10,15 +10,16 @@
/*
initialize velocity structure
*/
- v = (struct ggrd_vel *)calloc(1,sizeof(struct ggrd_vel));
- ggrd_init_vstruc(v);
- v->history = TRUE; /* expect history */
- v->use_age = TRUE; /* expect seafloor age files */
- v->age_bandlim = 900;
+ ggrd = (struct ggrd_master *)calloc(1,sizeof(struct ggrd_master));
+ ggrd_init_master(ggrd);
+ ggrd->v.history = TRUE; /* expect history */
+ ggrd->age_control = TRUE; /* expect seafloor age files */
+ ggrd->age_bandlim = 900;
/*
read in velocities
*/
- if(ggrd_read_vel_grids(v,1.0,FALSE,TRUE,"/home/walter/becker/data/plates/past/clb/hall/")){
+ if(ggrd_read_vel_grids(ggrd,1.0,FALSE,TRUE,
+ "/home/walter/becker/data/plates/past/clb/hall/")){
fprintf(stderr,"error opening grids\n");
exit(-1);
}
@@ -42,11 +43,11 @@
/*
interpolate
*/
- if(ggrd_find_vel_and_der(xloc,time,dtrange,v,order,calc_derivatives,
+ if(ggrd_find_vel_and_der(xloc,time,dtrange,ggrd,order,calc_derivatives,
TRUE,vr,vtheta,vphi))
exit(-1);
if(interpolate_seafloor_ages(xloc[HC_THETA],
- xloc[HC_PHI],time,v, &age))
+ xloc[HC_PHI],time,ggrd, &age))
exit(-1);
fprintf(stdout,"%11g %11g\t%11g %11g %11g\t%11g\n",lon,lat,vphi[0],-vtheta[0],vr[0],age);
Modified: mc/1D/hc/trunk/ggrd_velinterpol.c
===================================================================
--- mc/1D/hc/trunk/ggrd_velinterpol.c 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/ggrd_velinterpol.c 2008-02-27 05:54:14 UTC (rev 11265)
@@ -30,7 +30,7 @@
int ggrd_find_vel_and_der(GGRD_CPREC *xloc,
GGRD_CPREC time,GGRD_CPREC dtrange,
- struct ggrd_vel *v,int order,
+ struct ggrd_master *ggrd,int order,
hc_boolean icalc_der,
hc_boolean verbose,
GGRD_CPREC *dvr,GGRD_CPREC *dvtheta,
@@ -45,7 +45,7 @@
struct wgt weights[3];
- if(!v->vd.init){
+ if(!ggrd->v.vd.init){
//
// do some checks if called for the first time
//
@@ -54,38 +54,38 @@
order,GGRD_MAX_ORDER);
return(-1);
}
- v->vd.old_order = order;
- v->vd.orderp1 = order+1;
- if(v->n[HC_R] < order+1){
+ ggrd->v.vd.old_order = order;
+ ggrd->v.vd.orderp1 = order+1;
+ if(ggrd->v.n[HC_R] < order+1){
if(verbose)
fprintf(stderr,"ggrd_find_vel_and_der: WARNING: reducing r stencil to nl-1: %i\n",
- v->n[HC_R] -1 );
- v->vd.reduce_r_stencil = TRUE;
+ ggrd->v.n[HC_R] -1 );
+ ggrd->v.vd.reduce_r_stencil = TRUE;
}
- if(v->n[HC_PHI] < order+1){
+ if(ggrd->v.n[HC_PHI] < order+1){
fprintf(stderr,"ggrd_find_vel_and_der: need at least four lon levels\n");
fprintf(stderr,"ggrd_find_vel_and_der: using polynomial interpolation\n");
return(-1);
}
- if(v->n[HC_THETA] < order+1){
+ if(ggrd->v.n[HC_THETA] < order+1){
fprintf(stderr,"ggrd_find_vel_and_der: need at least four lat levels\n");
fprintf(stderr,"ggrd_find_vel_and_der: using polynomial interpolation\n");
return(-1);
}
/* test levels */
- for(i=1;i < v->n[HC_R];i++){
- if(v->rlevels[i] <= v->rlevels[i-1]){
+ for(i=1;i < ggrd->v.n[HC_R];i++){
+ if(ggrd->v.rlevels[i] <= ggrd->v.rlevels[i-1]){
fprintf(stderr,"ggrd_find_vel_and_der: error:\n");
fprintf(stderr,"ggrd_find_vel_and_der: rlevels have to be ascending\n");
fprintf(stderr,"ggrd_find_vel_and_der: i: %i r(i): %g r(i-1): %g\n",i,
- v->rlevels[i],v->rlevels[i-1]);
+ ggrd->v.rlevels[i],ggrd->v.rlevels[i-1]);
return(-1);
}
}
for(i=0;i < 3;i++){
if(i==HC_R)
- lorder = (v->vd.reduce_r_stencil)?(v->n[HC_R]-1):(order);
+ lorder = (ggrd->v.vd.reduce_r_stencil)?(ggrd->v.n[HC_R]-1):(order);
else
lorder = order;
if(lorder < 0){
@@ -96,20 +96,20 @@
/* loop through dimensions */
/* all directions will be interpolated up to the
default order + 1 */
- v->vd.istencil[i] = lorder + 1;
+ ggrd->v.vd.istencil[i] = lorder + 1;
// these are offsets for the interpolation routine
- v->vd.isshift[i] = (int)(v->vd.istencil[i]/2.0);
+ ggrd->v.vd.isshift[i] = (int)(ggrd->v.vd.istencil[i]/2.0);
}
//
// this is for derivatives, initialize once as zeroes
//
for(i=0;i < 1+GGRD_MAX_IORDER*3;i++)
- v->vd.ider[i] = 0.0;
- v->vd.init = TRUE;
+ ggrd->v.vd.ider[i] = 0.0;
+ ggrd->v.vd.init = TRUE;
} /* end of init loop */
- if(order != v->vd.old_order){
+ if(order != ggrd->v.vd.old_order){
fprintf(stderr,"ggrd_find_vel_and_der: error: order (%i) shouldn't change, old: %i\n",
- order,v->vd.old_order);
+ order,ggrd->v.vd.old_order);
return(-1);
}
//
@@ -157,109 +157,109 @@
// RADIAL COMPONENT
//
// default order-1 interpolation for radial direction
- v->vd.ixtracer[HC_R] = -1;
- ilim = v->n[HC_R] - 1;
+ ggrd->v.vd.ixtracer[HC_R] = -1;
+ ilim = ggrd->v.n[HC_R] - 1;
i=0;
- while ((v->vd.ixtracer[HC_R] == -1) && (i < ilim)){
- if(xloc[HC_R] <= v->rlevels[i]){
- v->vd.ixtracer[HC_R] = i;
+ while ((ggrd->v.vd.ixtracer[HC_R] == -1) && (i < ilim)){
+ if(xloc[HC_R] <= ggrd->v.rlevels[i]){
+ ggrd->v.vd.ixtracer[HC_R] = i;
break;
}
i++;
}
- if(v->vd.ixtracer[HC_R] == -1){ // no depth levels found, tracer is above surface
- v->vd.ixtracer[HC_R] = ilim; // assign last layer, x_r should be corrected by the RK routines
+ if(ggrd->v.vd.ixtracer[HC_R] == -1){ // no depth levels found, tracer is above surface
+ ggrd->v.vd.ixtracer[HC_R] = ilim; // assign last layer, x_r should be corrected by the RK routines
}
//
// pick indices of grid points for the radial stencil
//
- for(i=0;i < v->vd.istencil[HC_R];i++)
- igrid[HC_R][i] = v->vd.ixtracer[HC_R] - v->vd.isshift[HC_R] + i;
+ for(i=0;i < ggrd->v.vd.istencil[HC_R];i++)
+ igrid[HC_R][i] = ggrd->v.vd.ixtracer[HC_R] - ggrd->v.vd.isshift[HC_R] + i;
//
// make sure all grid points exist
//
ishift = igrid[HC_R][0];
if(ishift < 0)
- for(i=0;i < v->vd.istencil[HC_R];i++)
+ for(i=0;i < ggrd->v.vd.istencil[HC_R];i++)
igrid[HC_R][i] -= ishift;
// same for upper limit
- ishift = igrid[HC_R][v->vd.istencil[HC_R]-1] - ilim;
+ ishift = igrid[HC_R][ggrd->v.vd.istencil[HC_R]-1] - ilim;
if(ishift > 0)
- for(i=0;i < v->vd.istencil[HC_R];i++)
+ for(i=0;i < ggrd->v.vd.istencil[HC_R];i++)
igrid[HC_R][i] -= ishift;
// find values of r for each grid point
- for(j=HC_R*(v->vd.orderp1),i=0;i < v->vd.istencil[HC_R];i++)
- grid[j+i] = v->rlevels[igrid[HC_R][i]];
+ for(j=HC_R*(ggrd->v.vd.orderp1),i=0;i < ggrd->v.vd.istencil[HC_R];i++)
+ grid[j+i] = ggrd->v.rlevels[igrid[HC_R][i]];
//
// THETA COMPONENT
//
- ilim = v->n[HC_THETA]-1;
- v->vd.ixtracer[HC_THETA] = (int)(xloc[HC_THETA]/v->dtheta);
+ ilim = ggrd->v.n[HC_THETA]-1;
+ ggrd->v.vd.ixtracer[HC_THETA] = (int)(xloc[HC_THETA]/ggrd->v.dtheta);
// pick grid points
- for(i=0;i < v->vd.istencil[HC_THETA];i++)
- igrid[HC_THETA][i] = v->vd.ixtracer[HC_THETA] - v->vd.isshift[HC_THETA]+i;
+ for(i=0;i < ggrd->v.vd.istencil[HC_THETA];i++)
+ igrid[HC_THETA][i] = ggrd->v.vd.ixtracer[HC_THETA] - ggrd->v.vd.isshift[HC_THETA]+i;
//
// adust grid points to avoid wrap-around
//
ishift = igrid[HC_THETA][0];
if(ishift < 0){
- for(i=0;i < v->vd.istencil[HC_THETA];i++)
+ for(i=0;i < ggrd->v.vd.istencil[HC_THETA];i++)
igrid[HC_THETA][i] -= ishift;
}
// same for upper limit
- ishift = igrid[HC_THETA][v->vd.istencil[HC_THETA]-1] - ilim;
+ ishift = igrid[HC_THETA][ggrd->v.vd.istencil[HC_THETA]-1] - ilim;
if(ishift > 0)
- for(i=0;i < v->vd.istencil[HC_THETA];i++)
+ for(i=0;i < ggrd->v.vd.istencil[HC_THETA];i++)
igrid[HC_THETA][i] -= ishift;
//
// find values of theta: since given on dtheta/2 .... pi-dtheta/2
// theta_i = (i+0.5)*dtheta, i=0,1,...
//
- for(j=HC_THETA*(v->vd.orderp1),i=0;i < v->vd.istencil[HC_THETA];i++)
- grid[j+i] = (igrid[HC_THETA][i] + 0.5) * v->dtheta;
+ for(j=HC_THETA*(ggrd->v.vd.orderp1),i=0;i < ggrd->v.vd.istencil[HC_THETA];i++)
+ grid[j+i] = (igrid[HC_THETA][i] + 0.5) * ggrd->v.dtheta;
//
// now for phi
//
- ilim = v->n[HC_PHI]-1;
- v->vd.ixtracer[HC_PHI] = (int)(xloc[HC_PHI]/v->dphi+.5);
+ ilim = ggrd->v.n[HC_PHI]-1;
+ ggrd->v.vd.ixtracer[HC_PHI] = (int)(xloc[HC_PHI]/ggrd->v.dphi+.5);
//pick grid points
- for(i=0;i < v->vd.istencil[HC_PHI];i++){
- igrid[HC_PHI][i] = v->vd.ixtracer[HC_PHI] - v->vd.isshift[HC_PHI] + i;
+ for(i=0;i < ggrd->v.vd.istencil[HC_PHI];i++){
+ igrid[HC_PHI][i] = ggrd->v.vd.ixtracer[HC_PHI] - ggrd->v.vd.isshift[HC_PHI] + i;
//
// wrap around
//
if(igrid[HC_PHI][i] > ilim)
- igrid[HC_PHI][i] -= v->n[HC_PHI];
+ igrid[HC_PHI][i] -= ggrd->v.n[HC_PHI];
if(igrid[HC_PHI][i] < 0)
- igrid[HC_PHI][i] += v->n[HC_PHI];
+ igrid[HC_PHI][i] += ggrd->v.n[HC_PHI];
}
//
// find values of phi. phi_i = i * dphi
//
- for(j=HC_PHI*(v->vd.orderp1),i=0;i < v->vd.istencil[HC_PHI];i++)
- grid[j+i] = igrid[HC_PHI][i] * v->dphi;
+ for(j=HC_PHI*(ggrd->v.vd.orderp1),i=0;i < ggrd->v.vd.istencil[HC_PHI];i++)
+ grid[j+i] = igrid[HC_PHI][i] * ggrd->v.dphi;
#ifdef DEBUG
//
// check if all indices are ok
//
- for(i=0;i < v->vd.istencil[HC_R];i++){
- if((igrid[HC_R][i]< 0)||(igrid[HC_R][i] >= v->n[HC_R])){
+ for(i=0;i < ggrd->v.vd.istencil[HC_R];i++){
+ if((igrid[HC_R][i]< 0)||(igrid[HC_R][i] >= ggrd->v.n[HC_R])){
fprintf(stderr,"ggrd_find_vel_and_der: row %i r index %i error\n",i,igrid[HC_R][i]);
return(-1);
}
}
- for(i=0;i < v->vd.istencil[HC_THETA];i++){
- if((igrid[HC_THETA][i] < 0) || (igrid[HC_THETA][i] >= v->n[HC_THETA])){
+ for(i=0;i < ggrd->v.vd.istencil[HC_THETA];i++){
+ if((igrid[HC_THETA][i] < 0) || (igrid[HC_THETA][i] >= ggrd->v.n[HC_THETA])){
fprintf(stderr,"ggrd_find_vel_and_der: row %i theta index %i error \n",i,igrid[HC_THETA][i]);
return(-1);
}
}
- for(i=0;i < v->vd.istencil[HC_PHI];i++){
- if((igrid[HC_PHI][i] < 0)||(igrid[HC_PHI][i] >= v->n[HC_PHI])){
+ for(i=0;i < ggrd->v.vd.istencil[HC_PHI];i++){
+ if((igrid[HC_PHI][i] < 0)||(igrid[HC_PHI][i] >= ggrd->v.n[HC_PHI])){
fprintf(stderr,"ggrd_find_vel_and_der: row %i phi index %i error\n",
i,igrid[HC_PHI][i]);
return(-1);
@@ -272,20 +272,20 @@
if(verbose >= 2){ /* debugging output */
fprintf(stderr,"ggrd_velinterpol: x={%g, %g, %g} [%i (%i), %i (%i), %i(%i)]\n",
- xloc[HC_R],xloc[HC_THETA],xloc[HC_PHI],v->vd.ixtracer[HC_R],v->n[HC_R],
- v->vd.ixtracer[HC_THETA],v->n[HC_THETA],
- v->vd.ixtracer[HC_PHI],v->n[HC_PHI]);
+ xloc[HC_R],xloc[HC_THETA],xloc[HC_PHI],ggrd->v.vd.ixtracer[HC_R],ggrd->v.n[HC_R],
+ ggrd->v.vd.ixtracer[HC_THETA],ggrd->v.n[HC_THETA],
+ ggrd->v.vd.ixtracer[HC_PHI],ggrd->v.n[HC_PHI]);
for(i=0;i < 3;i++){
rnet = GGRD_TWOPI;
fprintf(stderr,"ggrd_velinterpol: dim: %i:",i);
- for(j=0;j < v->vd.istencil[i];j++){
- fprintf(stderr,"%.5f (%3i) ",grid[i*(v->vd.orderp1)+j],igrid[i][j]);
+ for(j=0;j < ggrd->v.vd.istencil[i];j++){
+ fprintf(stderr,"%.5f (%3i) ",grid[i*(ggrd->v.vd.orderp1)+j],igrid[i][j]);
/* find min distance to stencil point */
- vrloc = fabs(grid[i*(v->vd.orderp1)+j]-xloc[i]);
+ vrloc = fabs(grid[i*(ggrd->v.vd.orderp1)+j]-xloc[i]);
if(vrloc < rnet){rnet=vrloc;k=j;}
}
fprintf(stderr,"\tms: %.4f(%i)\n",
- (GGRD_CPREC)k/(GGRD_CPREC)(v->vd.istencil[i]-1),v->vd.istencil[i]);
+ (GGRD_CPREC)k/(GGRD_CPREC)(ggrd->v.vd.istencil[i]-1),ggrd->v.vd.istencil[i]);
}
}
#endif
@@ -295,7 +295,7 @@
// compute all the weights for each stencil
//
for(i=0;i < 3;i++){ /* loop through spatial dimension */
- ggrd_weights(xloc[i],(grid+i*(v->vd.orderp1)),v->vd.istencil[i],iorder,weights[i].w);
+ ggrd_weights(xloc[i],(grid+i*(ggrd->v.vd.orderp1)),ggrd->v.vd.istencil[i],iorder,weights[i].w);
}
//
// first calculate velocities only (idindex=0) or vel and derivatives
@@ -306,17 +306,17 @@
first velocities
*/
dvr[0] = dvtheta[0] = dvphi[0] = 0.0;
- for(i=0;i < v->vd.istencil[HC_R];i++){ // radial
- for(j=0; j < v->vd.istencil[HC_THETA];j++){ // theta
- for(k=0; k < v->vd.istencil[HC_PHI];k++){ // phi
+ for(i=0;i < ggrd->v.vd.istencil[HC_R];i++){ // radial
+ for(j=0; j < ggrd->v.vd.istencil[HC_THETA];j++){ // theta
+ for(k=0; k < ggrd->v.vd.istencil[HC_PHI];k++){ // phi
rnet = weights[HC_R].w[i][0];
rnet *= weights[HC_THETA].w[j][0];
rnet *= weights[HC_PHI].w[k][0];
- index = igrid[HC_R][i] * v->n[HC_TPPROD] +
- igrid[HC_THETA][j] * v->n[HC_PHI] +
+ index = igrid[HC_R][i] * ggrd->v.n[HC_TPPROD] +
+ igrid[HC_THETA][j] * ggrd->v.n[HC_PHI] +
igrid[HC_PHI][k];
ggrd_get_velocities(&vrloc,&vthetaloc,&vphiloc,index,
- v,time,dtrange);
+ ggrd,time,dtrange);
dvr[0] += rnet * vrloc;
dvtheta[0] += rnet * vthetaloc;
dvphi[0] += rnet * vphiloc;
@@ -333,25 +333,25 @@
// this is the derivative yes/no array
// set once, and switch off again below//
//
- v->vd.ider[m] = 1;
- for(i=0;i < v->vd.istencil[HC_R];i++){
- for(j=0; j < v->vd.istencil[HC_THETA];j++){
- for(k=0; k < v->vd.istencil[HC_PHI];k++){
- rnet = weights[HC_R].w[i][v->vd.ider[HC_R+1]];
- rnet *= weights[HC_THETA].w[j][v->vd.ider[HC_THETA+1]];
- rnet *= weights[HC_PHI].w[k][v->vd.ider[HC_PHI+1]];
- index = igrid[HC_R][i] * v->n[HC_TPPROD] +
- igrid[HC_THETA][j] * v->n[HC_PHI] +
+ ggrd->v.vd.ider[m] = 1;
+ for(i=0;i < ggrd->v.vd.istencil[HC_R];i++){
+ for(j=0; j < ggrd->v.vd.istencil[HC_THETA];j++){
+ for(k=0; k < ggrd->v.vd.istencil[HC_PHI];k++){
+ rnet = weights[HC_R].w[i][ggrd->v.vd.ider[HC_R+1]];
+ rnet *= weights[HC_THETA].w[j][ggrd->v.vd.ider[HC_THETA+1]];
+ rnet *= weights[HC_PHI].w[k][ggrd->v.vd.ider[HC_PHI+1]];
+ index = igrid[HC_R][i] * ggrd->v.n[HC_TPPROD] +
+ igrid[HC_THETA][j] * ggrd->v.n[HC_PHI] +
igrid[HC_PHI][k];
ggrd_get_velocities(&vrloc,&vthetaloc,&vphiloc,index,
- v,time,dtrange);
+ ggrd,time,dtrange);
dvr[m] += rnet * vrloc;
dvtheta[m] += rnet * vthetaloc;
dvphi[m] += rnet * vphiloc;
}
}
}
- v->vd.ider[m]=0; // reset derivative switxh to zero
+ ggrd->v.vd.ider[m]=0; // reset derivative switxh to zero
}
/* succesful return */
return 0;
@@ -370,36 +370,36 @@
//
void ggrd_get_velocities(GGRD_CPREC *vrloc,GGRD_CPREC *vthetaloc,
GGRD_CPREC *vphiloc,
- int index, struct ggrd_vel *v,
+ int index, struct ggrd_master *ggrd,
GGRD_CPREC time,GGRD_CPREC dtrange)
{
int index1,i1,i2;
GGRD_CPREC vf1,vf2;
- if((index < 0) || (index >= v->n[HC_NRNTNP]) ){
+ if((index < 0) || (index >= ggrd->v.n[HC_NRNTNP]) ){
HC_ERROR("ggrd_get_velocities","index out of bounds");
exit(-1);
}
- if(v->thist.nvtimes == 1){
+ if(ggrd->time_hist.nvtimes == 1){
// only one time-step, steady-state calculation
- *vrloc= v->vr[index];
- *vthetaloc = v->vt[index];
- *vphiloc= v->vp[index];
+ *vrloc= ggrd->v.vr[index];
+ *vthetaloc = ggrd->v.vt[index];
+ *vphiloc= ggrd->v.vp[index];
} else {
// interpolate in time
- ggrd_interpol_time(time,&v->thist,&i1,&i2,&vf1,&vf2,dtrange);
+ ggrd_interpol_time(time,&ggrd->time_hist,&i1,&i2,&vf1,&vf2,dtrange);
if(fabs(vf1) > 1e-7){
- index1 = i1 * v->n[HC_NRNTNP] + index;
- *vrloc= v->vr[index1] * vf1 ;
- *vthetaloc = v->vt[index1] * vf1;
- *vphiloc= v->vp[index1] * vf1 ;
+ index1 = i1 * ggrd->v.n[HC_NRNTNP] + index;
+ *vrloc= ggrd->v.vr[index1] * vf1 ;
+ *vthetaloc = ggrd->v.vt[index1] * vf1;
+ *vphiloc= ggrd->v.vp[index1] * vf1 ;
}else{
*vrloc = *vthetaloc = *vphiloc = 0.0;
}
if(fabs(vf2) > 1e-7){
- index1 = i2 * v->n[HC_NRNTNP] + index;
- *vrloc += v->vr[index1] * vf2;
- *vthetaloc += v->vt[index1] * vf2;
- *vphiloc += v->vp[index1] * vf2;
+ index1 = i2 * ggrd->v.n[HC_NRNTNP] + index;
+ *vrloc += ggrd->v.vr[index1] * vf2;
+ *vthetaloc += ggrd->v.vt[index1] * vf2;
+ *vphiloc += ggrd->v.vp[index1] * vf2;
}
}
}
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2008-02-27 05:54:14 UTC (rev 11265)
@@ -11,20 +11,20 @@
int ggrd_init_thist_from_file(struct ggrd_t *, char *, unsigned char, unsigned char);
void ggrd_gt_interpolate_z(double, float *, int, int *, int *, double *, double *, unsigned char, unsigned char *);
void ggrd_interpol_time(double, struct ggrd_t *, int *, int *, double *, double *, double);
-int interpolate_seafloor_ages(double, double, double, struct ggrd_vel *, double *);
+int interpolate_seafloor_ages(double, double, double, struct ggrd_master *, double *);
void ggrd_vecalloc(double **, int, char *);
void ggrd_vecrealloc(double **, int, char *);
float ggrd_gt_rms(float *, int);
float ggrd_gt_mean(float *, int);
/* ggrd_readgrds.c */
-void ggrd_init_vstruc(struct ggrd_vel *);
-int ggrd_read_vel_grids(struct ggrd_vel *, double, unsigned short, unsigned short, char *);
+void ggrd_init_vstruc(struct ggrd_master *);
+int ggrd_read_vel_grids(struct ggrd_master *, double, unsigned short, unsigned short, char *);
void ggrd_resort_and_check(double *, float *, double *, int, int, unsigned short, double, unsigned short, unsigned short, double, unsigned char *);
-void ggrd_read_depth_levels(struct ggrd_vel *, int **, char *, unsigned short);
+void ggrd_read_depth_levels(struct ggrd_master *, int **, char *, unsigned short);
/* ggrd_test.c */
/* ggrd_velinterpol.c */
-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);
+int ggrd_find_vel_and_der(double *, double, double, struct ggrd_master *, int, unsigned short, unsigned short, double *, double *, double *);
+void ggrd_get_velocities(double *, double *, double *, int, struct ggrd_master *, double, double);
void ggrd_weights(double, double *, int, int, double [(5 +1)][(1 +1)]);
/* grdinttester.c */
/* hc_extract_sh_layer.c */
Modified: mc/1D/hc/trunk/hc_extract_sh_layer.c
===================================================================
--- mc/1D/hc/trunk/hc_extract_sh_layer.c 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/hc_extract_sh_layer.c 2008-02-27 05:54:14 UTC (rev 11265)
@@ -63,7 +63,7 @@
/*
read in solution
*/
- in = hc_open(argv[1],"r","hc_extract_sh_layer");
+ in = ggrd_open(argv[1],"r","hc_extract_sh_layer");
hc_read_sh_solution(model,&sol,in,binary,verbose);
fclose(in);
nsol = model->nradp2 * 3;
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/hc_init.c 2008-02-27 05:54:14 UTC (rev 11265)
@@ -432,7 +432,7 @@
from bottom to top
*/
- in = hc_open(filename,"r","hc_assign_viscosity");
+ in = ggrd_open(filename,"r","hc_assign_viscosity");
hc_vecrealloc(&hc->rvisc,1,"hc_assign_viscosity");
hc_vecrealloc(&hc->visc,1,"hc_assign_viscosity");
hc->nvis = 0;mean = 0.0;mws = 0.0;
@@ -557,7 +557,7 @@
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");
+ in = ggrd_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");
@@ -590,7 +590,7 @@
*/
- in = hc_open(filename,"r","hc_assign_density");
+ in = ggrd_open(filename,"r","hc_assign_density");
if(verbose)
fprintf(stderr,"hc_assign_density: reading density anomalies in [%g%%] from %s\n",
100*HC_DENSITY_SCALING,filename);
@@ -884,7 +884,7 @@
if(verbose)
fprintf(stderr,"hc_assign_plate_velocities: expecting [cm/yr] pol/tor from %s\n",
filename);
- in = hc_open(filename,"r","hc_assign_plate_velocities");
+ in = ggrd_open(filename,"r","hc_assign_plate_velocities");
if(!sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer, &nset,
&zlabel,&ivec,in,FALSE,
pvel_in_binary,verbose)){
Modified: mc/1D/hc/trunk/hc_misc.c
===================================================================
--- mc/1D/hc/trunk/hc_misc.c 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/hc_misc.c 2008-02-27 05:54:14 UTC (rev 11265)
@@ -12,20 +12,6 @@
-/*
- open a file safely and give an error message if there was
- a problem
-*/
-FILE *hc_open(char *name, char *mode, char *program)
-{
- FILE *in;
- if((in=fopen(name,mode)) == NULL){
- fprintf(stderr,"%s: error: can not open file %s for mode %s access\n",
- program,name,mode);
- exit(-1);
- }
- return in;
-}
/* double vector allocation */
void hc_dvecalloc(double **x,int n,char *message)
{
@@ -208,150 +194,3 @@
}
*i += 1;
}
-/*
-
- calculate mean and standard deviation of x_i
- if hypoth is set, calc mean and stddev of sqrt(x_i^2 + y_i^2)
- if weighted is set, uses weights[n] for weighting the mean and so on
-
- this way of calculating the stddev is inaccurate but fast
-
-*/
-void hc_calc_mean_and_stddev(HC_PREC *x, HC_PREC *y,int n,HC_PREC *mean,
- HC_PREC *stddev,HC_PREC *rms,
- hc_boolean hypoth,
- hc_boolean weighted,HC_PREC *weight)
-{
- HC_PREC sum1=0.0,sum2=0.0,tmp,ws;
- int i;
- if(n <= 1){
- fprintf(stderr,"hc_calc_mean_and_stddev: error: n: %i\n",n);
- exit(-1);
- }
- ws=0.0;
- if(hypoth){// sqrt(x^2+y+2)
- for(i=0;i<n;i++){
- if(weighted){
- tmp = hypot(x[i],y[i]) * weight[i];
- ws += weight[i];
- }else{
- tmp = hypot(x[i],y[i]);
- ws += 1.0;
- }
- sum1 += tmp;sum2 += tmp * tmp;
- }
- }else{
- for(i=0;i<n;i++){
- if(weighted){
- tmp = x[i] * weight[i];
- ws += weight[i];
- }else{
- tmp = x[i];
- ws += 1.0;
- }
- sum1 += tmp;sum2 += tmp*tmp;
- }
- }
- // standard deviation
- tmp = (ws * sum2 - sum1 * sum1) / (ws*(ws-1.0));
- if(tmp > 0)
- *stddev = sqrt(tmp);
- else
- *stddev = 0.0;
- *rms = sqrt(sum2 / ws);// RMS
- *mean = sum1 / ws; // mean
-}
-/*
-
- sort array by order, returns index
-
- GIVE INPUT INDEX SHIFTED
- (ie. x-1)
- output will be 0...n-1
-
- based on numerical recipes
-
-
-*/
-
-
-#define HC_INDEXX_SWAP(a,b) itemp=(a);(a)=(b);(b)=itemp;
-#define HC_INDEXX_M 7
-#define HC_INDEXX_NSTACK 5000
-
-void hc_indexx(int n,HC_PREC *arr, int *indx)
-{
- int i,indxt,ir=n,itemp,j,k,l=1;
- int jstack=0,*istack;
- HC_PREC a;
- //
- istack=(int *)malloc(sizeof(int)*(HC_INDEXX_NSTACK+1));
- if(!istack)
- HC_MEMERROR("indexx");
- //
- for (j=1;j<=n;j++)
- indx[j]=j;
- for (;;) {
- if (ir-l < HC_INDEXX_M) {
- for (j=l+1;j<=ir;j++) {
- indxt=indx[j];
- a=arr[indxt];
- for (i=j-1;i>=1;i--) {
- if (arr[indx[i]] <= a) break;
- indx[i+1]=indx[i];
- }
- indx[i+1]=indxt;
- }
- if (jstack == 0) break;
- ir=istack[jstack--];
- l=istack[jstack--];
- } else {
- k=(l+ir) >> 1;
- HC_INDEXX_SWAP(indx[k],indx[l+1]);
- if (arr[indx[l+1]] > arr[indx[ir]]) {
- HC_INDEXX_SWAP(indx[l+1],indx[ir]);
- }
- if (arr[indx[l]] > arr[indx[ir]]) {
- HC_INDEXX_SWAP(indx[l],indx[ir]);
- }
- if (arr[indx[l+1]] > arr[indx[l]]) {
- HC_INDEXX_SWAP(indx[l+1],indx[l]);
- }
- i=l+1;
- j=ir;
- indxt=indx[l];
- a=arr[indxt];
- for (;;) {
- do i++; while (arr[indx[i]] < a);
- do j--; while (arr[indx[j]] > a);
- if (j < i) break;
- HC_INDEXX_SWAP(indx[i],indx[j]);
- }
- indx[l]=indx[j];
- indx[j]=indxt;
- jstack += 2;
- if (jstack > HC_INDEXX_NSTACK) {
- fprintf(stderr,
- "indexx: HC_INDEXX_NSTACK too small");
- exit(-1);
- }
- if (ir-i+1 >= j-l) {
- istack[jstack]=ir;
- istack[jstack-1]=i;
- ir=j-1;
- } else {
- istack[jstack]=j-1;
- istack[jstack-1]=l;
- l=i;
- }
- }
- }
- free(istack);
- // go back to normal numbering
- for(i=1;i<=n;i++)
- indx[i]--;
-}
-#undef HC_INDEXX_M
-#undef HC_INDEXX_NSTACK
-#undef HC_INDEXX_SWAP
-
Modified: mc/1D/hc/trunk/hc_output.c
===================================================================
--- mc/1D/hc/trunk/hc_output.c 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/hc_output.c 2008-02-27 05:54:14 UTC (rev 11265)
@@ -139,7 +139,7 @@
1,verbose);
/* depth file */
- dout = hc_open(dfilename,"w","hc_print_spatial_solution");
+ dout = ggrd_open(dfilename,"w","hc_print_spatial_solution");
if(verbose >= 2)
fprintf(stderr,"hc_print_spatial_solution: writing depth levels to %s\n",
dfilename);
@@ -172,7 +172,7 @@
if(binary){
/* binary output */
sprintf(filename,"%s.%i.bin",name,i+1);
- out = hc_open(filename,"w","hc_print_spatial_solution");
+ out = ggrd_open(filename,"w","hc_print_spatial_solution");
for(j=los=0;j < np;j++,los+=2){ /* loop through all points in layer */
fwrite((xy+los),sizeof(float),2,out);
for(k=0;k<3;k++)
@@ -184,7 +184,7 @@
}else{
/* ASCII output */
sprintf(filename,"%s.%i.dat",name,i+1);
- out = hc_open(filename,"w","hc_print_spatial_solution");
+ out = ggrd_open(filename,"w","hc_print_spatial_solution");
for(j=los=0;j < np;j++,los+=2){ /* loop through all points in layer */
for(k=0;k<3;k++)
value[k] = sol_x[os[k]] * fac[k];
@@ -341,7 +341,7 @@
/* number of output layers */
nl = hc->nrad + 2;
- out = hc_open(filename,"w","hc_print_poloidal_solution");
+ out = ggrd_open(filename,"w","hc_print_poloidal_solution");
for(l=1;l <= ll;l++){
for(m=0;m <= l;m++){
alim = (m==0)?(1):(2);
@@ -380,7 +380,7 @@
if(verbose)
fprintf(stderr,"hc_print_toiroidal_solution: writing toroidal solutions 1 and 2 as f(l,r) to %s\n",
filename);
- out = hc_open(filename,"w","hc_toroidal_solution");
+ out = ggrd_open(filename,"w","hc_toroidal_solution");
for(l=1;l <= ll;l++){
for(os=i=0;i < nl;i++,os+=lmaxp1)
fprintf(out,"%3i %16.7e %16.7e %16.7e\n",
Modified: mc/1D/hc/trunk/main.c
===================================================================
--- mc/1D/hc/trunk/main.c 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/main.c 2008-02-27 05:54:14 UTC (rev 11265)
@@ -167,7 +167,7 @@
if(p->verbose)
fprintf(stderr,"%s: writing spherical harmonics solution to %s\n",
argv[0],filename);
- out = hc_open(filename,"w","main");
+ out = ggrd_open(filename,"w","main");
hc_print_spectral_solution(model,sol_spectral,out,
p->solution_mode,
p->sol_binary_out,p->verbose);
@@ -179,7 +179,7 @@
sprintf(filename,"%s",HC_GEOID_FILE);
if(p->verbose)
fprintf(stderr,"%s: writing geoid to %s\n",argv[0],filename);
- out = hc_open(filename,"w","main");
+ out = ggrd_open(filename,"w","main");
hc_print_sh_scalar_field(geoid,out,FALSE,FALSE,p->verbose);
fclose(out);
}
Modified: mc/1D/hc/trunk/sh_test.c
===================================================================
--- mc/1D/hc/trunk/sh_test.c 2008-02-27 05:42:49 UTC (rev 11264)
+++ mc/1D/hc/trunk/sh_test.c 2008-02-27 05:54:14 UTC (rev 11265)
@@ -73,7 +73,7 @@
spherical harmonics routines operate
*/
sprintf(filename,"hc.%i.dat",type);
- out = hc_open(filename,"w",argv[0]);
+ out = ggrd_open(filename,"w",argv[0]);
sh_print_model_spatial_basis(&model,out,verbose);
fclose(out);
if(verbose)
@@ -87,7 +87,7 @@
read in data at pixel locations , allocates odata array
*/
sprintf(filename,"test.%i.data",type);
- in = hc_open(filename,"r",argv[0]);
+ in = ggrd_open(filename,"r",argv[0]);
sh_read_model_spatial_data(&model,&model.data,in,verbose);
fclose(in);
@@ -108,7 +108,7 @@
output of spatial basis of that expansion
*/
sprintf(filename,"testn.%i.data",type);
- out = hc_open(filename,"w",argv[0]);
+ out = ggrd_open(filename,"w",argv[0]);
sh_print_model_spatial_data(&model,ndata,out,verbose);
fclose(out);
for(i=0;i<3;i++){
More information about the cig-commits
mailing list