[cig-commits] r6759 - mc/3D/CitcomS/branches/compressible/lib

[tan2 at geodynamics.org]

Author: tan2
Date: 2007-05-02 17:02:06 -0700 (Wed, 02 May 2007)
New Revision: 6759

Modified:
   mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c
Log:
Incorporating the effect of density/expansivity/gravity profiles when computing geoid

Modified: mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c	2007-05-02 20:40:50 UTC (rev 6758)
+++ mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c	2007-05-03 00:02:06 UTC (rev 6759)
@@ -56,8 +56,6 @@
    if (E->parallel.me_loc[3]==E->parallel.nprocz-1)
      get_surf_stress(E,SXX,SYY,SZZ,SXY,SXZ,SZY);
 
-/*    topo_scaling1=1.0/((E->data.density-E->data.density_above)*E->data.grav_acc); */
-/*    topo_scaling2=1.0/((E->data.density_below-E->data.density)*E->data.grav_acc); */
 
    topo_scaling1 = topo_scaling2 = 1.0;
 
@@ -330,12 +328,6 @@
   (E->exchange_node_f)(E,divv,lev);
   (E->exchange_node_f)(E,vorv,lev);
 
-  /*    stress_scaling = 1.0e-6*E->data.ref_viscosity*E->data.therm_diff/ */
-  /*                       (E->data.radius_km*E->data.radius_km); */
-
-  /*    velo_scaling = 100.*365.*24.*3600.*1.0e-3*E->data.therm_diff/E->data.radius_km; */
-  /* cm/yr */
-
   stress_scaling = velo_scaling = 1.0;
 
   for(m=1;m<=E->sphere.caps_per_proc;m++)
@@ -455,6 +447,7 @@
 {
     int m,k,ll,mm,node,i,j,p,noz,snode;
     float *TT[NCS],radius,*geoid[2],dlayer,con1,con2,scaling2,scaling;
+    double buoy2rho;
     void sphere_expansion();
     void sum_across_depth_sph1();
 
@@ -475,13 +468,15 @@
 
     /* loop over each layer */
     for(k=1;k<E->lmesh.noz;k++)  {
-
+        buoy2rho = scaling2 * E->refstate.rho[k]
+            * E->refstate.expansivity[k] / E->refstate.gravity[k];
         for(m=1;m<=E->sphere.caps_per_proc;m++)
             for(i=1;i<=E->lmesh.noy;i++)
                 for(j=1;j<=E->lmesh.nox;j++)  {
                     node=k+(j-1)*E->lmesh.noz+(i-1)*E->lmesh.nox*E->lmesh.noz;
                     p = j + (i-1)*E->lmesh.nox;
-                    TT[m][p] = (E->buoyancy[m][node]+E->buoyancy[m][node+1])*0.5*scaling2;
+                    TT[m][p] = (E->buoyancy[m][node]+E->buoyancy[m][node+1])
+                        * 0.5 * buoy2rho;
                 }
 
         /* expand TT into spherical harmonics */
@@ -539,18 +534,14 @@
         (E->data.radius_km*E->data.radius_km*1e6);
 
     /* density contrast across surface */
-    den_contrast1 = (E->data.density-E->data.density_above);
+    den_contrast1 = (E->data.density-E->data.density_above)*E->refstate.rho[E->lmesh.noz];
     /* density contrast across CMB */
-    den_contrast2 = (E->data.density_below-E->data.density);
+    den_contrast2 = (E->data.density_below-E->data.density)*E->refstate.rho[1];
 
     /* scale for surface and CMB topo */
-    topo_scaling1 = stress_scaling / (den_contrast1*E->data.grav_acc);
-    topo_scaling2 = stress_scaling / (den_contrast2*E->data.grav_acc);
+    topo_scaling1 = stress_scaling / (den_contrast1*E->data.grav_acc*E->refstate.gravity[E->lmesh.noz]);
+    topo_scaling2 = stress_scaling / (den_contrast2*E->data.grav_acc*E->refstate.gravity[1]);
 
-    /* scale for geoid */
-    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]);
@@ -561,6 +552,10 @@
                 tpgt[j][i] *= topo_scaling1;
             }
 
+        /* scale for geoid */
+        scaling = 1.0e3 * 4.0 * M_PI * E->data.radius_km * E->data.grav_const
+            / (E->data.grav_acc * E->refstate.gravity[E->lmesh.noz]);
+
         /* 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++)   {
@@ -582,6 +577,11 @@
             for (i=0; i<E->sphere.hindice; i++) {
                 tpgb[j][i] *= topo_scaling2;
             }
+
+        /* scale for geoid */
+        scaling = 1.0e3 * 4.0 * M_PI * E->data.radius_km * E->data.grav_const
+            / (E->data.grav_acc * E->refstate.gravity[1]);
+
         /* 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++)   {