[cig-commits] r8015 - mc/3D/CitcomS/trunk/lib

[tan2 at geodynamics.org]

Author: tan2
Date: 2007-09-24 14:28:05 -0700 (Mon, 24 Sep 2007)
New Revision: 8015

Modified:
   mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
Log:
Multiply rho*cp to DT/dt

Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c	2007-09-24 19:55:35 UTC (rev 8014)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c	2007-09-24 21:28:05 UTC (rev 8015)
@@ -546,7 +546,7 @@
     double prod,sfn;
     struct Shape_function1 GM;
     struct Shape_function1_dA dGamma;
-    double temp,rho,heating;
+    double temp,rho,cp,heating;
     int nz;
 
     void get_global_1d_shape_fn();
@@ -609,6 +609,7 @@
 
     nz = ((el-1) % E->lmesh.elz) + 1;
     rho = 0.5 * (E->refstate.rho[nz] + E->refstate.rho[nz+1]);
+    cp = 0.5 * (E->refstate.heat_capacity[nz] + E->refstate.heat_capacity[nz+1]);
 
     if(E->control.disptn_number == 0)
         heating = rho * Q;
@@ -625,7 +626,8 @@
 	for(i=1;i<=vpts;i++)
 	  Eres[j] -=
 	    PG.vpt[GNVINDEX(j,i)] * dOmega.vpt[i]
-              * (dT[i] - heating + v1[i]*tx1[i] + v2[i]*tx2[i] + v3[i]*tx3[i])
+              * ((dT[i] + v1[i]*tx1[i] + v2[i]*tx2[i] + v3[i]*tx3[i])*rho*cp
+                 - heating )
               + diff * dOmega.vpt[i] * E->heating_latent[m][el]
               * (GNx.vpt[GNVXINDEX(0,j,i)]*tx1[i]*rtf[3][i] +
                  GNx.vpt[GNVXINDEX(1,j,i)]*tx2[i]*sint[i] +