[cig-commits] r5342 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Nov 22 11:41:44 PST 2006
Author: tan2
Date: 2006-11-22 11:41:43 -0800 (Wed, 22 Nov 2006)
New Revision: 5342
Modified:
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Output_h5.c
mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c
mc/3D/CitcomS/trunk/lib/Topo_gravity.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Refactored geoid calculation
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2006-11-22 02:06:37 UTC (rev 5341)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2006-11-22 19:41:43 UTC (rev 5342)
@@ -249,9 +249,11 @@
FILE *fp1;
compute_geoid(E, E->sphere.harm_geoid,
- E->sphere.harm_tpgt, E->sphere.harm_tpgb);
+ E->sphere.harm_geoid_from_bncy,
+ E->sphere.harm_geoid_from_tpgt,
+ E->sphere.harm_geoid_from_tpgb);
- if (E->parallel.me == 0) {
+ if (E->parallel.me == (E->parallel.nprocz-1)) {
sprintf(output_file, "%s.geoid.%d.%d", E->control.data_file,
E->parallel.me, cycles);
fp1 = output_open(output_file);
@@ -264,10 +266,15 @@
for (ll=0; ll<=E->output.llmax; ll++)
for(mm=0; mm<=ll; mm++) {
p = E->sphere.hindex[ll][mm];
- fprintf(fp1,"%d %d %.4e %.4e %.4e %.4e\n",
+ fprintf(fp1,"%d %d %.4e %.4e %.4e %.4e %.4e %.4e\n",
ll, mm,
- E->sphere.harm_geoid[0][p], E->sphere.harm_geoid[1][p],
- E->sphere.harm_tpgt[0][p], E->sphere.harm_tpgt[1][p]);
+ E->sphere.harm_geoid[0][p],
+ E->sphere.harm_geoid[1][p],
+ E->sphere.harm_geoid_from_tpgt[0][p],
+ E->sphere.harm_geoid_from_tpgt[1][p],
+ E->sphere.harm_geoid_from_bncy[0][p],
+ E->sphere.harm_geoid_from_bncy[1][p]);
+
}
fclose(fp1);
Modified: mc/3D/CitcomS/trunk/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_h5.c 2006-11-22 02:06:37 UTC (rev 5341)
+++ mc/3D/CitcomS/trunk/lib/Output_h5.c 2006-11-22 19:41:43 UTC (rev 5342)
@@ -136,7 +136,8 @@
extern void parallel_process_termination();
extern void heat_flux(struct All_variables *);
extern void get_STD_topo(struct All_variables *, float**, float**, float**, float**, int);
-extern void compute_geoid(struct All_variables *, float**, float**, float**);
+extern void compute_geoid(struct All_variables *, float**, float**,
+ float**, float**);
/****************************************************************************
@@ -188,9 +189,9 @@
/* allocate buffer */
if (E->output.stress == 1)
- E->hdf5.data = (float *)malloc((tensor3d->n) * sizeof(float));
+ E->hdf5.data = (float *)malloc((tensor3d->n) * sizeof(float));
else
- E->hdf5.data = (float *)malloc((vector3d->n) * sizeof(float));
+ E->hdf5.data = (float *)malloc((vector3d->n) * sizeof(float));
/* reuse buffer */
tensor3d->data = E->hdf5.data;
@@ -1041,10 +1042,12 @@
{
int ll;
int mm;
- float geoid_sin;
- float geoid_cos;
+ float total_sin;
+ float total_cos;
float tpgt_sin;
float tpgt_cos;
+ float bncy_sin;
+ float bncy_cos;
} *row;
@@ -1065,11 +1068,11 @@
H5T_NATIVE_INT);
status = H5Tinsert(datatype, "order", HOFFSET(struct HDF5_GEOID, mm),
H5T_NATIVE_INT);
- status = H5Tinsert(datatype, "geoid_sin",
- HOFFSET(struct HDF5_GEOID, geoid_sin),
+ status = H5Tinsert(datatype, "total_sin",
+ HOFFSET(struct HDF5_GEOID, total_sin),
H5T_NATIVE_FLOAT);
- status = H5Tinsert(datatype, "geoid_cos",
- HOFFSET(struct HDF5_GEOID, geoid_cos),
+ status = H5Tinsert(datatype, "total_cos",
+ HOFFSET(struct HDF5_GEOID, total_cos),
H5T_NATIVE_FLOAT);
status = H5Tinsert(datatype, "tpgt_sin",
HOFFSET(struct HDF5_GEOID, tpgt_sin),
@@ -1077,6 +1080,12 @@
status = H5Tinsert(datatype, "tpgt_cos",
HOFFSET(struct HDF5_GEOID, tpgt_cos),
H5T_NATIVE_FLOAT);
+ status = H5Tinsert(datatype, "bncy_sin",
+ HOFFSET(struct HDF5_GEOID, bncy_sin),
+ H5T_NATIVE_FLOAT);
+ status = H5Tinsert(datatype, "bncy_cos",
+ HOFFSET(struct HDF5_GEOID, bncy_cos),
+ H5T_NATIVE_FLOAT);
/* Create the dataspace */
dataspace = H5Screate_simple(rank, &dim, NULL);
@@ -1098,10 +1107,12 @@
set_attribute_string(dataset, "FIELD_0_NAME", "degree");
set_attribute_string(dataset, "FIELD_1_NAME", "order");
- set_attribute_string(dataset, "FIELD_2_NAME", "geoid_sin");
- set_attribute_string(dataset, "FIELD_3_NAME", "geoid_cos");
+ set_attribute_string(dataset, "FIELD_2_NAME", "total_sin");
+ set_attribute_string(dataset, "FIELD_3_NAME", "total_cos");
set_attribute_string(dataset, "FIELD_4_NAME", "tpgt_sin");
set_attribute_string(dataset, "FIELD_5_NAME", "tpgt_cos");
+ set_attribute_string(dataset, "FIELD_6_NAME", "bncy_sin");
+ set_attribute_string(dataset, "FIELD_7_NAME", "bncy_cos");
set_attribute_double(dataset, "FIELD_0_FILL", 0);
set_attribute_double(dataset, "FIELD_1_FILL", 0);
@@ -1109,13 +1120,17 @@
set_attribute_double(dataset, "FIELD_3_FILL", 0);
set_attribute_double(dataset, "FIELD_4_FILL", 0);
set_attribute_double(dataset, "FIELD_5_FILL", 0);
+ set_attribute_double(dataset, "FIELD_6_FILL", 0);
+ set_attribute_double(dataset, "FIELD_7_FILL", 0);
/* Create property list for independent dataset write */
dxpl_id = H5Pcreate(H5P_DATASET_XFER);
status = H5Pset_dxpl_mpio(dxpl_id, H5FD_MPIO_INDEPENDENT);
compute_geoid(E, E->sphere.harm_geoid,
- E->sphere.harm_tpgt, E->sphere.harm_tpgb);
+ E->sphere.harm_geoid_from_bncy,
+ E->sphere.harm_geoid_from_tpgt,
+ E->sphere.harm_geoid_from_tpgb);
if (E->parallel.me == 0) {
/* Prepare data */
@@ -1126,10 +1141,12 @@
for(mm = 0; mm <= ll; mm++) {
row[i].ll = ll;
row[i].mm = mm;
- row[i].geoid_sin = E->sphere.harm_geoid[0][i];
- row[i].geoid_cos = E->sphere.harm_geoid[1][i];
- row[i].tpgt_sin = E->sphere.harm_tpgt[0][i];
- row[i].tpgt_cos = E->sphere.harm_tpgt[1][i];
+ row[i].total_sin = E->sphere.harm_geoid[0][i];
+ row[i].total_cos = E->sphere.harm_geoid[1][i];
+ row[i].tpgt_sin = E->sphere.harm_geoid_from_tpgt[0][i];
+ row[i].tpgt_cos = E->sphere.harm_geoid_from_tpgt[1][i];
+ row[i].bncy_sin = E->sphere.harm_geoid_from_bncy[0][i];
+ row[i].bncy_cos = E->sphere.harm_geoid_from_bncy[1][i];
i ++;
}
Modified: mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c 2006-11-22 02:06:37 UTC (rev 5341)
+++ mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c 2006-11-22 19:41:43 UTC (rev 5342)
@@ -20,19 +20,20 @@
i=0;
for (ll=0;ll<=E->output.llmax;ll++)
- for (mm=0;mm<=ll;mm++) {
- E->sphere.hindex[ll][mm] = i;
- i++;
- }
+ for (mm=0;mm<=ll;mm++) {
+ E->sphere.hindex[ll][mm] = i;
+ i++;
+ }
E->sphere.hindice = i;
/* spherical harmonic coeff (0=cos, 1=sin)
for surface topo, cmb topo and geoid */
for (i=0;i<=1;i++) {
- E->sphere.harm_tpgt[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
- E->sphere.harm_tpgb[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
- E->sphere.harm_geoid[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
+ E->sphere.harm_geoid[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
+ E->sphere.harm_geoid_from_bncy[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
+ E->sphere.harm_geoid_from_tpgt[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
+ E->sphere.harm_geoid_from_tpgb[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
}
compute_sphereh_table(E);
@@ -55,40 +56,40 @@
struct Shape_function1_dA dGamma;
for (i=0;i<E->sphere.hindice;i++) {
- sphc[i] = 0.0;
- sphs[i] = 0.0;
+ sphc[i] = 0.0;
+ sphs[i] = 0.0;
}
for (m=1;m<=E->sphere.caps_per_proc;m++)
- for (es=1;es<=E->lmesh.snel;es++) {
- el = es*E->lmesh.elz;
+ for (es=1;es<=E->lmesh.snel;es++) {
+ el = es*E->lmesh.elz;
- get_global_1d_shape_fn(E,el,&M,&dGamma,1,m);
+ get_global_1d_shape_fn(E,el,&M,&dGamma,1,m);
- for (ll=0;ll<=E->output.llmax;ll++)
- for (mm=0; mm<=ll; mm++) {
+ for (ll=0;ll<=E->output.llmax;ll++)
+ for (mm=0; mm<=ll; mm++) {
- p = E->sphere.hindex[ll][mm];
+ p = E->sphere.hindex[ll][mm];
- for(nint=1;nint<=onedvpoints[E->mesh.nsd];nint++) {
- for(d=1;d<=onedvpoints[E->mesh.nsd];d++) {
- j = E->sien[m][es].node[d];
- sphc[p] += TG[m][E->sien[m][es].node[d]]
- * E->sphere.tablesplm[m][j][p]
- * E->sphere.tablescosf[m][j][mm]
- * E->M.vpt[GMVINDEX(d,nint)]
- * dGamma.vpt[GMVGAMMA(1,nint)];
- sphs[p] += TG[m][E->sien[m][es].node[d]]
- * E->sphere.tablesplm[m][j][p]
- * E->sphere.tablessinf[m][j][mm]
- * E->M.vpt[GMVINDEX(d,nint)]
- * dGamma.vpt[GMVGAMMA(1,nint)];
- }
- }
+ for(nint=1;nint<=onedvpoints[E->mesh.nsd];nint++) {
+ for(d=1;d<=onedvpoints[E->mesh.nsd];d++) {
+ j = E->sien[m][es].node[d];
+ sphc[p] += TG[m][E->sien[m][es].node[d]]
+ * E->sphere.tablesplm[m][j][p]
+ * E->sphere.tablescosf[m][j][mm]
+ * E->M.vpt[GMVINDEX(d,nint)]
+ * dGamma.vpt[GMVGAMMA(1,nint)];
+ sphs[p] += TG[m][E->sien[m][es].node[d]]
+ * E->sphere.tablesplm[m][j][p]
+ * E->sphere.tablessinf[m][j][mm]
+ * E->M.vpt[GMVINDEX(d,nint)]
+ * dGamma.vpt[GMVGAMMA(1,nint)];
+ }
+ }
- } /* end for ll and mm */
+ } /* end for ll and mm */
- }
+ }
sum_across_surf_sph1(E,sphc,sphs);
@@ -108,33 +109,33 @@
for(m=1;m<=E->sphere.caps_per_proc;m++) {
- E->sphere.tablesplm[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
- E->sphere.tablescosf[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
- E->sphere.tablessinf[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
+ E->sphere.tablesplm[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
+ E->sphere.tablescosf[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
+ E->sphere.tablessinf[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
- for (i=1;i<=E->lmesh.nsf;i++) {
- E->sphere.tablesplm[m][i]= (double *)malloc((E->sphere.hindice+3)*sizeof(double));
- E->sphere.tablescosf[m][i]= (double *)malloc((E->output.llmax+3)*sizeof(double));
- E->sphere.tablessinf[m][i]= (double *)malloc((E->output.llmax+3)*sizeof(double));
- }
+ for (i=1;i<=E->lmesh.nsf;i++) {
+ E->sphere.tablesplm[m][i]= (double *)malloc((E->sphere.hindice+3)*sizeof(double));
+ E->sphere.tablescosf[m][i]= (double *)malloc((E->output.llmax+3)*sizeof(double));
+ E->sphere.tablessinf[m][i]= (double *)malloc((E->output.llmax+3)*sizeof(double));
+ }
}
for(m=1;m<=E->sphere.caps_per_proc;m++) {
- for (j=1;j<=E->lmesh.nsf;j++) {
- node = j*E->lmesh.noz;
- f=E->sx[m][2][node];
- t=E->sx[m][1][node];
- for (mm=0;mm<=E->output.llmax;mm++) {
- E->sphere.tablescosf[m][j][mm] = cos( (double)(mm)*f );
- E->sphere.tablessinf[m][j][mm] = sin( (double)(mm)*f );
- }
+ for (j=1;j<=E->lmesh.nsf;j++) {
+ node = j*E->lmesh.noz;
+ f=E->sx[m][2][node];
+ t=E->sx[m][1][node];
+ for (mm=0;mm<=E->output.llmax;mm++) {
+ E->sphere.tablescosf[m][j][mm] = cos( (double)(mm)*f );
+ E->sphere.tablessinf[m][j][mm] = sin( (double)(mm)*f );
+ }
- for (ll=0;ll<=E->output.llmax;ll++)
- for (mm=0;mm<=ll;mm++) {
- p = E->sphere.hindex[ll][mm];
- E->sphere.tablesplm[m][j][p] = modified_plgndr_a(ll,mm,t) ;
- }
- }
+ for (ll=0;ll<=E->output.llmax;ll++)
+ for (mm=0;mm<=ll;mm++) {
+ p = E->sphere.hindex[ll][mm];
+ E->sphere.tablesplm[m][j][p] = modified_plgndr_a(ll,mm,t) ;
+ }
+ }
}
return;
@@ -154,46 +155,46 @@
void parallel_process_termination();
for (m=1;m<=E->sphere.caps_per_proc;m++) {
- ia[1] = 1;
- ia[2] = E->lmesh.noz*E->lmesh.nox-E->lmesh.noz+1;
- ia[3] = E->lmesh.nno-E->lmesh.noz+1;
- ia[4] = ia[3]-E->lmesh.noz*(E->lmesh.nox-1);
+ ia[1] = 1;
+ ia[2] = E->lmesh.noz*E->lmesh.nox-E->lmesh.noz+1;
+ ia[3] = E->lmesh.nno-E->lmesh.noz+1;
+ ia[4] = ia[3]-E->lmesh.noz*(E->lmesh.nox-1);
- for (i=1;i<=4;i++) {
- xx[1][i] = E->x[m][1][ia[i]]/E->sx[m][3][ia[1]];
- xx[2][i] = E->x[m][2][ia[i]]/E->sx[m][3][ia[1]];
- xx[3][i] = E->x[m][3][ia[i]]/E->sx[m][3][ia[1]];
- }
+ for (i=1;i<=4;i++) {
+ xx[1][i] = E->x[m][1][ia[i]]/E->sx[m][3][ia[1]];
+ xx[2][i] = E->x[m][2][ia[i]]/E->sx[m][3][ia[1]];
+ xx[3][i] = E->x[m][3][ia[i]]/E->sx[m][3][ia[1]];
+ }
- get_angle_sphere_cap(xx,angle);
+ get_angle_sphere_cap(xx,angle);
- for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
- E->sphere.angle[m][i] = angle[i];
+ for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
+ E->sphere.angle[m][i] = angle[i];
- E->sphere.area[m] = area_sphere_cap(angle);
+ E->sphere.area[m] = area_sphere_cap(angle);
- for (lev=E->mesh.levmax;lev>=E->mesh.levmin;lev--)
- for (es=1;es<=E->lmesh.SNEL[lev];es++) {
- el = (es-1)*E->lmesh.ELZ[lev]+1;
- for (i=1;i<=4;i++)
- ia[i] = E->IEN[lev][m][el].node[i];
+ for (lev=E->mesh.levmax;lev>=E->mesh.levmin;lev--)
+ for (es=1;es<=E->lmesh.SNEL[lev];es++) {
+ el = (es-1)*E->lmesh.ELZ[lev]+1;
+ for (i=1;i<=4;i++)
+ ia[i] = E->IEN[lev][m][el].node[i];
- for (i=1;i<=4;i++) {
- xx[1][i] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[2][i] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[3][i] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
- }
+ for (i=1;i<=4;i++) {
+ xx[1][i] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[2][i] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[3][i] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
+ }
- get_angle_sphere_cap(xx,angle);
+ get_angle_sphere_cap(xx,angle);
- for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
- E->sphere.angle1[lev][m][i][es] = angle[i];
+ for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
+ E->sphere.angle1[lev][m][i][es] = angle[i];
- E->sphere.area1[lev][m][es] = area_sphere_cap(angle);
+ E->sphere.area1[lev][m][es] = area_sphere_cap(angle);
-/* fprintf(E->fp_out,"lev%d %d %.6e %.6e %.6e %.6e %.6e\n",lev,es,angle[1],angle[2],angle[3],angle[4],E->sphere.area1[lev][m][es]); */
+/* fprintf(E->fp_out,"lev%d %d %.6e %.6e %.6e %.6e %.6e\n",lev,es,angle[1],angle[2],angle[3],angle[4],E->sphere.area1[lev][m][es]); */
- } /* end for lev and es */
+ } /* end for lev and es */
} /* end for m */
@@ -266,22 +267,22 @@
double xx[4],angle[5],angle1[5];
for (i=1;i<=4;i++)
- ia[i] = E->IEN[lev][m][el].node[i];
+ ia[i] = E->IEN[lev][m][el].node[i];
es = (el-1)/E->lmesh.ELZ[lev]+1;
for (i=1;i<=4;i++) {
- xx[1] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[2] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[3] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
- angle[i] = get_angle(x,xx); /* get angle bet (i,j) and other four*/
- angle1[i]= E->sphere.angle1[lev][m][i][es];
+ xx[1] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[2] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[3] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
+ angle[i] = get_angle(x,xx); /* get angle bet (i,j) and other four*/
+ angle1[i]= E->sphere.angle1[lev][m][i][es];
}
area1 = area_of_sphere_triag(angle[1],angle[2],angle1[1])
- + area_of_sphere_triag(angle[2],angle[3],angle1[2])
- + area_of_sphere_triag(angle[3],angle[4],angle1[3])
- + area_of_sphere_triag(angle[4],angle[1],angle1[4]);
+ + area_of_sphere_triag(angle[2],angle[3],angle1[2])
+ + area_of_sphere_triag(angle[3],angle[4],angle1[3])
+ + area_of_sphere_triag(angle[4],angle[1],angle1[4]);
return (area1);
}
@@ -298,17 +299,17 @@
double y1[4],y2[4],get_angle();;
for (i=1;i<=4;i++) { /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
- for (j=1;j<=3;j++) {
- ii=(i==4)?1:(i+1);
- y1[j] = xx[j][i];
- y2[j] = xx[j][ii];
- }
- angle[i] = get_angle(y1,y2);
+ for (j=1;j<=3;j++) {
+ ii=(i==4)?1:(i+1);
+ y1[j] = xx[j][i];
+ y2[j] = xx[j][ii];
+ }
+ angle[i] = get_angle(y1,y2);
}
for (j=1;j<=3;j++) {
- y1[j] = xx[j][1];
- y2[j] = xx[j][3];
+ y1[j] = xx[j][1];
+ y2[j] = xx[j][3];
}
angle[5] = get_angle(y1,y2); /* angle5 for betw 1 and 3: diagonal */
@@ -326,8 +327,8 @@
const double two = 2.0;
dist=sqrt( (x[1]-xx[1])*(x[1]-xx[1])
- + (x[2]-xx[2])*(x[2]-xx[2])
- + (x[3]-xx[3])*(x[3]-xx[3]) )*pt5;
+ + (x[2]-xx[2])*(x[2]-xx[2])
+ + (x[3]-xx[3])*(x[3]-xx[3]) )*pt5;
angle = asin(dist)*two;
return (angle);
Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2006-11-22 02:06:37 UTC (rev 5341)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2006-11-22 19:41:43 UTC (rev 5342)
@@ -404,21 +404,27 @@
/* ===================================================================
=================================================================== */
-void compute_geoid(E, harm_geoid, harm_tpgt, harm_tpgb)
+void compute_geoid(E, harm_geoid, harm_geoid_from_bncy,
+ harm_geoid_from_tpgt, harm_geoid_from_tpgb)
+
struct All_variables *E;
- float *harm_geoid[2], *harm_tpgt[2], *harm_tpgb[2];
+ float *harm_geoid[2], *harm_geoid_from_bncy[2];
+ float *harm_geoid_from_tpgt[2], *harm_geoid_from_tpgb[2];
{
- int ll, mm, p;
+ int i, p;
- for (ll=1;ll<=E->output.llmax;ll++)
- for (mm=0;mm<=ll;mm++) {
- p = E->sphere.hindex[ll][mm];
- harm_geoid[0][p] = 0.0;
- harm_geoid[1][p] = 0.0;
- }
+ geoid_from_buoyancy(E, harm_geoid_from_bncy);
+ geoid_from_topography(E, harm_geoid_from_tpgt, harm_geoid_from_tpgb);
- geoid_from_buoyancy(E, harm_geoid);
- geoid_from_topography(E, harm_geoid, harm_tpgt, harm_tpgb);
+ if (E->parallel.me == (E->parallel.nprocz-1)) {
+ for (i = 0; i < 2; i++)
+ for (p = 0; p <= E->sphere.hindice; p++) {
+ harm_geoid[i][p] = harm_geoid_from_bncy[i][p]
+ + harm_geoid_from_tpgt[i][p]
+ + harm_geoid_from_tpgb[i][p];
+ }
+ }
+
}
@@ -493,30 +499,27 @@
-static void geoid_from_topography(E, harm_geoid, harm_tpgt, harm_tpgb)
+static void geoid_from_topography(E, geoid_tpgt, geoid_tpgb)
struct All_variables *E;
- float *harm_geoid[2], *harm_tpgt[2], *harm_tpgb[2];
+ float *geoid_tpgt[2], *geoid_tpgb[2];
{
-
- float *geoid[2],con1,con2,scaling,den_contrast1,den_contrast2,stress_scaling,topo_scaling1,topo_scaling2;
+ float con1,con2,scaling,den_contrast1,den_contrast2,stress_scaling,topo_scaling1,topo_scaling2;
int i,j,k,ll,mm,s;
+ float *tpgt[2], *tpgb[2];
void sum_across_depth_sph1();
void sphere_expansion();
- geoid[0] = (float *)malloc((E->sphere.hindice+2)*sizeof(float));
- geoid[1] = (float *)malloc((E->sphere.hindice+2)*sizeof(float));
-
- for (i=0;i<E->sphere.hindice;i++) {
- geoid[0][i] = 0.0;
- geoid[1][i] = 0.0;
+ for (i=0; i<2; i++) {
+ tpgt[i] = (float *)malloc((E->sphere.hindice+2)*sizeof(float));
+ tpgb[i] = (float *)malloc((E->sphere.hindice+2)*sizeof(float));
}
stress_scaling = E->data.ref_viscosity*E->data.therm_diff/
(E->data.radius_km*E->data.radius_km*1e6);
- /* density across surface */
+ /* density contrast across surface */
den_contrast1 = (E->data.density-E->data.density_above);
- /* density across CMB */
+ /* density contrast across CMB */
den_contrast2 = (E->data.density_below-E->data.density);
/* scale for surface and CMB topo */
@@ -527,96 +530,64 @@
scaling = 1.0e3*4.0*M_PI*E->data.radius_km*
E->data.grav_const/E->data.grav_acc;
+ if (E->parallel.me_loc[3] == E->parallel.nprocz-1) {
+ /* expand surface topography into sph. harm. */
+ sphere_expansion(E, E->slice.tpg, tpgt[0], tpgt[1]);
- if (E->parallel.me_loc[3]==E->parallel.nprocz-1) {
- /* dimensionalize surface topography and
- expand into spherical harmonics */
- for(j=1;j<=E->sphere.caps_per_proc;j++)
- for(i=1;i<=E->lmesh.nsf;i++) {
- E->slice.tpg[j][i] = E->slice.tpg[j][i]*topo_scaling1;
+ /* dimensionalize surface topography */
+ for (j=0; j<2; j++)
+ for (i=0; i<E->sphere.hindice; i++) {
+ tpgt[j][i] *= topo_scaling1;
}
- sphere_expansion(E,E->slice.tpg,E->sphere.harm_tpgt[0],
- E->sphere.harm_tpgt[1]);
- }
- //print_field_spectral_regular(E,TL,E->sphere.harm_tpgt[0],E->sphere.harm_tpgt[1],E->parallel.nprocz-1,0,"tpgt");
-
-
-
- if (E->parallel.me_loc[3]==0) {
- /* dimensionalize CMB topography and
- expand into spherical harmonics */
- for(j=1;j<=E->sphere.caps_per_proc;j++)
- for(i=1;i<=E->lmesh.nsf;i++) {
- E->slice.tpgb[j][i] = E->slice.tpgb[j][i]*topo_scaling2;
+ /* compute geoid due to surface topo, skip degree-0 and 1 term */
+ for (j=0; j<2; j++)
+ for (ll=2; ll<=E->output.llmax; ll++) {
+ con1 = den_contrast1 * scaling / (2.0*ll + 1.0);
+ for (mm=0; mm<=ll; mm++) {
+ i = E->sphere.hindex[ll][mm];
+ geoid_tpgt[j][i] = tpgt[j][i] * con1;
}
- sphere_expansion(E,E->slice.tpgb,E->sphere.harm_tpgb[0],
- E->sphere.harm_tpgb[1]);
+ }
}
- //print_field_spectral_regular(E,TL,E->sphere.harm_tpgb[0],E->sphere.harm_tpgb[1],0,0,"tpgb");
+ if (E->parallel.me_loc[3] == 0) {
+ /* expand bottom topography into sph. harm. */
+ sphere_expansion(E, E->slice.tpgb, tpgb[0], tpgb[1]);
-
-
- if (E->parallel.me_loc[3] == E->parallel.nprocz-1)
- for (ll=2;ll<=E->output.llmax;ll++) {
- con1 = scaling/(2.0*ll+1.0);
- for (mm=0;mm<=ll;mm++) {
- i = E->sphere.hindex[ll][mm];
- geoid[0][i]+= E->sphere.harm_tpgt[0][i]*con1*den_contrast1;
- geoid[1][i]+= E->sphere.harm_tpgt[1][i]*con1*den_contrast1;
+ /* dimensionalize bottom topography */
+ for (j=0; j<2; j++)
+ for (i=0; i<E->sphere.hindice; i++) {
+ tpgb[j][i] *= topo_scaling2;
}
- }
-
-
-
-
- if (E->parallel.me_loc[3] == 0)
- for (ll=2;ll<=E->output.llmax;ll++) {
- con1 = scaling/(2.0*ll+1.0);
- con2 = pow(E->sphere.ri,((double)(ll+2)));
- for (mm=0;mm<=ll;mm++) {
- i = E->sphere.hindex[ll][mm];
- geoid[0][i]+= E->sphere.harm_tpgb[0][i]*con1*con2*den_contrast2;
- geoid[1][i]+= E->sphere.harm_tpgb[1][i]*con1*con2*den_contrast2;
+ /* compute geoid due to bottom topo, skip degree-0 and 1 term */
+ for (j=0; j<2; j++)
+ for (ll=2; ll<=E->output.llmax; ll++) {
+ con1 = den_contrast2 * scaling / (2.0*ll + 1.0);
+ con2 = con1 * pow(E->sphere.ri, ((double)(ll+2)));
+ for (mm=0; mm<=ll; mm++) {
+ i = E->sphere.hindex[ll][mm];
+ geoid_tpgb[j][i] = tpgb[j][i] * con2;
}
}
+ }
+
/* accumulate geoid to the surface (top processors) */
- sum_across_depth_sph1(E,geoid[0],geoid[1]);
+ sum_across_depth_sph1(E, geoid_tpgb[0], geoid_tpgb[1]);
- // add to the geoid coeff
- for (ll=2;ll<=E->output.llmax;ll++)
- for (mm=0;mm<=ll;mm++) {
- i = E->sphere.hindex[ll][mm];
- harm_geoid[0][i]+= geoid[0][i];
- harm_geoid[1][i]+= geoid[1][i];
- }
+ for (j=0; j<2; j++) {
+ free ((void *) tpgt[j]);
+ free ((void *) tpgb[j]);
+ }
-
-#if 0
- E->sphere.harm_geoid[0][0]=0.0;
- E->sphere.harm_geoid[1][0]=0.0;
- E->sphere.harm_geoid[0][1]=0.0;
- E->sphere.harm_geoid[1][1]=0.0;
- E->sphere.harm_geoid[0][2]=0.0;
- E->sphere.harm_geoid[1][2]=0.0; /* zero the 00 10 11 */
-
- inv_sphere_harmonics(E,E->sphere.harm_geoid[0],E->sphere.harm_geoid[1],TL,E->parallel.nprocz-1);
-
- print_field_spectral_regular(E,TL,E->sphere.harm_geoid[0],E->sphere.harm_geoid[1],E->parallel.nprocz-1,1,"geoid");
-#endif
-
-
-
- free ((void *)geoid[0]);
- free ((void *)geoid[1]);
-
return;
}
+
+
/* ===================================================================
Consistent boundary flux method for stress ... Zhong,Gurnis,Hulbert
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2006-11-22 02:06:37 UTC (rev 5341)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2006-11-22 19:41:43 UTC (rev 5342)
@@ -378,9 +378,10 @@
int max_connections;
int *hindex[100];
int hindice;
- float *harm_tpgt[2];
- float *harm_tpgb[2];
float *harm_geoid[2];
+ float *harm_geoid_from_bncy[2];
+ float *harm_geoid_from_tpgt[2];
+ float *harm_geoid_from_tpgb[2];
double **tablenplm;
double **tablencosf;
More information about the cig-commits
mailing list