[cig-commits] r11757 - in mc/3D/CitcomS/trunk: bin lib
becker at geodynamics.org
becker at geodynamics.org
Sun Apr 6 14:58:45 PDT 2008
Author: becker
Date: 2008-04-06 14:58:45 -0700 (Sun, 06 Apr 2008)
New Revision: 11757
Modified:
mc/3D/CitcomS/trunk/bin/Citcom.c
mc/3D/CitcomS/trunk/lib/Advection_diffusion.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/Instructions.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Problem_related.c
mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
mc/3D/CitcomS/trunk/lib/ggrd_handling.h
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
Log:
GMT/Netcdf grd input can now deal with velocity boundary conditions,
material dependence, and local Rayleigh number in surface layers.
Modified: mc/3D/CitcomS/trunk/bin/Citcom.c
===================================================================
--- mc/3D/CitcomS/trunk/bin/Citcom.c 2008-04-06 20:27:32 UTC (rev 11756)
+++ mc/3D/CitcomS/trunk/bin/Citcom.c 2008-04-06 21:58:45 UTC (rev 11757)
@@ -58,7 +58,11 @@
void vcopy();
void construct_mat_group();
void read_velocity_boundary_from_file();
+ void read_rayleigh_from_file();
void read_mat_from_file();
+#ifdef USE_GGRD
+ void read_rayleigh_from_file();
+#endif
void open_time();
void output_finalize();
void PG_timestep_init();
@@ -209,14 +213,15 @@
if ((E->monitor.solution_cycles % E->control.checkpoint_frequency)==0) {
output_checkpoint(E);
}
-
/* updating time-dependent material group
* if mat_control is 0, the material group has already been
* initialized in initial_conditions() */
if(E->control.mat_control==1)
read_mat_from_file(E);
-
-
+#ifdef USE_GGRD
+ if(E->control.ggrd.ray_control)
+ read_rayleigh_from_file(E);
+#endif
if(E->control.vbcs_file==1)
read_velocity_boundary_from_file(E);
/*
Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2008-04-06 20:27:32 UTC (rev 11756)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2008-04-06 21:58:45 UTC (rev 11757)
@@ -30,10 +30,10 @@
Brooks PhD thesis (Caltech) which refers back to Hughes, Liu and Brooks. */
#include <sys/types.h>
-#include <math.h>
+
#include "element_definitions.h"
#include "global_defs.h"
-
+#include <math.h>
#include "advection_diffusion.h"
#include "parsing.h"
Modified: mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c 2008-04-06 20:27:32 UTC (rev 11756)
+++ mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c 2008-04-06 21:58:45 UTC (rev 11757)
@@ -188,7 +188,8 @@
}
#endif
break;
-
+ /* mode 4 is rayleigh control for GGRD, see below */
+
} /* end switch */
@@ -272,8 +273,7 @@
#endif
break;
- case 3: /* read element materials */
-
+ case 3: /* read element materials and Ray */
#ifdef USE_GGRD
if(E->control.ggrd.mat_control){ /* use netcdf grids */
ggrd_read_mat_from_file(E, 1);
@@ -325,6 +325,14 @@
} /* end of branch if allowing for ggrd handling */
#endif
break;
+ case 4: /* material control */
+#ifdef USE_GGRD
+ if(E->control.ggrd.ray_control)
+ ggrd_read_ray_from_file(E, 1);
+#else
+ myerror(E,"input_from_files: mode 4 only for GGRD");
+#endif
+ break;
} /* end switch */
} /* end for m */
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2008-04-06 20:27:32 UTC (rev 11756)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2008-04-06 21:58:45 UTC (rev 11757)
@@ -385,8 +385,10 @@
/* end init */
}
if(timedep || (!E->control.ggrd.mat_control_init)){
+ age = find_age_in_MY(E);
+ if(E->parallel.me == 0)
+ fprintf(stderr,"ggrd_read_mat_from_ggrd_file: assigning at age %g\n",age);
if(timedep){
- 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);
interpolate = 1;
@@ -423,14 +425,14 @@
/* material */
if(!ggrd_grdtrack_interpolate_tp(rout[1],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]);
+ rout[1],rout[2]);
parallel_process_termination();
}
if(interpolate){
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]);
+ rout[1],rout[2]);
parallel_process_termination();
}
/* average smoothly between the two tectonic stages */
@@ -459,12 +461,133 @@
} /* end elz loop */
} /* end m loop */
} /* end assignment loop */
-
+ if((!timedep) && (!E->control.ggrd.mat_control_init)){ /* forget the grid */
+ ggrd_grdtrack_free_gstruc(E->control.ggrd.mat);
+ }
E->control.ggrd.mat_control_init = 1;
} /* end mat control */
+/*
+
+read in Rayleigh number prefactor from file, this will get assigned if
+
+layer <= E->control.ggrd.ray_control
+
+
+*/
+void ggrd_read_ray_from_file(struct All_variables *E, int is_global)
+{
+ MPI_Status mpi_stat;
+ 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,nox1,noz1,noy1,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];
+ 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;
+ 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
+ any function can do that
+
+ we could only use the surface processors, but maybe the rayleigh
+ number is supposed to be changed at large depths
+ */
+ if(!E->control.ggrd.time_hist.init){
+ 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;
+ }
+ timedep = (E->control.ggrd.time_hist.nvtimes > 1)?(1):(0);
+ if(!E->control.ggrd.ray_control_init){
+ /* init step */
+ if(E->parallel.me==0)
+ fprintf(stderr,"ggrd_read_ray_from_file: initializing from %s\n",E->control.ggrd.ray_file);
+ if(is_global) /* decide on GMT flag */
+ sprintf(gmt_string,GGRD_GMT_GLOBAL_STRING); /* global */
+ else
+ sprintf(gmt_string,"");
+ 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);
+ E->control.ggrd.ray = (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(!timedep) /* constant */
+ sprintf(tfilename,"%s",E->control.ggrd.ray_file);
+ else
+ sprintf(tfilename,"%s/%i/rayleigh.grd",E->control.ggrd.ray_file,i+1);
+ if(ggrd_grdtrack_init_general(FALSE,tfilename,&char_dummy,
+ gmt_string,(E->control.ggrd.ray+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_ray_from_file: last processor done with ggrd ray init\n");
+ }
+ E->control.surface_rayleigh = (float *)malloc(sizeof(float)*(E->lmesh.nsf+2));
+ if(!E->control.surface_rayleigh)
+ myerror(E,"ggrd rayleigh mem error");
+ }
+ if(timedep || (!E->control.ggrd.ray_control_init)){
+ if(timedep){
+ 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);
+ interpolate = 1;
+ }else{
+ interpolate = 0;i1 = 0;
+ }
+ if(E->parallel.me == 0)
+ fprintf(stderr,"ggrd_read_ray_from_ggrd_file: assigning at time %g\n",age);
+ for (m=1;m <= E->sphere.caps_per_proc;m++) {
+ /* loop through all surface nodes */
+ for (j=1;j <= E->lmesh.nsf;j++) {
+ node = j * E->lmesh.noz ;
+ rout[1] = (double)E->sx[m][1][node];
+ rout[2] = (double)E->sx[m][2][node];
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.ray+i1),&indbl,FALSE)){
+ fprintf(stderr,"ggrd_read_ray_from_ggrd_file: interpolation error at %g, %g\n",
+ rout[1],rout[2]);
+ parallel_process_termination();
+ }
+ //fprintf(stderr,"%i %i %g %g %g\n",j,E->lmesh.nsf,rout[1],rout[2],indbl);
+ if(interpolate){
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],
+ (E->control.ggrd.ray+i2),&indbl2,FALSE)){
+ fprintf(stderr,"ggrd_read_ray_from_ggrd_file: interpolation error at %g, %g\n",
+ rout[1],rout[2]);
+ parallel_process_termination();
+ }
+ /* average smoothly between the two tectonic stages */
+ vip = f1*indbl+f2*indbl2;
+ }else{
+ vip = indbl;
+ }
+ E->control.surface_rayleigh[j] = vip;
+ } /* end node loop */
+ } /* end cap loop */
+ } /* end assign loop */
+ if((!timedep) && (!E->control.ggrd.ray_control_init)){ /* forget the grid */
+ ggrd_grdtrack_free_gstruc(E->control.ggrd.ray);
+ }
+ E->control.ggrd.ray_control_init = 1;
+} /* end ray control */
+
+
+
+
void ggrd_read_vtop_from_file(struct All_variables *E, int is_global)
{
MPI_Status mpi_stat;
@@ -491,7 +614,8 @@
if this file is not found, will use constant velocities
*/
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));
+ 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;
}
@@ -500,9 +624,10 @@
if(!E->control.ggrd.vtop_control_init){
if(E->parallel.me==0)
fprintf(stderr,"ggrd_read_vtop_from_file: initializing ggrd velocities\n");
- if(is_global) /* decide on GMT flag */
- sprintf(gmt_string,GGRD_GMT_XPERIODIC_STRING); /* global */
- else
+ if(is_global){ /* decide on GMT flag */
+ //sprintf(gmt_string,GGRD_GMT_XPERIODIC_STRING); /* periodic */
+ sprintf(gmt_string,GGRD_GMT_GLOBAL_STRING); /* global */
+ }else
sprintf(gmt_string,"");
/*
@@ -541,12 +666,13 @@
}
/* end init */
}
- if((E->control.ggrd.time_hist.nvtimes > 1)||(!E->control.ggrd.vtop_control_init)){
+ if((E->control.ggrd.time_hist.nvtimes > 1)||
+ (!E->control.ggrd.vtop_control_init)){
+ age = find_age_in_MY(E);
if(timedep){
/*
interpolate by time
*/
- age = find_age_in_MY(E);
if(age < 0){ /* Opposite of other method */
interpolate = 0;
/* present day should be last file*/
@@ -565,8 +691,11 @@
interpolate = 0; /* single timestep, use single file */
i1 = 0;
if(E->parallel.me == 0)
- fprintf(stderr,"ggrd_read_vtop_from_file: constant velocity BC at age %g\n",age);
+ fprintf(stderr,"ggrd_read_vtop_from_file: constant velocity BC \n");
}
+ if(E->parallel.me==0)
+ fprintf(stderr,"ggrd_read_vtop_from_file: assigning velocities BC, timedep: %i time: %g\n",
+ timedep,age);
/*
loop through all elements and assign
*/
@@ -617,7 +746,10 @@
} /* end top proc branch */
} /* end cap loop */
} /* end assignment branch */
-
+ if((!timedep)&&(!E->control.ggrd.vtop_control_init)){ /* forget the grids */
+ ggrd_grdtrack_free_gstruc(E->control.ggrd.svt);
+ ggrd_grdtrack_free_gstruc(E->control.ggrd.svp);
+ }
E->control.ggrd.vtop_control_init = 1;
}
@@ -627,6 +759,59 @@
{
myerror(E,"not implemented yet");
} /* end age control */
+
+/* adjust Ra in top boundary layer */
+void ggrd_adjust_tbl_rayleigh(struct All_variables *E,
+ double **buoy)
+{
+ int m,snode,node,i;
+ 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)
+ 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
+ */
+ for(m=1;m <= E->sphere.caps_per_proc;m++){
+ for(snode=1;snode <= E->lmesh.nsf;snode++){ /* loop through surface nodes */
+ if(fabs(E->control.surface_rayleigh[snode]-1.0)>1e-6){
+ for(i=1;i <= E->lmesh.noz;i++){ /* go through depth layers */
+ node = (snode-1)*E->lmesh.noz + i; /* global node number */
+ if(layers(E,m,node) <= E->control.ggrd.ray_control){
+ /*
+ node is in top layers
+ */
+ /* depth factor, cos^2 tapered */
+ xloc=1.0 + ((1 - E->sx[m][3][node]) -
+ E->viscosity.zbase_layer[E->control.ggrd.ray_control])/
+ E->viscosity.zbase_layer[E->control.ggrd.ray_control];
+ fac = cos(xloc*1.5707963267);fac *= fac; /* cos^2
+ tapering,
+ factor
+ decrease from
+ 1 at surface
+ to zero at
+ boundary */
+ bnew = buoy[m][node] * E->control.surface_rayleigh[snode]; /* modified rayleigh */
+ /* debugging */
+ //fprintf(stderr,"z: %11g tl: %i zm: %11g fac: %11g sra: %11g bnew: %11g bold: %11g\n",
+ // (1 - E->sx[m][3][node])*6371,E->control.ggrd.ray_control,
+ // E->viscosity.zbase_layer[E->control.ggrd.ray_control]*6371,
+ // fac,E->control.surface_rayleigh[snode],(fac * bnew + (1-fac)*buoy[m][node]),buoy[m][node]);
+ buoy[m][node] = fac * bnew + (1-fac)*buoy[m][node];
+ }
+ }
+ }
+ }
+ }
+
+}
+
+
+
#endif
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2008-04-06 20:27:32 UTC (rev 11756)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2008-04-06 21:58:45 UTC (rev 11757)
@@ -205,9 +205,12 @@
general_stokes_solver_setup(E);
+#ifdef USE_GGRD
+ if(E->control.ggrd.ray_control)
+ read_rayleigh_from_file(E);
+#endif
+
(E->next_buoyancy_field_init)(E);
-
-
if (E->parallel.me==0) fprintf(stderr,"time=%f\n",
CPU_time0()-E->monitor.cpu_time_at_start);
@@ -395,16 +398,19 @@
/* for layers */
/*
-
- these boundaries are a little wacko
-
-
+ the default boundaries are a little off
*/
input_float("z_cmb",&(E->viscosity.zcmb),"0.45",m); /* does this ever get used? */
input_float("z_lmantle",&(E->viscosity.zlm),"0.45",m);
input_float("z_410",&(E->viscosity.z410),"0.225",m); /* 0.06434, more like it */
input_float("z_lith",&(E->viscosity.zlith),"0.225",m); /* 0.0157, more like it */
+ /* to do: make layers more flexible and avoid duplication */
+ E->viscosity.zbase_layer[1] = E->viscosity.zlith;
+ E->viscosity.zbase_layer[2] = E->viscosity.z410;
+ E->viscosity.zbase_layer[3] = E->viscosity.zlm;
+ E->viscosity.zbase_layer[4] = E->viscosity.zcmb;
+
/* the start age and initial subduction history */
input_float("start_age",&(E->control.start_age),"0.0",m);
input_int("reset_startage",&(E->control.reset_startage),"0",m);
@@ -488,7 +494,15 @@
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;
+ /*
+
+ Surface layer Rayleigh number control, similar to above
+ */
+ input_int("ggrd_rayleigh_control",
+ &(E->control.ggrd.ray_control),"0",m);
+ input_string("ggrd_rayleigh_file",
+ E->control.ggrd.ray_file,"",m); /* file to read prefactors from */
/*
surface velocity control, similar to material control above
@@ -530,6 +544,7 @@
input_boolean("precond",&(E->control.precondition),"off",m);
+
input_int("mg_cycle",&(E->control.mg_cycle),"2,0,nomax",m);
input_int("down_heavy",&(E->control.down_heavy),"1,0,nomax",m);
input_int("up_heavy",&(E->control.up_heavy),"1,0,nomax",m);
@@ -540,6 +555,7 @@
input_int("piterations",&(E->control.p_iterations),"100,0,nomax",m);
input_float("rayleigh",&(E->control.Atemp),"essential",m);
+
input_float("dissipation_number",&(E->control.disptn_number),"0.0",m);
input_float("gruneisen",&(tmp),"0.0",m);
/* special case: if tmp==0, set gruneisen as inf */
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2008-04-06 20:27:32 UTC (rev 11756)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2008-04-06 21:58:45 UTC (rev 11757)
@@ -150,42 +150,49 @@
/* thermal buoyancy */
for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.nno;i++) {
- int nz = ((i-1) % E->lmesh.noz) + 1;
+ for(i=1;i<=E->lmesh.nno;i++) {
+ int nz = ((i-1) % E->lmesh.noz) + 1;
/* We don't need to substract adiabatic T profile from T here,
* since the horizontal average of buoy will be removed.
*/
buoy[m][i] = temp * E->refstate.rho[nz]
- * E->refstate.thermal_expansivity[nz] * E->T[m][i];
- }
-
+ * E->refstate.thermal_expansivity[nz] * E->T[m][i];
+ }
+
/* chemical buoyancy */
if(E->control.tracer &&
(E->composition.ichemical_buoyancy)) {
- for(j=0;j<E->composition.ncomp;j++) {
- /* TODO: how to scale chemical buoyancy wrt reference density? */
- temp2 = E->composition.buoyancy_ratio[j] * temp;
+ for(j=0;j<E->composition.ncomp;j++) {
+ /* TODO: how to scale chemical buoyancy wrt reference density? */
+ temp2 = E->composition.buoyancy_ratio[j] * temp;
for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.nno;i++)
- buoy[m][i] -= temp2 * E->composition.comp_node[m][j][i];
- }
+ for(i=1;i<=E->lmesh.nno;i++)
+ buoy[m][i] -= temp2 * E->composition.comp_node[m][j][i];
+ }
}
-
+#ifdef USE_GGRD
+ /* surface layer Rayleigh modification? */
+ if(E->control.ggrd.ray_control)
+ ggrd_adjust_tbl_rayleigh(E,buoy);
+#endif
/* phase change buoyancy */
phase_change_apply_410(E, buoy);
phase_change_apply_670(E, buoy);
phase_change_apply_cmb(E, buoy);
- /* convert density to buoyancy */
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.noz;i++)
- for(j=0;j<E->lmesh.nox*E->lmesh.noy;j++) {
- int n = j*E->lmesh.noz + i;
- buoy[m][n] *= E->refstate.gravity[i];
- }
+ /* convert density to buoyancy */
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.noz;i++)
+ for(j=0;j<E->lmesh.nox*E->lmesh.noy;j++) {
+ int n = j*E->lmesh.noz + i;
+ buoy[m][n] *= E->refstate.gravity[i];
+ }
+
+
+
remove_horiz_ave2(E,buoy);
-
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Problem_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Problem_related.c 2008-04-06 20:27:32 UTC (rev 11756)
+++ mc/3D/CitcomS/trunk/lib/Problem_related.c 2008-04-06 21:58:45 UTC (rev 11757)
@@ -44,7 +44,14 @@
(E->solver.read_input_files_for_timesteps)(E,1,1); /* read velocity(1) and output(1) */
return;
}
-
+#ifdef USE_GGRD
+void read_rayleigh_from_file(E)
+ struct All_variables *E;
+{
+ (E->solver.read_input_files_for_timesteps)(E,4,1); /* read Rayleigh number for top layers */
+ return;
+}
+#endif
/*=======================================================================
construct material array
=========================================================================*/
Modified: mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c 2008-04-06 20:27:32 UTC (rev 11756)
+++ mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c 2008-04-06 21:58:45 UTC (rev 11757)
@@ -188,6 +188,9 @@
}
#endif
break;
+ /* mode 4 is rayleigh control for GGRD, see below */
+
+
} /* end switch */
@@ -320,6 +323,15 @@
} /* end of branch if allowing for ggrd handling */
#endif
break;
+ case 4: /* material control */
+#ifdef USE_GGRD
+ if(E->control.ggrd.ray_control)
+ ggrd_read_ray_from_file(E, 0);
+#else
+ myerror(E,"input_from_files: mode 4 only for GGRD");
+#endif
+ break;
+
} /* end switch */
return;
Modified: mc/3D/CitcomS/trunk/lib/ggrd_handling.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/ggrd_handling.h 2008-04-06 20:27:32 UTC (rev 11756)
+++ mc/3D/CitcomS/trunk/lib/ggrd_handling.h 2008-04-06 21:58:45 UTC (rev 11757)
@@ -38,10 +38,11 @@
void ggrd_reg_temp_init(struct All_variables *);
void ggrd_temp_init_general(struct All_variables *,int);
void ggrd_read_mat_from_file(struct All_variables *, int );
+void ggrd_read_ray_from_file(struct All_variables *, int );
void ggrd_read_vtop_from_file(struct All_variables *, int );
void ggrd_read_age_from_file(struct All_variables *, int );
void xyz2rtp(float ,float ,float ,float *);
void xyz2rtpd(float ,float ,float ,double *);
float find_age_in_MY();
+void ggrd_adjust_tbl_rayleigh(struct All_variables *,double **);
-
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2008-04-06 20:27:32 UTC (rev 11756)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2008-04-06 21:58:45 UTC (rev 11757)
@@ -33,6 +33,9 @@
to functions across the whole filespace of CITCOM.
#include this file everywhere !
*/
+#ifdef USE_GGRD
+#include "hc.h"
+#endif
#include <assert.h>
#include <stdio.h>
@@ -41,11 +44,6 @@
-#ifdef USE_GGRD
-#include "hc.h"
-#endif
-
-
#ifdef USE_HDF5
#include "hdf5.h"
#endif
@@ -509,6 +507,7 @@
int mat_control;
#ifdef USE_GGRD
struct ggrd_master ggrd;
+ float *surface_rayleigh;
#endif
double accuracy;
double vaccuracy;
Modified: mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2008-04-06 20:27:32 UTC (rev 11756)
+++ mc/3D/CitcomS/trunk/lib/viscosity_descriptions.h 2008-04-06 21:58:45 UTC (rev 11757)
@@ -61,6 +61,7 @@
float zlm;
float z410;
float zlith;
+ float zbase_layer[40]; /* make more flexible */
int FREEZE;
float freeze_thresh;
@@ -125,9 +126,13 @@
char VISC_OPT[20];
int layers; /* number of layers with properties .... */
- float layer_depth[40];
- float layer_visc[40];
+ /* those were unused?, see above with zbase_layer
+ which should take over
+ */
+ //float layer_depth[40];
+ //float layer_visc[40];
+
int SLABLVZ; /* slab structure imposed on top of 3 layer structure */
int slvzd1,slvzd2,slvzd3; /* layer thicknesses (nodes) */
int slvzD1,slvzD2; /* slab posn & length */
More information about the cig-commits
mailing list