[cig-commits] r11268 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Tue Feb 26 23:42:44 PST 2008
Author: becker
Date: 2008-02-26 23:42:44 -0800 (Tue, 26 Feb 2008)
New Revision: 11268
Added:
mc/3D/CitcomS/trunk/lib/ggrd_handling.h
Modified:
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/Lith_age.c
mc/3D/CitcomS/trunk/lib/Output_gzdir.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c
mc/3D/CitcomS/trunk/lib/Tracer_setup.c
Log:
Built surface velocity BC support via Netcdf grd files back in. This
needs more testing, but it would be good to have the framework in
place before new additions are made. Also, age control is not yet
implemented, but will be soon.
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-27 07:24:31 UTC (rev 11267)
+++ mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c 2008-02-27 07:42:44 UTC (rev 11268)
@@ -29,6 +29,9 @@
#include <sys/types.h>
#include "element_definitions.h"
#include "global_defs.h"
+#ifdef USE_GGRD
+#include "ggrd_handling.h"
+#endif
/*=======================================================================
Calculate ages (MY) for opening input files -> material, ages, velocities
@@ -94,6 +97,9 @@
switch (action) { /* set up files to open */
case 1: /* read velocity boundary conditions */
+#ifdef USE_GGRD
+ if(!E->control.ggrd.vtop_control){
+#endif
sprintf(output_file1,"%s%0.0f.%d",E->control.velocity_boundary_file,newage1,cap);
sprintf(output_file2,"%s%0.0f.%d",E->control.velocity_boundary_file,newage2,cap);
fp1=fopen(output_file1,"r");
@@ -116,9 +122,15 @@
else
fprintf(E->fp,"Velocity: File2 = No file inputted (negative age)\n");
}
+#ifdef USE_GGRD
+ }
+#endif
break;
case 2: /* read ages for lithosphere temperature assimilation */
+#ifdef USE_GGRD
+ if(!E->control.ggrd.age_control){
+#endif
sprintf(output_file1,"%s%0.0f.%d",E->control.lith_age_file,newage1,cap);
sprintf(output_file2,"%s%0.0f.%d",E->control.lith_age_file,newage2,cap);
fp1=fopen(output_file1,"r");
@@ -141,6 +153,9 @@
else
fprintf(E->fp,"Age: File2 = No file inputted (negative age)\n");
}
+#ifdef USE_GGRD
+ }
+#endif
break;
case 3: /* read element materials */
@@ -181,6 +196,11 @@
switch (action) { /* Read the contents of files and average */
case 1: /* velocity boundary conditions */
+#ifdef USE_GGRD
+ if(E->control.ggrd.vtop_control){
+ ggrd_read_vtop_from_file(E, 1);
+ }else{
+#endif
nnn=nox*noy;
for(i=1;i<=dims;i++) {
VB1[i]=(float*) malloc ((nnn+1)*sizeof(float));
@@ -222,9 +242,17 @@
free ((void *) VB1[i]);
free ((void *) VB2[i]);
}
+#ifdef USE_GGRD
+ }
+#endif
break;
case 2: /* ages for lithosphere temperature assimilation */
+#ifdef USE_GGRD
+ if(E->control.ggrd.age_control){
+ ggrd_read_age_from_file(E, 1);
+ }else{
+#endif
for(i=1;i<=noy;i++)
for(j=1;j<=nox;j++) {
node=j+(i-1)*nox;
@@ -239,6 +267,9 @@
}
fclose(fp1);
if (pos_age) fclose(fp2);
+#ifdef USE_GGRD
+ } /* end of branch if allowing for ggrd handling */
+#endif
break;
case 3: /* read element materials */
Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2008-02-27 07:24:31 UTC (rev 11267)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c 2008-02-27 07:42:44 UTC (rev 11268)
@@ -42,17 +42,8 @@
#ifdef USE_GGRD
#include "hc.h" /* ggrd and hc packages */
+#include "ggrd_handling.h"
-void ggrd_init_tracer_flavors(struct All_variables *);
-int layers_r(struct All_variables *,float );
-void myerror(struct All_variables *,char *);
-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();
-
/*
assign tracer flavor based on its depth (within top n layers),
@@ -78,7 +69,7 @@
*/
if (E->parallel.nprocxy == 12){
/* use GMT's geographic boundary conditions */
- sprintf(gmt_bc,"-Lg");
+ sprintf(gmt_bc,GGRD_GMT_GLOBAL_STRING);
}else{ /* regional */
sprintf(gmt_bc,"");
}
@@ -147,11 +138,11 @@
void ggrd_full_temp_init(struct All_variables *E)
{
- ggrd_temp_init_general(E,"-L");
+ ggrd_temp_init_general(E,1);
}
void ggrd_reg_temp_init(struct All_variables *E)
{
- ggrd_temp_init_general(E,"");
+ ggrd_temp_init_general(E,0);
}
@@ -162,15 +153,19 @@
*/
-void ggrd_temp_init_general(struct All_variables *E,char *gmtflag)
+void ggrd_temp_init_general(struct All_variables *E,int is_global)
{
MPI_Status mpi_stat;
int mpi_rc;
int mpi_inmsg, mpi_success_message = 1;
double temp1,tbot,tgrad,tmean,tadd,rho_prem;
-
+ char gmt_string[10];
int i,j,k,m,node,noxnoz,nox,noy,noz;
+ if(is_global) /* decide on GMT flag */
+ sprintf(gmt_string,GGRD_GMT_GLOBAL_STRING); /* global */
+ else
+ sprintf(gmt_string,"");
noy=E->lmesh.noy;
nox=E->lmesh.nox;
@@ -178,7 +173,7 @@
noxnoz = nox * noz;
if(E->parallel.me == 0)
- fprintf(stderr,"ggrd_temp_init_general: using GMT grd files for temperatures, gmtflag: %s\n",gmtflag);
+ fprintf(stderr,"ggrd_temp_init_general: using GMT grd files for temperatures, gmtflag: %s\n",gmt_string);
/*
@@ -209,7 +204,7 @@
*/
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.dfile,gmt_string,
E->control.ggrd.temp_init.d,(E->parallel.me == 0),
FALSE))
myerror(E,"grd init error");
@@ -314,41 +309,27 @@
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_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,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];
+ 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];
-
-
-
- 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;
-
+ 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
*/
@@ -360,16 +341,15 @@
E->control.ggrd.time_hist.file,TRUE,(E->parallel.me == 0));
E->control.ggrd.time_hist.init = 1;
}
+ /* time dependent? */
+ timedep = (E->control.ggrd.time_hist.nvtimes > 1)?(1):(0);
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 */
+ sprintf(gmt_string,GGRD_GMT_GLOBAL_STRING); /* global */
else
sprintf(gmt_string,"");
/*
@@ -385,7 +365,7 @@
*/
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)
+ if(!timedep) /* constant */
sprintf(tfilename,"%s",E->control.ggrd.mat_file);
else
sprintf(tfilename,"%s/%i/weak.grd",E->control.ggrd.mat_file,i+1);
@@ -404,7 +384,16 @@
/* end init */
}
- if((E->control.ggrd.time_hist.nvtimes > 1)||(!E->control.ggrd.mat_control_init)){
+ if(timedep || (!E->control.ggrd.mat_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;
+ }
/*
loop through all elements and assign
*/
@@ -430,38 +419,24 @@
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)){
+ xyz2rtpd(xloc[1],xloc[2],xloc[3],rout);
+ /* 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]);
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)){
+ }
+ 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]);
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);
+ }else{
+ vip = indbl;
}
/* limit the input scaling? */
if(vip < 1e-5)
@@ -486,8 +461,172 @@
} /* end assignment loop */
E->control.ggrd.mat_control_init = 1;
+} /* end mat control */
+
+
+
+void ggrd_read_vtop_from_file(struct All_variables *E, int is_global)
+{
+ MPI_Status mpi_stat;
+ int mpi_rc,interpolate,timedep;
+ int mpi_inmsg, mpi_success_message = 1;
+ int m,el,i,k,i1,i2,ind,nodel;
+ int nox1,noz1,noy1,lev,lselect,idim,nox1noz1;
+ char gmt_string[10],char_dummy;
+
+ double vin1[2],vin2[2],age,f1,f2,vp,vt,vscale,rout[3];
+ char tfilename1[1000],tfilename2[1000];
+
+ const int dims=E->mesh.nsd;
+
+ nox1 = E->lmesh.nox;noz1=E->lmesh.noz;noy1=E->lmesh.noy;
+ nox1noz1 = nox1*noz1;
+ lev = E->mesh.levmax;
+
+ /* velocity scaling, assuming input is cm/yr */
+ vscale = E->data.scalev * E->data.timedir;
+
+ /*
+ if we have not initialized the time history structure, do it now
+ 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));
+ E->control.ggrd.time_hist.init = 1;
+ }
+
+ timedep = (E->control.ggrd.time_hist.nvtimes > 1)?(1):(0);
+
+ 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
+ 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 velocity file(s)
+ */
+ E->control.ggrd.svt = (struct ggrd_gt *)calloc(E->control.ggrd.time_hist.nvtimes,sizeof(struct ggrd_gt));
+ E->control.ggrd.svp = (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(tfilename1,"%s/vt.grd",E->control.ggrd.vtop_dir);
+ sprintf(tfilename2,"%s/vp.grd",E->control.ggrd.vtop_dir);
+ } else { /* f(t) */
+ sprintf(tfilename1,"%s/%i/vt.grd",E->control.ggrd.vtop_dir,i+1);
+ sprintf(tfilename2,"%s/%i/vp.grd",E->control.ggrd.vtop_dir,i+1);
+ }
+ if(ggrd_grdtrack_init_general(FALSE,tfilename1,&char_dummy,
+ gmt_string,(E->control.ggrd.svt+i),(E->parallel.me == 0),FALSE))
+ myerror(E,"ggrd init error vt");
+ if(ggrd_grdtrack_init_general(FALSE,tfilename2,&char_dummy,
+ gmt_string,(E->control.ggrd.svp+i),(E->parallel.me == 0),FALSE))
+ myerror(E,"ggrd init error vp");
+ }
+ 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_vtop_from_file: last processor done with ggrd vtop BC init, %i timesteps\n",
+ E->control.ggrd.time_hist.nvtimes);
+ }
+ /* end init */
+ }
+ if((E->control.ggrd.time_hist.nvtimes > 1)||(!E->control.ggrd.vtop_control_init)){
+ 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*/
+ i1 = E->control.ggrd.time_hist.nvtimes - 1;
+ if(E->parallel.me == 0)
+ fprintf(stderr,"ggrd_read_vtop_from_file: using present day vtop for age = %g\n",age);
+ }else{
+ /* */
+ ggrd_interpol_time(age,&E->control.ggrd.time_hist,&i1,&i2,&f1,&f2,
+ E->control.ggrd.time_hist.vstage_transition);
+ interpolate = 1;
+ if(E->parallel.me == 0)
+ fprintf(stderr,"ggrd_read_vtop_from_file: interpolating vtop for age = %g\n",age);
+ }
+ }else{
+ 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);
+ }
+ /*
+ loop through all elements and assign
+ */
+ for (m=1;m <= E->sphere.caps_per_proc;m++) {
+ if(E->parallel.me_loc[3] == E->parallel.nprocz-1 ) { /* if on top */
+ for(k=1;k <= noy1;k++) {/* loop through surface nodes */
+ for(i=1;i <= nox1;i++) {
+ nodel = (k-1)*nox1noz1 + (i-1)*noz1+noz1;
+
+ rout[1] = E->sx[m][1][nodel]; /* theta,phi coordinates */
+ rout[2] = E->sx[m][2][nodel];
+
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.svt+i1),
+ vin1,FALSE)){
+ fprintf(stderr,"ggrd_read_vtop_from_ggrd_file: interpolation error at %g, %g\n",
+ rout[1],rout[2]);
+ parallel_process_termination();
+ }
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.svp+i1),
+ (vin1+1),FALSE)){
+ fprintf(stderr,"ggrd_read_vtop_from_ggrd_file: interpolation error at %g, %g\n",
+ rout[1],rout[2]);
+ parallel_process_termination();
+ }
+ if(interpolate){ /* second time */
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.svt+i2),vin2,FALSE)){
+ fprintf(stderr,"ggrd_read_mat_from_ggrd_file: interpolation error at %g, %g\n",
+ rout[1],rout[2]);
+ parallel_process_termination();
+ }
+ if(!ggrd_grdtrack_interpolate_tp(rout[1],rout[2],(E->control.ggrd.svp+i2),(vin2+1),FALSE)){
+ fprintf(stderr,"ggrd_read_mat_from_ggrd_file: interpolation error at %g, %g\n",
+ rout[1],rout[2]);
+ parallel_process_termination();
+ }
+ vt = (f1*vin1[0] + f2*vin2[0])*vscale;
+ vp = (f1*vin1[1] + f2*vin2[1])*vscale;
+ }else{
+ vt = vin1[0]*vscale;
+ vp = vin1[1]*vscale;
+ }
+ /* assign as boundary condition */
+ E->sphere.cap[m].VB[1][nodel] = vt; /* theta */
+ E->sphere.cap[m].VB[2][nodel] = vp; /* phi */
+ E->sphere.cap[m].VB[3][nodel] = 0.0; /* r */
+ }
+ } /* end surface node loop */
+ } /* end top proc branch */
+ } /* end cap loop */
+ } /* end assignment branch */
+
+ E->control.ggrd.vtop_control_init = 1;
}
+
+
+void ggrd_read_age_from_file(struct All_variables *E, int is_global)
+{
+ myerror(E,"not implemented yet");
+} /* end age control */
#endif
Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2008-02-27 07:24:31 UTC (rev 11267)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2008-02-27 07:42:44 UTC (rev 11268)
@@ -43,6 +43,9 @@
#ifdef USE_GZDIR
void restart_tic_from_gzdir_file(struct All_variables *);
#endif
+#ifdef USE_GGRD
+#include "ggrd_handling.h"
+#endif
void tic_input(struct All_variables *E)
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2008-02-27 07:24:31 UTC (rev 11267)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2008-02-27 07:42:44 UTC (rev 11268)
@@ -444,7 +444,14 @@
input_string("mat_file",E->control.mat_file,"",m);
#ifdef USE_GGRD
+
+
/*
+
+ note that this part of the code might override mat_control, file_vbcs,
+
+ MATERIAL CONTROL
+
usage:
(a)
@@ -460,21 +467,44 @@
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
+ 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
+
+ -->age is positive, and forward marching in time decreases the age<--
- -60 -30
- -30 -15
- -15 0
+ 0 15
+ 15 30
+ 30 60
*/
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 */
+ /* this is controlling velocities, material, and age */
+ /* time history file, if not specified, will use constant VBCs and material grids */
+ input_string("ggrd_time_hist_file",
+ E->control.ggrd.time_hist.file,"",m);
+ /* if > 0, will use top E->control.ggrd.mat_control layers and assign a prefactor for the viscosity */
+ input_int("ggrd_mat_control",&(E->control.ggrd.mat_control),"0",m);
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 velocity control, similar to material control above
+
+ if time-dependent, will look for ggrd_vtop_file/i/v?.grd
+ if constant, will look for ggrd_vtop_file/v?.grd
+
+ where vp/vt.grd are Netcdf GRD files with East and South velocities in cm/yr
+
+
+ */
+ input_int("ggrd_vtop_control",&(E->control.ggrd.vtop_control),"0",m);
+ input_string("ggrd_vtop_dir",E->control.ggrd.vtop_dir,"",m); /* file to read prefactors from */
+ if(E->control.ggrd.vtop_control) /* this will override mat_control setting */
+ E->control.vbcs_file = 1;
+
#endif
input_int("nodex",&(E->mesh.nox),"essential",m);
Modified: mc/3D/CitcomS/trunk/lib/Lith_age.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Lith_age.c 2008-02-27 07:24:31 UTC (rev 11267)
+++ mc/3D/CitcomS/trunk/lib/Lith_age.c 2008-02-27 07:42:44 UTC (rev 11268)
@@ -48,6 +48,13 @@
E->control.temperature_bound_adj = 0;
input_int("lith_age",&(E->control.lith_age),"0",m);
+#ifdef USE_GGRD
+ input_int("ggrd_age_control",&(E->control.ggrd.age_control),"0",m); /* if > 0, will use top E->control.ggrd.mat_control layers and assign a prefactor for the viscosity */
+ if(E->control.ggrd.age_control){
+ E->control.lith_age = 1;
+ }
+#endif
+
input_float("mantle_temp",&(E->control.lith_age_mantle_temp),"1.0",m);
if (E->control.lith_age) {
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2008-02-27 07:24:31 UTC (rev 11267)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2008-02-27 07:42:44 UTC (rev 11268)
@@ -901,7 +901,7 @@
for(n=0; n<E->composition.ncomp; n++)
gzprintf(fp1," %.4e", E->Have.C[n][j]);
}
- fprintf(fp1,"\n");
+ gzprintf(fp1,"\n");
}
gzclose(fp1);
}
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2008-02-27 07:24:31 UTC (rev 11267)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2008-02-27 07:42:44 UTC (rev 11268)
@@ -57,6 +57,7 @@
void *safe_malloc (size_t );
void myerror(struct All_variables *,char *);
void xyz2rtp(float ,float ,float ,float *);
+void xyz2rtpd(float ,float ,float ,double *);
void get_r_spacing_fine(double *,struct All_variables *);
void get_r_spacing_at_levels(double *,struct All_variables *);
@@ -422,6 +423,15 @@
rout[1] = atan2(sqrt(tmp1),z); /* theta */
rout[2] = atan2(y,x); /* phi */
}
+void xyz2rtpd(float x,float y,float z,double *rout)
+{
+ double tmp1,tmp2;
+ tmp1 = (double)x*(double)x + (double)y*(double)y;
+ tmp2 = tmp1 + (double)z*(double)z;
+ rout[0] = sqrt(tmp2); /* r */
+ rout[1] = atan2(sqrt(tmp1),(double)z); /* theta */
+ rout[2] = atan2((double)y,(double)x); /* phi */
+}
/* compute base vectors for conversion of polar to cartesian vectors
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-27 07:24:31 UTC (rev 11267)
+++ mc/3D/CitcomS/trunk/lib/Regional_read_input_from_files.c 2008-02-27 07:42:44 UTC (rev 11268)
@@ -29,6 +29,9 @@
#include <sys/types.h>
#include "element_definitions.h"
#include "global_defs.h"
+#ifdef USE_GGRD
+#include "ggrd_handling.h"
+#endif
/*=======================================================================
Calculate ages (MY) for opening input files -> material, ages, velocities
@@ -65,6 +68,7 @@
nox=E->mesh.nox;
noy=E->mesh.noy;
noz=E->mesh.noz;
+
nox1=E->lmesh.nox;
noz1=E->lmesh.noz;
noy1=E->lmesh.noy;
@@ -78,7 +82,8 @@
age=find_age_in_MY(E);
- if (age < 0.0) { /* age is negative -> use age=0 for input files */
+ if (age < 0.0) { /* age is negative -> use age=0 for input
+ files */
intage = 0;
newage2 = newage1 = 0.0;
pos_age = 0;
@@ -91,8 +96,10 @@
}
switch (action) { /* set up files to open */
-
case 1: /* read velocity boundary conditions */
+#ifdef USE_GGRD
+ if(!E->control.ggrd.vtop_control){ /* regular input */
+#endif
sprintf(output_file1,"%s%0.0f",E->control.velocity_boundary_file,newage1);
sprintf(output_file2,"%s%0.0f",E->control.velocity_boundary_file,newage2);
fp1=fopen(output_file1,"r");
@@ -115,9 +122,15 @@
else
fprintf(E->fp,"Velocity: File2 = No file inputted (negative age)\n");
}
+#ifdef USE_GGRD
+ }
+#endif
break;
case 2: /* read ages for lithosphere temperature assimilation */
+#ifdef USE_GGRD
+ if(!E->control.ggrd.age_control){ /* regular input */
+#endif
sprintf(output_file1,"%s%0.0f",E->control.lith_age_file,newage1);
sprintf(output_file2,"%s%0.0f",E->control.lith_age_file,newage2);
fp1=fopen(output_file1,"r");
@@ -139,6 +152,9 @@
else
fprintf(E->fp,"Age: File2 = No file inputted (negative age)\n");
}
+#ifdef USE_GGRD
+ }
+#endif
break;
case 3: /* read element materials */
@@ -179,6 +195,11 @@
switch (action) { /* Read the contents of files and average */
case 1: /* velocity boundary conditions */
+#ifdef USE_GGRD
+ if(E->control.ggrd.vtop_control){
+ ggrd_read_vtop_from_file(E, 0);
+ }else{
+#endif
nnn=nox*noy;
for(i=1;i<=dims;i++) {
VB1[i]=(float*) malloc ((nnn+1)*sizeof(float));
@@ -218,9 +239,19 @@
free ((void *) VB1[i]);
free ((void *) VB2[i]);
}
+
+
+#ifdef USE_GGRD
+ } /* end of branch if allowing for ggrd handling */
+#endif
break;
case 2: /* ages for lithosphere temperature assimilation */
+#ifdef USE_GGRD
+ if(E->control.ggrd.age_control){
+ ggrd_read_age_from_file(E, 0);
+ }else{
+#endif
for(i=1;i<=noy;i++)
for(j=1;j<=nox;j++) {
node=j+(i-1)*nox;
@@ -235,6 +266,9 @@
}
fclose(fp1);
if (pos_age) fclose(fp2);
+#ifdef USE_GGRD
+ } /* end of branch if allowing for ggrd handling */
+#endif
break;
case 3: /* read element materials */
Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2008-02-27 07:24:31 UTC (rev 11267)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2008-02-27 07:42:44 UTC (rev 11268)
@@ -46,8 +46,9 @@
#include "composition_related.h"
#ifdef USE_GGRD
-void ggrd_init_tracer_flavors(struct All_variables *);
+#include "ggrd_handling.h"
#endif
+
#ifdef USE_GZDIR
int open_file_zipped(char *, FILE **,struct All_variables *);
void gzip_file(char *);
Added: mc/3D/CitcomS/trunk/lib/ggrd_handling.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/ggrd_handling.h (rev 0)
+++ mc/3D/CitcomS/trunk/lib/ggrd_handling.h 2008-02-27 07:42:44 UTC (rev 11268)
@@ -0,0 +1,47 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<LicenseText>
+ *
+ * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
+ * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
+ * Copyright (C) 1994-2005, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+/*
+
+routines that deal with GMT/netcdf grd I/O as supported through
+the ggrd subroutines of the hc package
+
+*/
+void ggrd_init_tracer_flavors(struct All_variables *);
+int layers_r(struct All_variables *,float );
+void myerror(struct All_variables *,char *);
+void ggrd_full_temp_init(struct All_variables *);
+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_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();
+
+
More information about the cig-commits
mailing list