[cig-commits] r11177 - in mc/3D/CitcomS/trunk: bin lib

becker at geodynamics.org becker at geodynamics.org
Sat Feb 16 16:14:07 PST 2008


Author: becker
Date: 2008-02-16 16:14:06 -0800 (Sat, 16 Feb 2008)
New Revision: 11177

Modified:
   mc/3D/CitcomS/trunk/bin/Citcom.c
   mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c
   mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
   mc/3D/CitcomS/trunk/lib/Initial_temperature.c
   mc/3D/CitcomS/trunk/lib/Instructions.c
   mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
   mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
   mc/3D/CitcomS/trunk/lib/convection_variables.h
   mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Modified Citcom.c such that E->mat gets assigned based on element
depth (as was default before), even if mat_control is not selected.

Added (back in?) ggrd netcdf type VIP prefactor material control via
ggrd_mat_control flags (as documented in Instructions.c)

Renamed some of the ggrd stuff to make it more modular.





Modified: mc/3D/CitcomS/trunk/bin/Citcom.c
===================================================================
--- mc/3D/CitcomS/trunk/bin/Citcom.c	2008-02-17 00:11:30 UTC (rev 11176)
+++ mc/3D/CitcomS/trunk/bin/Citcom.c	2008-02-17 00:14:06 UTC (rev 11177)
@@ -210,13 +210,18 @@
 	output_checkpoint(E);
     }
 
-    if(E->control.mat_control==1)
-      read_mat_from_file(E);
-    /*
-      else
-      construct_mat_group(E);
+    /* 
+       make sure that E->mat always gets initialized as a function of depth
+       the previous version would not call construct_mat_group if E->control.mat_control != 1
+       
+       i intend to restore default behavior here
+       
+       TWB 
+
     */
+    initialize_material(E);
 
+
     if(E->control.vbcs_file==1)
       read_velocity_boundary_from_file(E);
     /*

Modified: mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c	2008-02-17 00:11:30 UTC (rev 11176)
+++ mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c	2008-02-17 00:14:06 UTC (rev 11177)
@@ -144,7 +144,9 @@
 	break;
 
       case 3:  /* read element materials */
-
+#ifdef USE_GGRD
+	if(!E->control.ggrd.mat_control){
+#endif
 	sprintf(output_file1,"%s%0.0f.%d",E->control.mat_file,newage1,cap);
 	sprintf(output_file2,"%s%0.0f.%d",E->control.mat_file,newage2,cap);
 	fp1=fopen(output_file1,"r");
@@ -167,7 +169,9 @@
 	  else
 	    fprintf(E->fp,"Mat: File2 = No file inputted (negative age)\n");
 	}
-
+#ifdef USE_GGRD
+	}
+#endif
 	break;
 
       } /* end switch */
@@ -239,6 +243,11 @@
 
       case 3:  /* read element materials */
 
+#ifdef USE_GGRD
+	if(E->control.ggrd.mat_control){ /* use netcdf grids */
+	  ggrd_read_mat_from_file(E, 1);
+	}else{
+#endif
         VIP1 = (float*) malloc ((emax+1)*sizeof(float));
         VIP2 = (float*) malloc ((emax+1)*sizeof(float));
         LL1 = (int*) malloc ((emax+1)*sizeof(int));
@@ -281,7 +290,9 @@
          free ((void *) VIP2);
          free ((void *) LL1);
          free ((void *) LL2);
-
+#ifdef USE_GGRD
+	} /* end of branch if allowing for ggrd handling */
+#endif
 	break;
 
       } /* end switch */

Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c	2008-02-17 00:11:30 UTC (rev 11176)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c	2008-02-17 00:14:06 UTC (rev 11177)
@@ -37,6 +37,7 @@
 #include "parsing.h"
 #include "parallel_related.h"
 #include "composition_related.h"
+#include "element_definitions.h"
 
 #ifdef USE_GGRD
 
@@ -48,6 +49,9 @@
 void ggrd_full_temp_init(struct All_variables *);
 void ggrd_reg_temp_init(struct All_variables *);
 void ggrd_temp_init_general(struct All_variables *,char *);
+void ggrd_read_mat_from_file(struct All_variables *, int );
+void xyz2rtp(float ,float ,float ,float *);
+float find_age_in_MY();
 
 /* 
 
@@ -117,7 +121,7 @@
 	/* interpolate from grid */
 	if(!ggrd_grdtrack_interpolate_tp((double)theta,(double)phi,
 					 ggrd_ict,&indbl,FALSE)){
-	  snprintf(error,255,"read_mat_from_ggrd_file: interpolation error at theta: %g phi: %g",
+	  snprintf(error,255,"ggrd_init_tracer_flavors: interpolation error at theta: %g phi: %g",
 		   theta,phi);
 	  myerror(E,error);
 	}
@@ -195,18 +199,18 @@
 		      0, MPI_COMM_WORLD, &mpi_stat);
   }
   
-  if(E->convection.ggrd_tinit_scale_with_prem){/* initialize PREM */
-    if(prem_read_model(E->convection.prem.model_filename,
-		       &E->convection.prem, (E->parallel.me==0)))
+  if(E->control.ggrd.temp_init.scale_with_prem){/* initialize PREM */
+    if(prem_read_model(E->control.ggrd.temp_init.prem.model_filename,
+		       &E->control.ggrd.temp_init.prem, (E->parallel.me == 0)))
       myerror(E,"PREM init error");
   }
   /* 
      initialize the GMT grid files 
   */
-  E->convection.ggrd_tinit_d[0].init = FALSE;
-  if(ggrd_grdtrack_init_general(TRUE,E->convection.ggrd_tinit_gfile, 
-				E->convection.ggrd_tinit_dfile,gmtflag, 
-				E->convection.ggrd_tinit_d,(E->parallel.me == 0),
+  E->control.ggrd.temp_init.d[0].init = FALSE;
+  if(ggrd_grdtrack_init_general(TRUE,E->control.ggrd.temp_init.gfile, 
+				E->control.ggrd.temp_init.dfile,gmtflag, 
+				E->control.ggrd.temp_init.d,(E->parallel.me == 0),
 				FALSE))
     myerror(E,"grd init error");
   if(E->parallel.me <  E->parallel.nproc-1){
@@ -233,7 +237,7 @@
   /* 
      mean temp is (top+bot)/2 + offset 
   */
-  tmean = (tbot + E->control.TBCtopval)/2.0 +  E->convection.ggrd_tinit_offset;
+  tmean = (tbot + E->control.TBCtopval)/2.0 +  E->control.ggrd.temp_init.offset;
 
 
   for(m=1;m <= E->sphere.caps_per_proc;m++)
@@ -246,35 +250,34 @@
 	  /* 
 	     get interpolated velocity anomaly 
 	  */
-	  if(!ggrd_grdtrack_interpolate_rtp((double)E->sx[m][3][node],
-					    (double)E->sx[m][1][node],
+	  if(!ggrd_grdtrack_interpolate_rtp((double)E->sx[m][3][node],(double)E->sx[m][1][node],
 					    (double)E->sx[m][2][node],
-					    E->convection.ggrd_tinit_d,&tadd,
+					    E->control.ggrd.temp_init.d,&tadd,
 					    FALSE))
 	    myerror(E,"ggrd_temp_init_general");
-	  if(E->convection.ggrd_tinit_scale_with_prem){
+	  if(E->control.ggrd.temp_init.scale_with_prem){
 	    /* 
 	       get the PREM density at r for additional scaling  
 	    */
-	    prem_get_rho(&rho_prem,(double)E->sx[m][3][node],&E->convection.prem);
+	    prem_get_rho(&rho_prem,(double)E->sx[m][3][node],&E->control.ggrd.temp_init.prem);
 	    if(rho_prem < 3200.0)
 	      rho_prem = 3200.0; /* we don't want the density of water */
 	    /* 
 	       assign temperature 
 	    */
-	    E->T[m][node] = tmean + tadd * E->convection.ggrd_tinit_scale * 
+	    E->T[m][node] = tmean + tadd * E->control.ggrd.temp_init.scale * 
 	      rho_prem / E->data.density;
 	  }else{
 	    /* no PREM scaling */
-	    E->T[m][node] = tmean + tadd * E->convection.ggrd_tinit_scale;
+	    E->T[m][node] = tmean + tadd * E->control.ggrd.temp_init.scale;
 	  }
 
-	  if(E->convection.ggrd_tinit_limit_trange){
+	  if(E->control.ggrd.temp_init.limit_trange){
 	    /* limit to 0 < T < 1 ?*/
 	    E->T[m][node] = min(max(E->T[m][node], 0.0),1.0);
 	  }
 	  //fprintf(stderr,"z: %11g T: %11g\n",E->sx[m][3][node],E->T[m][node]);
-	  if(E->convection.ggrd_tinit_override_tbc){
+	  if(E->control.ggrd.temp_init.override_tbc){
 	    if((k == 1) && (E->mesh.bottbc == 1)){ /* bottom TBC */
 	      E->sphere.cap[m].TB[1][node] =  E->T[m][node];
 	      E->sphere.cap[m].TB[2][node] =  E->T[m][node];
@@ -296,15 +299,195 @@
      free the structure, not needed anymore since T should now
      change internally
   */
-  ggrd_grdtrack_free_gstruc(E->convection.ggrd_tinit_d);
+  ggrd_grdtrack_free_gstruc(E->control.ggrd.temp_init.d);
   /* 
      end temperature/density from GMT grd init
   */
   temperatures_conform_bcs(E);
 }
 
+/* 
 
 
+read in material, i.e. viscosity prefactor from ggrd file, this will get assigned if 
+
+layer <=  E->control.ggrd.mat_control
+
+
+
+ */
+void ggrd_read_mat_from_file(struct All_variables *E, int is_global)
+{
+  MPI_Status mpi_stat;
+  int mpi_rc;
+  int mpi_inmsg, mpi_success_message = 1;
+  int m,el,i,j,k,inode,i1,i2,elxlz,elxlylz,ind;
+  int llayer,nox,noy,noz,nox1,noz1,noy1,lev,lselect,idim,elx,ely,elz;
+  char gmt_string[10],char_dummy;
+  double indbl,indbl2,age,f1,f2,vip;
+  float rout[3],xloc[4];
+  char tfilename[1000];
+
+  const int dims=E->mesh.nsd;
+  const int ends = enodes[dims];
+
+  
+
+
+  nox=E->mesh.nox;
+  noy=E->mesh.noy;
+  noz=E->mesh.noz;
+  nox1=E->lmesh.nox;
+  noz1=E->lmesh.noz;
+  noy1=E->lmesh.noy;
+  elx=E->lmesh.elx;
+  elz=E->lmesh.elz;
+  ely=E->lmesh.ely;
+
+  elxlz = elx * elz;
+  elxlylz = elxlz * ely;
+
+  lev=E->mesh.levmax;
+  
+  /* 
+     if we have not initialized the time history structure, do it now
+  */
+  if(!E->control.ggrd.time_hist.init){
+    /* 
+       init times, if available
+    */
+    ggrd_init_thist_from_file(&E->control.ggrd.time_hist,
+			      E->control.ggrd.time_hist.file,TRUE,(E->parallel.me == 0));
+    E->control.ggrd.time_hist.init = 1;
+  }
+  if(!E->control.ggrd.mat_control_init){
+
+    /* assign the general depth dependent material group */
+    construct_mat_group(E);
+
+    if(E->parallel.me==0)
+      fprintf(stderr,"ggrd_read_mat_from_file: initializing ggrd materials\n");
+
+    if(is_global)		/* decide on GMT flag */
+      sprintf(gmt_string,"-Lg"); /* global */
+    else
+      sprintf(gmt_string,"");
+    /* 
+       
+    initialization steps
+    
+    */
+    if(E->parallel.me > 0)	/* wait for previous processor */
+      mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), 
+			0, MPI_COMM_WORLD, &mpi_stat);
+    /* 
+       read in the material file(s)
+    */
+    E->control.ggrd.mat = (struct  ggrd_gt *)calloc(E->control.ggrd.time_hist.nvtimes,sizeof(struct ggrd_gt));
+    for(i=0;i < E->control.ggrd.time_hist.nvtimes;i++){
+      if(E->control.ggrd.time_hist.nvtimes == 1)
+	sprintf(tfilename,"%s",E->control.ggrd.mat_file);
+      else
+	sprintf(tfilename,"%s/%i/weak.grd",E->control.ggrd.mat_file,i+1);
+      /* 2D file init */
+      if(ggrd_grdtrack_init_general(FALSE,tfilename,&char_dummy,
+				    gmt_string,(E->control.ggrd.mat+i),(E->parallel.me == 0),FALSE))
+	myerror(E,"ggrd init error");
+    }
+    if(E->parallel.me <  E->parallel.nproc-1){ /* tell the next proc to go ahead */
+      mpi_rc = MPI_Send(&mpi_success_message, 1, 
+			MPI_INT, (E->parallel.me+1), 0, MPI_COMM_WORLD);
+    }else{
+      fprintf(stderr,"ggrd_read_mat_from_file: last processor done with ggrd mat init\n");
+      fprintf(stderr,"ggrd_read_mat_from_file: WARNING: assuming a regular grid geometry\n");
+    }
+    
+    /* end init */
+  }
+  if((E->control.ggrd.time_hist.nvtimes > 1)||(!E->control.ggrd.mat_control_init)){
+    /* 
+       loop through all elements and assign
+    */
+    for (m=1;m <= E->sphere.caps_per_proc;m++) {
+      for (j=1;j <= elz;j++)  {	/* this assumes a regular grid sorted as in (1)!!! */
+	if(E->mat[m][j] <= E->control.ggrd.mat_control ){
+	  /* 
+	     lithosphere or asthenosphere 
+	  */
+	  for (k=1;k <= ely;k++){
+	    for (i=1;i <= elx;i++)   {
+	      /* eq.(1) */
+	      el = j + (i-1) * elz + (k-1)*elxlz;
+	      /* 
+		 find average horizontal coordinate 
+		 
+		 (DO WE HAVE THIS STORED ALREADY, E.G. FROM PRESSURE
+		 EVAL FORM FUNCTION???)
+	      */
+	      xloc[1] = xloc[2] = xloc[3] = 0.0;
+	      for(inode=1;inode <= 4;inode++){
+		ind = E->ien[m][el].node[inode];
+		xloc[1] += E->x[m][1][ind];xloc[2] += E->x[m][2][ind];xloc[3] += E->x[m][3][ind];
+	      }
+	      xloc[1]/=4.;xloc[2]/=4.;xloc[3]/=4.;
+	      xyz2rtp(xloc[1],xloc[2],xloc[3],rout);
+	      /* interpolate material */
+	      if(E->control.ggrd.time_hist.nvtimes == 1){
+		if(!ggrd_grdtrack_interpolate_tp((double)rout[1],(double)rout[2],E->control.ggrd.mat,
+						 &vip,FALSE)){
+		  fprintf(stderr,"ggrd_read_mat_from_ggrd_file: interpolation error at %g, %g\n",
+			  xloc[1],xloc[2]);
+		  parallel_process_termination();
+		}
+	      }else{		
+		/* 
+		   interpolate by time 
+		*/
+		age = find_age_in_MY(E);
+		/*  */
+		ggrd_interpol_time(age,&E->control.ggrd.time_hist,&i1,&i2,&f1,&f2,
+				   E->control.ggrd.time_hist.vstage_transition);
+		/*  */
+		if(!ggrd_grdtrack_interpolate_tp((double)rout[1],(double)rout[2],(E->control.ggrd.mat+i1),&indbl,FALSE)){
+		  fprintf(stderr,"ggrd_read_mat_from_ggrd_file: interpolation error at %g, %g\n",
+			  xloc[1],xloc[2]);
+		  parallel_process_termination();
+		}
+		if(!ggrd_grdtrack_interpolate_tp((double)rout[1],(double)rout[2],(E->control.ggrd.mat+i2),&indbl2,FALSE)){
+		  fprintf(stderr,"ggrd_read_mat_from_ggrd_file: interpolation error at %g, %g\n",
+			  xloc[1],xloc[2]);
+		  parallel_process_termination();
+		}
+		/* average smoothly between the two tectonic stages */
+
+		vip = exp((f1*log(indbl)+f2*log(indbl2)));
+		//fprintf(stderr,"%g %i %i %g %g %g %g -> %g\n",age, i1,i2,f1,f2,indbl,indbl2,vip);
+	      }
+	      /* limit the input scaling? */
+	      if(vip < 1e-5)
+	      	vip = 1e-5;
+	      if(vip > 1e5)
+	      	vip = 1e5;
+	      E->VIP[m][el] = vip;
+	    }
+	  }
+	}else{
+	  /* outside the lithosphere */
+	  for (k=1;k <= ely;k++){
+	    for (i=1;i <= elx;i++)   {
+	      el = j + (i-1) * elz + (k-1)*elxlz;
+	      /* no scaling else */
+	      E->VIP[m][el] = 1.0;
+	    }
+	  }
+	}
+      }	/* end elz loop */
+    } /* end m loop */
+  } /* end assignment loop */
+
+  E->control.ggrd.mat_control_init = 1;
+}
+
 #endif
 
 

Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c	2008-02-17 00:11:30 UTC (rev 11176)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c	2008-02-17 00:14:06 UTC (rev 11177)
@@ -150,19 +150,19 @@
 #ifdef USE_GGRD
     /* read in some more parameters */
     /* scale the anomalies with PREM densities */
-    input_boolean("ggrd_tinit_scale_with_prem",&(E->convection.ggrd_tinit_scale_with_prem),"off",E->parallel.me);
+    input_boolean("ggrd_tinit_scale_with_prem",&(E->control.ggrd.temp_init.scale_with_prem),"off",E->parallel.me);
     /* limit T to 0...1 */
-    input_boolean("ggrd_tinit_limit_trange",&(E->convection.ggrd_tinit_limit_trange),"on",E->parallel.me);
+    input_boolean("ggrd_tinit_limit_trange",&(E->control.ggrd.temp_init.limit_trange),"on",E->parallel.me);
    /* scaling factor for the grids */
-    input_double("ggrd_tinit_scale",&(E->convection.ggrd_tinit_scale),"1.0",E->parallel.me); /* scale */
+    input_double("ggrd_tinit_scale",&(E->control.ggrd.temp_init.scale),"1.0",E->parallel.me); /* scale */
     /* temperature offset factor */
-    input_double("ggrd_tinit_offset",&(E->convection.ggrd_tinit_offset),"0.0",E->parallel.me); /* offset */
+    input_double("ggrd_tinit_offset",&(E->control.ggrd.temp_init.offset),"0.0",E->parallel.me); /* offset */
     /* grid name, without the .i.grd suffix */
-    input_string("ggrd_tinit_gfile",E->convection.ggrd_tinit_gfile,"",E->parallel.me); /* grids */
-    input_string("ggrd_tinit_dfile",E->convection.ggrd_tinit_dfile,"",E->parallel.me); /* depth.dat layers of grids*/
+    input_string("ggrd_tinit_gfile",E->control.ggrd.temp_init.gfile,"",E->parallel.me); /* grids */
+    input_string("ggrd_tinit_dfile",E->control.ggrd.temp_init.dfile,"",E->parallel.me); /* depth.dat layers of grids*/
     /* override temperature boundary condition? */
-    input_boolean("ggrd_tinit_override_tbc",&(E->convection.ggrd_tinit_override_tbc),"off",E->parallel.me);
-    input_string("ggrd_tinit_prem_file",E->convection.prem.model_filename,"", E->parallel.me); /* PREM model filename */
+    input_boolean("ggrd_tinit_override_tbc",&(E->control.ggrd.temp_init.override_tbc),"off",E->parallel.me);
+    input_string("ggrd_tinit_prem_file",E->control.ggrd.temp_init.prem.model_filename,"", E->parallel.me); /* PREM model filename */
 #else
     fprintf(stderr,"tic_method 4 only works for USE_GGRD compiled code\n");
     parallel_process_termination();

Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c	2008-02-17 00:11:30 UTC (rev 11176)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c	2008-02-17 00:14:06 UTC (rev 11177)
@@ -443,6 +443,40 @@
   input_int("mat_control",&(E->control.mat_control),"0",m);
   input_string("mat_file",E->control.mat_file,"",m);
 
+#ifdef USE_GGRD
+  /* 
+     usage:
+     (a) 
+
+     ggrd_mat_control=2
+     ggrd_mat_file="weak.grd"
+
+     read in time-constant prefactors from weak.grd netcdf file that apply to top two E->mat layers
+
+     (b)
+
+     ggrd_mat_control=2
+     ggrd_mat_file="mythist"
+     ggrd_time_hist_file="mythist/times.dat"
+
+
+     time-dependent, will look for n files named
+     mythist/i/weak.grd where i = 1...n and n is the number of times as specified in ggrd_time_hist_file
+     which has time in Ma for n stages like so
+     
+     -60 -30
+     -30 -15
+     -15 0
+     
+  */
+  ggrd_init_master(&E->control.ggrd);
+  input_string("ggrd_time_hist_file",E->control.ggrd.time_hist.file,"",m); /* time history file, if not specified, will use constant VBCs and material grids */
+  input_int("ggrd_mat_control",&(E->control.ggrd.mat_control),"0",m); /* if > 0, will use top  E->control.ggrd.mat_control layers and assign a prefactor for the viscosity */
+  input_string("ggrd_mat_file",E->control.ggrd.mat_file,"",m); /* file to read prefactors from */
+  if(E->control.ggrd.mat_control) /* this will override mat_control setting */
+    E->control.mat_control = 1;
+#endif
+
   input_int("nodex",&(E->mesh.nox),"essential",m);
   input_int("nodez",&(E->mesh.noz),"essential",m);
   input_int("nodey",&(E->mesh.noy),"essential",m);

Modified: mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c	2008-02-17 00:11:30 UTC (rev 11176)
+++ mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c	2008-02-17 00:14:06 UTC (rev 11177)
@@ -142,7 +142,9 @@
         break;
 
       case 3:  /* read element materials */
-
+#ifdef USE_GGRD
+	if(!E->control.ggrd.mat_control){
+#endif
         sprintf(output_file1,"%s%0.0f.0",E->control.mat_file,newage1);
         sprintf(output_file2,"%s%0.0f.0",E->control.mat_file,newage2);
         fp1=fopen(output_file1,"r");
@@ -166,6 +168,10 @@
             fprintf(E->fp,"Mat: File2 = No file inputted (negative age)\n");
         }
 
+#ifdef USE_GGRD
+	}
+#endif
+	break;
     } /* end switch */
 
 
@@ -232,6 +238,11 @@
         break;
 
       case 3:  /* read element materials */
+#ifdef USE_GGRD
+	if(E->control.ggrd.mat_control){
+	  ggrd_read_mat_from_file(E, 0);
+	}else{
+#endif
 
         VIP1 = (float*) malloc ((emax+1)*sizeof(float));
         VIP2 = (float*) malloc ((emax+1)*sizeof(float));
@@ -271,7 +282,10 @@
          free ((void *) VIP2);
          free ((void *) LL1);
          free ((void *) LL2);
-
+#ifdef USE_GGRD
+	} /* end of branch if allowing for ggrd handling */
+#endif
+	break;
     } /* end switch */
 
    return;

Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2008-02-17 00:11:30 UTC (rev 11176)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2008-02-17 00:14:06 UTC (rev 11177)
@@ -444,12 +444,12 @@
                         zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
                     }
 
-                    if(E->control.mat_control==0)
+                    if(E->control.mat_control==0){
                         EEta[m][ (i-1)*vpts + jj ] = tempa*
                             exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
                                  - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
 
-                    if(E->control.mat_control==1) {
+                    }else{
                      /* visc1 = E->VIP[m][i];
                         visc2 = 2.0/(1./visc1 + 1.);
                         tempa_exp = tempa*
@@ -501,7 +501,7 @@
 		exp( E->viscosity.E[l]*(E->viscosity.T[l] - temp) +
 		     zzz *  E->viscosity.Z[l]);
 	      //fprintf(stderr,"N0 %11g T %11g T0 %11g E %11g z %11g km Z %11g mat: %i log10(eta): %11g\n",
-	      //      tempa,temp,E->viscosity.T[l],E->viscosity.E[l], zzz *6371 ,E->viscosity.Z[l],l+1,log10(EEta[m][ (i-1)*vpts + jj ]));
+	      //tempa,temp,E->viscosity.T[l],E->viscosity.E[l], zzz *6371 ,E->viscosity.Z[l],l+1,log10(EEta[m][ (i-1)*vpts + jj ]));
 	    }
 	  }
         break;

Modified: mc/3D/CitcomS/trunk/lib/convection_variables.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/convection_variables.h	2008-02-17 00:11:30 UTC (rev 11176)
+++ mc/3D/CitcomS/trunk/lib/convection_variables.h	2008-02-17 00:14:06 UTC (rev 11177)
@@ -51,18 +51,8 @@
 	}  heat_sources;
 
 
-#ifdef USE_GGRD
-  /* for temperature init from grd files */
-  int ggrd_tinit,ggrd_tinit_scale_with_prem;
-  int ggrd_tinit_override_tbc,ggrd_tinit_limit_trange;
-  double ggrd_tinit_scale,ggrd_tinit_offset,ggrd_vstage_transition;
-  char ggrd_tinit_gfile[1000];
-  char ggrd_tinit_dfile[1000];
-  struct ggrd_gt ggrd_tinit_d[1];
-  struct ggrd_t ggrd_time_hist;
-  struct prem_model prem; 
-#endif
 
+
 } convection;
 
 

Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h	2008-02-17 00:11:30 UTC (rev 11176)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h	2008-02-17 00:14:06 UTC (rev 11177)
@@ -511,6 +511,9 @@
     int side_sbcs;
     int vbcs_file;
     int mat_control;
+#ifdef USE_GGRD
+  struct ggrd_master ggrd;
+#endif
     double accuracy;
     double vaccuracy;
     char velocity_boundary_file[1000];



More information about the cig-commits mailing list