Hopefully fixed merge nightmare.

bin/Citcom.c

@@ -272,7 +272,9 @@ int main(argc,argv)
 
 
     if (E->parallel.me == 0)  {
-      fprintf(E->fp,"CPU total = %g & CPU = %g for step %d time = %.4e dt = %.4e  maxT = %.4e sub_iteration%d\n",CPU_time0()-start_time,CPU_time0()-time,E->monitor.solution_cycles,E->monitor.elapsed_time,E->advection.timestep,E->monitor.T_interior,E->advection.last_sub_iterations);
+      fprintf(E->fp,"CPU total = %g & CPU = %g for step %d time = %.4e dt = %.4e  maxT = %.4e sub_iteration%d\n",
+	      CPU_time0()-start_time,CPU_time0()-time,E->monitor.solution_cycles,E->monitor.elapsed_time,
+	      E->advection.timestep,E->monitor.T_interior,E->advection.last_sub_iterations);
 
       time = CPU_time0();
     }

lib/Advection_diffusion.c

@@ -64,7 +64,7 @@ static void element_residual(struct All_variables *E, int el,
 static void filter(struct All_variables *E);
 static void process_heating(struct All_variables *E, int psc_pass);
 
-/* ============================================
+/* =========================================
    Generic adv-diffusion for temperature field.
    ============================================ */
 
@@ -302,13 +302,19 @@ void PG_timestep_solve(struct All_variables *E)
       }
   }
 
+
+#ifdef USE_GGRD
+  if((E->control.lith_age) && (!E->control.ggrd.age_control)){
+#else
   if(E->control.lith_age) {
-      if(E->parallel.me==0) fprintf(stderr,"PG_timestep_solve - lith_age adjust\n");
-      lith_age_conform_tbc(E);
-      assimilate_lith_conform_bcs(E);
+#endif
+    /* this is the old mode, it didn't quite work for me, so I'll not
+       use it with ggrd */
+    if(E->parallel.me==0) fprintf(stderr,"PG_timestep_solve - lith_age adjust\n");
+    lith_age_conform_tbc(E);
+    assimilate_lith_conform_bcs(E);
+
   }
-
-
   return;
 }
 

lib/BC_util.c

@@ -32,6 +32,7 @@ void internal_horizontal_bc(struct All_variables *,float *[],int,int,float,unsig
 void myerror(struct All_variables *,char *);
 int layers(struct All_variables *,int,int);
 
+#include "lith_age.h"
 
 #ifdef USE_GGRD
 #include "ggrd_handling.h"
@@ -105,18 +106,26 @@ void strip_bcs_from_residual(E,Res,level)
 void temperatures_conform_bcs(E)
      struct All_variables *E;
 {
-  void temperatures_conform_bcs2(struct All_variables *);
-  void assimilate_lith_conform_bcs2(struct All_variables *);
-
   if(E->control.lith_age) {
+#ifdef USE_GGRD
+    if(E->control.ggrd.age_control){ 
+      if(!E->control.ggrd.vtop_control_init) 
+	myerror(E,"temperature_conform_bcs: error, ggrd age control was not initialized");
+      set_lith_age_for_t_and_tbc(E,TRUE); /* reassign for T and TBC */
+      temperatures_conform_bcs2(E);	  
+    }
+#endif    
     /*
     This sequence now moved to end of PG_time_step_solve
     lith_age_conform_tbc(E);
     assimilate_lith_conform_bcs(E);
+
+    this didn't work for me
+
     */
-    }
-  else
+  }else{
     temperatures_conform_bcs2(E);
+  }
   return;
 }
 

lib/Full_boundary_conditions.c

@@ -89,7 +89,7 @@ void full_velocity_boundary_conditions(E)
         horizontal_bc(E,E->sphere.cap[j].VB,noz,2,0.0,SBY,0,lv,j);
 
 #ifdef USE_GGRD
-	/* Ggrd velocity control */
+	/* Ggrd velocity (and age) control (only for the top processors */
 	if((lv==E->mesh.gridmax) && E->control.ggrd.vtop_control)
 	  ggrd_read_vtop_from_file(E,TRUE);
 #endif
@@ -119,24 +119,25 @@ void full_velocity_boundary_conditions(E)
       if(E->control.side_sbcs)
 	apply_side_sbc(E);
 
-/* if(E->control.verbose) { */
-/*  for (j=1;j<=E->sphere.caps_per_proc;j++) */
-/*    for (node=1;node<=E->lmesh.nno;node++) */
-/*       fprintf(E->fp_out,"m=%d VB== %d %g %g %g flag %u %u %u\n",j,node,E->sphere.cap[j].VB[1][node],E->sphere.cap[j].VB[2][node],E->sphere.cap[j].VB[3][node],E->node[j][node]&VBX,E->node[j][node]&VBY,E->node[j][node]&VBZ); */
-/*  fflush(E->fp_out); */
-/* } */
-
-  /* If any imposed internal velocity structure it goes here */
+      /* if(E->control.verbose) { */
+      /*  for (j=1;j<=E->sphere.caps_per_proc;j++) */
+      /*    for (node=1;node<=E->lmesh.nno;node++) */
+      /*       fprintf(E->fp_out,"m=%d VB== %d %g %g %g flag %u %u %u\n",j,node,E->sphere.cap[j].VB[1][node],E->sphere.cap[j].VB[2][node],E->sphere.cap[j].VB[3][node],E->node[j][node]&VBX,E->node[j][node]&VBY,E->node[j][node]&VBZ); */
+      /*  fflush(E->fp_out); */
+      /* } */
 
+      /* If any imposed internal velocity structure it goes here */
+
 
       /*
 	apply stress or velocity boundary conditions, read from file
 	settings are to be implemented in those routines (will only do
 	anything at present, if E->mesh.toplayerbc != 0
       */
       assign_internal_bc(E);
-
-   return; }
+
+      return;
+}
 
 /* ========================================== */
 
@@ -175,24 +176,19 @@ void full_temperature_boundary_conditions(E)
       horizontal_bc(E,E->sphere.cap[j].TB,1,3,E->control.TBCbotval,TBZ,0,lev,j);
       horizontal_bc(E,E->sphere.cap[j].TB,1,3,E->control.TBCbotval,FBZ,1,lev,j);
     }
-
     if(E->control.lith_age_time==1)  {
-
       /* set the regions in which to use lithosphere files to determine temperature
 	 note that this is called if the lithosphere age in inputted every time step
 	 OR it is only maintained in the boundary regions */
       lith_age_temperature_bound_adj(E,lev);
     }
-
-
   }     /* end for j */
 
   temperatures_conform_bcs(E);
   E->temperatures_conform_bcs = temperatures_conform_bcs;
 
-   return; }
-
-
+  return; 
+}
 /*  =========================================================  */
 
 void horizontal_bc(struct All_variables *E,float *BC[],int ROW,int dirn,float value,
@@ -237,11 +233,7 @@ void horizontal_bc(struct All_variables *E,float *BC[],int ROW,int dirn,float va
   return;
 }
 
-/* 
-
-   for latitude dependent surface temperature 
-
- */
+/* latitudinally dependent horizontal BC */
 
 void horizontal_bc_lat_dep(struct All_variables *E,float *BC[],int ROW,int dirn,
 			   unsigned int mask,char onoff,int level,int m)
@@ -284,6 +276,7 @@ void horizontal_bc_lat_dep(struct All_variables *E,float *BC[],int ROW,int dirn,
 
   return;
 }
+
 /* temp as f(latitude) */
 float lat_dep_temp(struct All_variables *E, int node)
 {
@@ -307,7 +300,7 @@ void temperature_apply_periodic_bcs(E)
  fprintf(E->fp,"Periodic temperature boundary conditions\n");
 
   return;
-  }
+}
 
 
 

lib/Ggrd_handling.c

@@ -282,7 +282,7 @@ void ggrd_temp_init_general(struct All_variables *E,int is_geographic)
 	for(k=1;k <= noz;k++)  {
 	  /* node numbers */
 	  node=k+(j-1)*noz+(i-1)*noxnoz;
-
+	 
 	  /*
 	     get interpolated velocity anomaly
 	  */
@@ -367,14 +367,10 @@ void ggrd_temp_init_general(struct All_variables *E,int is_geographic)
   ggrd_grdtrack_free_gstruc(E->control.ggrd.temp.d);
 
   if(E->control.ggrd.age_control){
-    /* age TBC control overrides this */
-
-    /* set the TBC boundary condition and set T as well within TBC*/
-    set_lith_age_for_t_and_tbc(E);
-    /* set the boundary flags */
+    /* set boundary flags for the nodels within the age defined
+       thermal boundary layer */
     lith_age_temperature_bound_adj(E,E->mesh.gridmax);
   }
-
   /*
      end temperature/density from GMT grd init
   */
@@ -406,7 +402,7 @@ void ggrd_read_mat_from_file(struct All_variables *E, int is_geographic)
   int mpi_rc,timedep,interpolate;
   int mpi_inmsg, mpi_success_message = 1;
   int m,el,i,j,k,inode,i1,i2,elxlz,elxlylz,ind;
-  int llayer,nox,noy,noz,level,lselect,idim,elx,ely,elz;
+  int llayer,level,lselect,idim,elx,ely,elz;
   char gmt_string[10],char_dummy;
   double indbl,indbl2,age,f1,f2,vip,rout[3],xloc[4];
   char tfilename[1000];
@@ -416,7 +412,6 @@ void ggrd_read_mat_from_file(struct All_variables *E, int is_geographic)
   FILE *in;
   ggrd_boolean use_nearneighbor;
 
-  nox=E->mesh.nox;noy=E->mesh.noy;noz=E->mesh.noz;
   elx=E->lmesh.elx;elz=E->lmesh.elz;ely=E->lmesh.ely;
   elxlz = elx * elz;
   elxlylz = elxlz * ely;
@@ -511,8 +506,7 @@ void ggrd_read_mat_from_file(struct All_variables *E, int is_geographic)
     if(E->parallel.me == 0)
       fprintf(stderr,"ggrd_read_mat_from_file: assigning at age %g\n",age);
     if(timedep){
-      if(!ggrd_interpol_time(age,&E->control.ggrd.time_hist,&i1,&i2,&f1,&f2,
-			     E->control.ggrd.time_hist.vstage_transition))
+      if(!ggrd_interpol_time(age,&E->control.ggrd.time_hist,&i1,&i2,&f1,&f2))
 	myerror(E,"interpolation error");
       interpolate = 1;
     }else{
@@ -681,7 +675,7 @@ void ggrd_read_ray_from_file(struct All_variables *E, int is_geographic)
   int mpi_rc,timedep,interpolate;
   int mpi_inmsg, mpi_success_message = 1;
   int m,el,i,j,k,node,i1,i2,elxlz,elxlylz,ind;
-  int llayer,nox,noy,noz,lev,lselect,idim,elx,ely,elz;
+  int llayer,lev,lselect,idim,elx,ely,elz;
   char gmt_string[10],char_dummy;
   double indbl,indbl2,age,f1,f2,vip,rout[3],xloc[4];
   char tfilename[1000];
@@ -691,7 +685,6 @@ void ggrd_read_ray_from_file(struct All_variables *E, int is_geographic)
   const int dims=E->mesh.nsd;
   const int ends = enodes[dims];
   /* dimensional ints */
-  nox=E->mesh.nox;noy=E->mesh.noy;noz=E->mesh.noz;
   elx=E->lmesh.elx;elz=E->lmesh.elz;ely=E->lmesh.ely;
   elxlz = elx * elz;
   elxlylz = elxlz * ely;
@@ -745,8 +738,7 @@ void ggrd_read_ray_from_file(struct All_variables *E, int is_geographic)
   if(timedep || (!E->control.ggrd.ray_control_init)){
     if(timedep){
       age = find_age_in_MY(E);
-      if(!ggrd_interpol_time(age,&E->control.ggrd.time_hist,&i1,&i2,&f1,&f2,
-			     E->control.ggrd.time_hist.vstage_transition))
+      if(!ggrd_interpol_time(age,&E->control.ggrd.time_hist,&i1,&i2,&f1,&f2))
 	myerror(E,"interpolation error");
       interpolate = 1;
     }else{
@@ -945,8 +937,12 @@ void ggrd_read_vtop_from_file(struct All_variables *E, int is_geographic)
       ggrd_init_thist_from_file(&E->control.ggrd.time_hist,E->control.ggrd.time_hist.file,
 				TRUE,(E->parallel.me == 0));
       if(E->control.ggrd_smooth_stages){
+
+	if(E->control.ggrd.time_hist.interpol_time_lin)
+	  myerror(E,"smoothing the stages and linear interpolation does not make sense");
+
 	/* smooth transition between stages, else this is 0.1 by default */
-	E->control.ggrd.time_hist.vstage_transition = 0;
+	E->control.ggrd.time_hist.vstage_transition = 0.0;
 	for(i=1;i<E->control.ggrd.time_hist.nvtimes;i++)
 	  E->control.ggrd.time_hist.vstage_transition += 
 	    (E->control.ggrd.time_hist.vtimes[i*3+1]-E->control.ggrd.time_hist.vtimes[(i-1)*3+1]);
@@ -956,8 +952,12 @@ void ggrd_read_vtop_from_file(struct All_variables *E, int is_geographic)
     }
     timedep = (E->control.ggrd.time_hist.nvtimes > 1)?(1):(0);
     if(verbose && !E->control.ggrd.vtop_control_init ){
-      fprintf(stderr,"ggrd_read_vtop_from_file: time dependence: %i, nsteps: %i stage_tran: %g\n",
-	      timedep,E->control.ggrd.time_hist.nvtimes,E->control.ggrd.time_hist.vstage_transition);
+      if(E->control.ggrd.time_hist.interpol_time_lin)
+	fprintf(stderr,"ggrd_read_vtop_from_file: time dependence: %i, nsteps: %i linear interpolation\n",
+		timedep,E->control.ggrd.time_hist.nvtimes);
+      else
+	fprintf(stderr,"ggrd_read_vtop_from_file: time dependence: %i, nsteps: %i stage_tran: %g\n",
+		timedep,E->control.ggrd.time_hist.nvtimes,E->control.ggrd.time_hist.vstage_transition);
 
     }
     if(!E->control.ggrd.vtop_control_init){
@@ -1132,13 +1132,18 @@ void ggrd_read_vtop_from_file(struct All_variables *E, int is_geographic)
 	    fprintf(stderr,"ggrd_read_vtop_from_file: using present day vtop for age = %g\n",age);
 	}else{
 	  /*  */
-	  if(!ggrd_interpol_time(age,&E->control.ggrd.time_hist,&i1,&i2,&f1,&f2,
-				 E->control.ggrd.time_hist.vstage_transition))
+	  if(!ggrd_interpol_time(age,&E->control.ggrd.time_hist,&i1,&i2,&f1,&f2))
 	    myerror(E,"interpolation error");
 	  interpolate = 1;
-	  if(verbose)
-	    fprintf(stderr,"ggrd_read_vtop_from_file: interpolating for age  %g (%.3f of %g Ma, %.3f of %g Ma)\n",
-		    age,f1,E->control.ggrd.time_hist.vtimes[i1*3+1],f2,E->control.ggrd.time_hist.vtimes[i2*3+1]);
+	  if(verbose){
+	    if(E->control.ggrd.time_hist.interpol_time_lin){
+	      fprintf(stderr,"ggrd_read_vtop_from_file: lin interp for age  %g (%.3f of %g Ma, %.3f of %g Ma)\n",
+		      age,f1,E->control.ggrd.time_hist.tl[i1],f2,E->control.ggrd.time_hist.tl[i2]);
+	    }else{
+	      fprintf(stderr,"ggrd_read_vtop_from_file: stage interp for age  %g (%.3f of %g Ma, %.3f of %g Ma)\n",
+		      age,f1,E->control.ggrd.time_hist.vtimes[i1*3+1],f2,E->control.ggrd.time_hist.vtimes[i2*3+1]);
+	    }
+	  }
 	}
 
       }else{
@@ -1149,8 +1154,10 @@ void ggrd_read_vtop_from_file(struct All_variables *E, int is_geographic)
       }
 
       if(verbose)
-	fprintf(stderr,"ggrd_read_vtop_from_file: assigning %s BC, timedep: %i time: %g Ma ndtime: %g\n",
-		(use_vel)?("velocities"):("tractions"),	timedep,age,E->monitor.elapsed_time);
+	fprintf(stderr,"ggrd_read_vtop_from_file: assigning %s BC%s, tdep: %i time: %g Ma ndtime: %g\n",
+		(use_vel)?("velocities"):("tractions"),	
+		(E->control.ggrd.age_control)?(" and ages"):(""),
+		timedep,age,E->monitor.elapsed_time);
 
       /* if mixed BCs are allowed, need to reassign the boundary
 	 condition */
@@ -1505,7 +1512,11 @@ void ggrd_read_vtop_from_file(struct All_variables *E, int is_geographic)
 
 
   } /* end top processor or allow internal branch branch */
-
+  if(! E->control.ggrd.vtop_control_init && verbose)
+    fprintf(stderr,"vtop_init from ggrd done\n");
+  parallel_process_sync(E);
+  if(E->control.ggrd.age_control)
+    E->control.ggrd.age_control_init = TRUE;
   E->control.ggrd.vtop_control_init = TRUE;
 
 }
@@ -1566,10 +1577,10 @@ void ggrd_adjust_tbl_rayleigh(struct All_variables *E,
   double xloc,fac,bnew;
   if(!E->control.ggrd.ray_control_init)
     myerror(E,"ggrd rayleigh not initialized, but in adjust tbl");
-  if(E->parallel.me == 0)
+  if(E->parallel.me == 0){
     fprintf(stderr,"ggrd__adjust_tbl_rayleigh: adjusting Rayleigh in top %i layers\n",
 	    E->control.ggrd.ray_control);
-
+  }
   /* 
      need to scale buoy with the material determined rayleigh numbers
   */

lib/Initial_temperature.c

@@ -658,8 +658,7 @@ static void construct_tic_from_input(struct All_variables *E)
 #ifdef USE_GGRD
         ggrd_temp_init_general(E,1);
 #else
-        fprintf(stderr,"tic_method 4 only works for USE_GGRD compiled code\n");
-        parallel_process_termination();
+        myerror(E,"tic_method 4 only works for USE_GGRD compiled code\n");
 #endif
         break;