[cig-commits] r8044 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Thu Sep 27 15:43:59 PDT 2007
Author: becker
Date: 2007-09-27 15:43:58 -0700 (Thu, 27 Sep 2007)
New Revision: 8044
Added:
mc/1D/hc/trunk/grdinttester.c
Modified:
mc/1D/hc/trunk/.gmtcommands
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/calc_vel_and_plot
mc/1D/hc/trunk/ggrd.h
mc/1D/hc/trunk/ggrd_grdtrack_util.c
mc/1D/hc/trunk/ggrd_grdtrack_util.h
mc/1D/hc/trunk/ggrd_struc.h
mc/1D/hc/trunk/ggrd_test.c
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_extract_sh_layer.c
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/hc_input.c
mc/1D/hc/trunk/hc_output.c
mc/1D/hc/trunk/hc_polsol.c
mc/1D/hc/trunk/hc_solve.c
mc/1D/hc/trunk/hc_torsol.c
mc/1D/hc/trunk/main.c
mc/1D/hc/trunk/make_tar
mc/1D/hc/trunk/rick_sh_c.c
Log:
Made several modifications that were needed to keep up to date
with GMT 4 development.
Modified: mc/1D/hc/trunk/.gmtcommands
===================================================================
--- mc/1D/hc/trunk/.gmtcommands 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/.gmtcommands 2007-09-27 22:43:58 UTC (rev 8044)
@@ -1,7 +1,7 @@
# GMT common arguments shelf
--B1/:v at -r@- [cm/yr]:
+-B.5/:v at -r@- [cm/yr]:
-JH180/7i
--JX7/3.5i
+-JQ180/7
+-JX180/7
-R0/360/-90/90
EOF
-EOF
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/Makefile 2007-09-27 22:43:58 UTC (rev 8044)
@@ -9,19 +9,19 @@
#
#
LD = $(CC)
-#CFLAGS = $(CFLAGS_DEBUG)
+CFLAGS = $(CFLAGS_DEBUG)
#
# EDIT HERE FOR GMT VERSION
#
#
# for GMT3.4.5, use the next two lines
-GGRD_INC_FLAGS = -I$(GMTHOME)/include -I$(NETCDFHOME)/include
-GGRD_LIBS_LINKLINE = -lggrd -lgmt -lnetcdf
+#GGRD_INC_FLAGS = -I$(GMTHOME)/include -I$(NETCDFHOME)/include
+#GGRD_LIBS_LINKLINE = -lggrd -lgmt -lnetcdf
#
# for GMT version >= 4.1.2, uncomment the next two lines
-#GGRD_INC_FLAGS = -I$(GMTHOME)/include -I$(NETCDFHOME)/include -DUSE_GMT4
-#GGRD_LIBS_LINKLINE = -lggrd -lgmt -lpsl -lnetcdf
+GGRD_INC_FLAGS = -I$(GMTHOME)/include -I$(NETCDFHOME)/include -DUSE_GMT4
+GGRD_LIBS_LINKLINE = -lggrd -lgmt -lpsl -lnetcdf
#
#
#
@@ -87,7 +87,7 @@
$(PREM_DEFINES)
GGRD_LIB_FLAGS = -L$(GMTHOME)/lib -L$(NETCDFHOME)/lib
GGRD_LIBS = $(ODIR)/libggrd.a $(ODIR)/libggrd.dfast.a $(ODIR)/libggrd.dbg.a
-GGRD_INCS = $(PREM_INCS)
+GGRD_INCS = $(PREM_INCS) ggrd_grdtrack_util.h ggrd.h ggrd_struc.h
#
#
@@ -149,7 +149,7 @@
debug_libs: $(HC_LIBS_DEBUG)
-really_all: proto all debug_libs hc.dbg hcplates ggrd_test
+really_all: proto all debug_libs hc.dbg hcplates ggrd_test grdinttester
proto: hc_auto_proto.h
@@ -194,6 +194,10 @@
$(LD) $(LIB_FLAGS) $(ODIR)/ggrd_test.o -o $(BDIR)/ggrd_test \
$(GGRD_LIBS_LINKLINE) -lhc -lrick -lm $(LDFLAGS)
+grdinttester: $(LIBS) $(INCS) $(ODIR)/grdinttester.o
+ $(LD) $(LIB_FLAGS) $(ODIR)/grdinttester.o -o $(BDIR)/grdinttester \
+ $(GGRD_LIBS_LINKLINE) -lhc -lrick -lm $(LDFLAGS)
+
hc_extract_sh_layer: $(LIBS) $(INCS) $(ODIR)/hc_extract_sh_layer.o
$(LD) $(LIB_FLAGS) $(ODIR)/hc_extract_sh_layer.o \
-o $(BDIR)/hc_extract_sh_layer \
Modified: mc/1D/hc/trunk/calc_vel_and_plot
===================================================================
--- mc/1D/hc/trunk/calc_vel_and_plot 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/calc_vel_and_plot 2007-09-27 22:43:58 UTC (rev 8044)
@@ -4,9 +4,10 @@
#
compute=${1-1} # produce a new solution
-plots=${2-"1 2"} # plot output modes
+plots=${2-"1 2 3"} # plot output modes
# 1: geoid
# 2: plate velocities
+ # 3: tractions
#
# options
@@ -27,11 +28,14 @@
# viscosity file
#
vf=example_data/visc.C2 # file name
+#vf=example_data/visc.A # file name
#vf=example_data/visc.dat # file name
verb=-vvv # verbosity level
+layers="21 15 6 2" # layers for plotting
+
#
# for temporary files
#
@@ -40,20 +44,31 @@
-if [ $compute -eq 1 ];then
- rm vel.sol.bin 2> /dev/null
- hc $verb -dens $dm $dt -ds $dfac $tbc -vf $vf 2> hc.log
- cat hc.log
- if [ ! -s vel.sol.bin ];then
- echo $0: error, no output produced
- exit
+for plot in $plots;do
+
+
+ if [[ $plot -eq 1 || $plot -eq 2 ]];then # velocities or geoid
+ f=vel
+
+ hc_opt=""
+ dfile=vdepth
+ elif [ $plot -eq 3 ];then # radial tractions
+ f=rtrac
+ hc_opt="-rtrac" #
+ dfile=sdepth
fi
-fi
+ if [ $compute -eq 1 ];then
+ rm $f.sol.bin 2> /dev/null
+ hc $verb -dens $dm $dt $hc_opt -ds $dfac $tbc -vf $vf 2> hc.log
+ if [ ! -s $f.sol.bin ];then
+ echo $0: error, no output $f.sol.bin produced
+ exit
+ fi
+ cat hc.log
+ fi
-for plot in $plots;do
-
#
# take care of GMT crap
#
@@ -77,75 +92,117 @@
#
# plotting setting
#
- inc=1;greg=-R0/359/-89.5/89.5 # grid range
+ inc=1;
+ greg=-R0/359/-89.5/89.5 # grid range
proj=-JH180/7i;preg=-R0/360/-90/90 # plotting range
if [ $plot -eq 1 ];then # plot the geoid
- ofile=geoid.ps
+ ofile=geoid.ps # output PS filename
echo $0: plotting geoid to $ofile
- makecpt -T-120/120/10 -Chaxby > $tmpn.cpt
+ makecpt -T-150/150/10 -Chaxby > $tmpn.cpt # makes colormap
+
+ # go from spherical harmonics to spatial
cat geoid.ab | sh_syn 0 0 $inc 2> /dev/null | \
- xyz2grd $greg -I$inc -G$tmpn.grd
- grdimage $tmpn.grd $preg $proj -P \
- -C$tmpn.cpt -K > $ofile
- pscoast $proj $preg -A70000 -K -O -Dc -W2 >> $ofile
- psscale -D50 -D3.5i/-.3i/3i/.25ih -E -C$tmpn.cpt \
- -B100/:"[m]": -O >> $ofile
+ xyz2grd $greg -I$inc -G$tmpn.grd # and convert to netcdf GRD
+
+ # make postscript plot from grd file
+ grdimage $tmpn.grd $preg $proj -P -C$tmpn.cpt -K > $ofile # first: -K
+ # coastline
+ pscoast $proj $preg -A70000 -K -O -Dc -W2 >> $ofile # all intermediat: -O -K
+ # colorbar
+ psscale -N50 -D3.5i/-.3i/3i/.25ih -E -C$tmpn.cpt -B50/:"[m]": -O >> $ofile # end: -O
+
+ # ps file done
- elif [ $plot -eq 2 ];then # plot velocities
+ elif [[ $plot -eq 2 || $plot -eq 3 ]];then # plot velocities or tractions
- # extract all depths
- hc_extract_sh_layer vel.sol.bin 2 4 | gawk '{print($2)}' > vdepth.dat
- nl=`wc -l vdepth.dat | gawk '{print($1)}'`
+ # extract all depths , use second colume, print depth in [km]
+ hc_extract_sh_layer $f.sol.bin 2 4 | gawk '{print($2)}' > $dfile.dat
+ # count number of lines (assign command output to variable)
+ nl=`wc -l $dfile.dat | gawk '{print($1)}'`
+
+ if [ $plot -eq 2 ];then
# colormaps for radial flow
- makecpt -T-1.5/1.5/.125 -Chaxby > $tmpn.cpt
+ makecpt -T-1.5/1.5/.125 -Chaxby > $tmpn.cpt
+ scsp=.5 # scale spacing
+ type="v";units="cm/yr"
+ else # tractions
+ scsp=50
+ type="t";units="MPa"
+ makecpt -T-150/150/10 -Chaxby > $tmpn.cpt
+ fi
# vector spacing
vinc=10;vreg=-R0/350/-85/85
+ #
# loop through some layers
- layer=2
- while [ $layer -le $nl ];do
+ #
+ lc=1
+ for layer in $layers;do # loop
+
+ if [ $layer -gt $nl ];then # check
+ echo $0: error, layer $layer out of bounds $nl
+ exit
+ fi
+ #
# extract depth
- z=`gawk '{if(NR==l)print($1)}' l=$layer vdepth.dat `
+ z=`gawk '{if(NR==l)printf("%.0f",$1)}' l=$layer $dfile.dat `
- ofile=vel.$layer.ps # output file name
+ ofile=$f.$layer.ps # output file name
- echo $0: plotting velocities at $z to $ofile
+ echo $0: plotting $f at $z to $ofile
#
# extract radial velocity of this layer
#
- hc_extract_sh_layer vel.sol.bin $layer 1 0 | sh_syn 0 0 $inc 2> /dev/null | \
+ # extract layer $layer as spherical harmonics , radial component
+ # convert to spatial
+ # make netcdf GRD
+ hc_extract_sh_layer $f.sol.bin $layer 1 0 | sh_syn 0 0 $inc 2> /dev/null | \
xyz2grd $greg -I$inc -G$tmpn.grd
#
# extract vx and vy velocities
- hc_extract_sh_layer vel.sol.bin $layer 2 0 | sh_syn 0 0 $vinc 2> /dev/null > $tmpn.dat
- gawk '{print($1,$2,$4)}' $tmpn.dat | xyz2grd -G$tmpn.vx.grd $vreg -I$vinc # vx = vphi
- gawk '{print($1,$2,-$3)}' $tmpn.dat | xyz2grd -G$tmpn.vy.grd $vreg -I$vinc # vy = -vtheta
+ #
+ hc_extract_sh_layer $f.sol.bin $layer 2 0 | sh_syn 0 0 $vinc 2> /dev/null > $tmpn.dat
+ #
+ gawk '{print($1,$2,$4)}' $tmpn.dat | xyz2grd -G$tmpn.vx.grd $vreg -I$vinc # vx = vphi, tphi
+ gawk '{print($1,$2,-$3)}' $tmpn.dat | xyz2grd -G$tmpn.vy.grd $vreg -I$vinc # vy = -vtheta,-ttheta
+ # computing vector length
+ grdmath $tmpn.vx.grd $tmpn.vy.grd R2 SQRT = $tmpn.abs.grd
+
+ if [ $lc -eq 1 ];then # first time around
+ #
+ # get scale from multiples of mean of first layer
+ hmean=`grdinfo -C -L2 $tmpn.abs.grd | gawk '{printf("%.1f",$12)}'` # compute mean
+ scale=`echo $hmean | gawk '{print($1*5)}'` # scale
+ fi
+ #
# actually plot velocities
-
- grdimage $tmpn.grd $preg $proj -P -C$tmpn.cpt -K > $ofile
+ #
+ grdimage $tmpn.grd $preg $proj -P -C$tmpn.cpt -K > $ofile # make a colormap
+ # coastline
pscoast $proj $preg -A70000 -K -O -Dc -W2 >> $ofile
+
+ #
echo 0.05 0 14 0 0 ML "z = $z km" | pstext -O -K -N -R0/1/0/1 -JX7/3.5i >> $ofile
+ echo 0.75 0 14 0 0 ML "@~\341@~$type at -h@-@~\361@~ = $hmean $units" | \
+ pstext -O -K -N -R0/1/0/1 -JX7/3.5i >> $ofile
+ # plot a vector field
grdvector $tmpn.vx.grd $tmpn.vy.grd \
- -T $reg $proj -Q0.015i/0.12i/0.045in.2 -S15 -O -K -G128/0/0 -W0.5 >> $ofile
+ -T $reg $proj -Q0.015i/0.12i/0.045in.2 -S$scale -O -K -G128/0/0 -W0.5 >> $ofile
+ # scale
psscale -D50 -D3.5i/-.3i/3i/.25ih -E -C$tmpn.cpt \
- -B1/:"v at -r@- [cm/yr]": -O >> $ofile
+ -B$scsp/:"$type at -r@- [$units]": -O >> $ofile
-
-
- ((layer=layer+6))
+ ((lc=lc+1))
done
-
-
-
fi
if [ -s $gf ];then # resset GMT
gmtset PAPER_MEDIA $old_media
Modified: mc/1D/hc/trunk/ggrd.h
===================================================================
--- mc/1D/hc/trunk/ggrd.h 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/ggrd.h 2007-09-27 22:43:58 UTC (rev 8044)
@@ -1,4 +1,7 @@
#ifndef __READ_GGRD_HEADER__
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
//
// max order of interpolation - 1
@@ -40,7 +43,10 @@
#define GGRD_TWOPI 6.283185307179586476925286766559005768394
#endif
-
+/* 180/pi */
+#ifndef GGRD_PIF
+#define GGRD_PIF 57.295779513082320876798154814105
+#endif
/*
modes
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -12,6 +12,7 @@
*/
#include "gmt.h"
+#include "gmt_bcr.h"
#include "ggrd_grdtrack_util.h"
#ifndef ONEEIGHTYOVERPI
#define ONEEIGHTYOVERPI 57.295779513082320876798154814105
@@ -19,6 +20,7 @@
#include <math.h>
#include <string.h>
+
/*
wrapper
@@ -59,21 +61,14 @@
g->east = g->west = g->south = g->north = 0.0;
g->is_three = is_three;
g->init = FALSE;
-#ifdef USE_GMT4
+
if(ggrd_grdtrack_init(&g->west,&g->east,&g->south,&g->north,&g->f,&g->mm,
grdfile,&g->grd,&g->edgeinfo,
gmt_edgeinfo_string,&geographic_in,
pad,is_three,depth_file,&g->z,&g->nz,
- bilinear,verbose,change_z_sign,&g->bcr))
+ bilinear,verbose,change_z_sign,
+ g->loc_bcr))
return 2;
-#else
- if(ggrd_grdtrack_init(&g->west,&g->east,&g->south,&g->north,&g->f,&g->mm,
- grdfile,&g->grd,&g->edgeinfo,
- gmt_edgeinfo_string,&geographic_in,
- pad,is_three,depth_file,&g->z,&g->nz,
- bilinear,verbose,change_z_sign))
- return 2;
-#endif
/*
check bandlimited maximum
@@ -110,7 +105,8 @@
//g->zlevels_are_negative = (g->zlevels_are_negative)?(FALSE):(TRUE);
if(verbose){
fprintf(stderr,"ggrd_grdtrack_init_general: initialized from %s, %s, bcflag: %s.\n",
- grdfile,(is_three)?("3-D"):("1-D"), gmt_edgeinfo_string);
+ grdfile,(is_three)?("3-D"):("1-D"),
+ gmt_edgeinfo_string);
if(is_three){
fprintf(stderr,"ggrd_grdtrack_init_general: depth file %s, %i levels, %s.\n",
depth_file,g->nz,
@@ -163,8 +159,8 @@
undefined and FALSE else
*/
ggrd_boolean ggrd_grdtrack_interpolate_rtp(double r,double t,double p,
- struct ggrd_gt *g,double *value,
- ggrd_boolean verbose)
+ struct ggrd_gt *g,double *value,
+ ggrd_boolean verbose)
{
double x[3];
ggrd_boolean result;
@@ -189,13 +185,10 @@
x[2] = (1.0-r) * 6371.0; /* depth in [km] */
if(g->zlevels_are_negative) /* adjust for depth */
x[2] = -x[2];
-#ifdef USE_GMT4
- result = ggrd_grdtrack_interpolate(x,TRUE,g->grd,g->f,g->edgeinfo,g->mm,g->z,g->nz,
- value,verbose,&g->bcr);
-#else
- result = ggrd_grdtrack_interpolate(x,TRUE,g->grd,g->f,g->edgeinfo,g->mm,g->z,g->nz,
- value,verbose);
-#endif
+ result = ggrd_grdtrack_interpolate(x,TRUE,g->grd,g->f,
+ g->edgeinfo,g->mm,g->z,
+ g->nz,value,verbose,
+ g->loc_bcr);
return result;
}
/*
@@ -205,9 +198,11 @@
undefined and FALSE else
this mean lon lat z
*/
-ggrd_boolean ggrd_grdtrack_interpolate_xyz(double x,double y,double z,
- struct ggrd_gt *g,double *value,
- ggrd_boolean verbose)
+ggrd_boolean ggrd_grdtrack_interpolate_xyz(double x,double y,
+ double z,
+ struct ggrd_gt *g,
+ double *value,
+ ggrd_boolean verbose)
{
double xloc[3];
ggrd_boolean result;
@@ -227,14 +222,10 @@
xloc[2] = z; /* depth, z*/
if(g->zlevels_are_negative) /* adjust for depth */
xloc[2] = -xloc[2];
-#ifdef USE_GMT4
- result = ggrd_grdtrack_interpolate(xloc,TRUE,g->grd,g->f,g->edgeinfo,g->mm,g->z,g->nz,
- value,verbose,&g->bcr);
-#else
- result = ggrd_grdtrack_interpolate(xloc,TRUE,g->grd,g->f,g->edgeinfo,g->mm,g->z,g->nz,
- value,verbose);
-
-#endif
+ result = ggrd_grdtrack_interpolate(xloc,TRUE,g->grd,g->f,
+ g->edgeinfo,g->mm,g->z,
+ g->nz,value,verbose,
+ g->loc_bcr);
return result;
}
@@ -272,13 +263,9 @@
x[0]-=360.0;
x[1] = 90.0 - t * ONEEIGHTYOVERPI; /* lat */
x[2] = 1.0;
-#ifdef USE_GMT4
- result = ggrd_grdtrack_interpolate(x,FALSE,g->grd,g->f,g->edgeinfo,g->mm,g->z,g->nz,
- value,verbose,&g->bcr);
-#else
- result = ggrd_grdtrack_interpolate(x,FALSE,g->grd,g->f,g->edgeinfo,g->mm,g->z,g->nz,
- value,verbose);
-#endif
+ result = ggrd_grdtrack_interpolate(x,FALSE,g->grd,g->f,
+ g->edgeinfo,g->mm,g->z,g->nz,
+ value,verbose,g->loc_bcr);
return result;
}
@@ -307,13 +294,11 @@
x[0] = xin;
x[1] = yin;
x[2] = 0.0;
-#ifdef USE_GMT4
- result = ggrd_grdtrack_interpolate(x,FALSE,g->grd,g->f,g->edgeinfo,g->mm,g->z,g->nz,
- value,verbose,&g->bcr);
-#else
- result = ggrd_grdtrack_interpolate(x,FALSE,g->grd,g->f,g->edgeinfo,g->mm,g->z,g->nz,
- value,verbose);
-#endif
+ result = ggrd_grdtrack_interpolate(x,FALSE,g->grd,g->f,
+ g->edgeinfo,
+ g->mm,g->z,g->nz,
+ value,verbose,
+ g->loc_bcr);
return result;
}
@@ -404,15 +389,16 @@
sign of the
depth
levels to go from depth (>0) to z (<0) */
- struct GMT_BCR *bcr)
+ struct GMT_BCR *loc_bcr)
#else
int ggrd_grdtrack_init(double *west, double *east,double *south, double *north,
float **f,int *mm,char *grdfile,struct GRD_HEADER **grd,
struct GMT_EDGEINFO **edgeinfo,char *edgeinfo_string,
ggrd_boolean *geographic_in,int *pad,ggrd_boolean three_d,
- char *dfile, float **z,int *nz,
+ char *dfile, float **z,int *nz,
ggrd_boolean bilinear,ggrd_boolean verbose,
- ggrd_boolean change_depth_sign)
+ ggrd_boolean change_depth_sign,
+ struct BCR *loc_bcr)
#endif
{
FILE *din;
@@ -428,8 +414,11 @@
/*
init first edgeinfo
*/
-
+#ifdef USE_GMT4
+ GMT_io_init (); /* Init the table i/o structure */
+#endif
GMT_boundcond_init (*edgeinfo);
+
/* check if geographic */
if (strlen(edgeinfo_string)>2){
GMT_boundcond_parse (*edgeinfo, (edgeinfo_string+2));
@@ -664,13 +653,13 @@
/*
Initialize bcr structure, this can be the same for
- all grids as long as they have the same dimesnions
+ all grids as long as they have the same dimensions
*/
#ifdef USE_GMT4
- GMT_bcr_init ((*grd+i), pad, bilinear,1.0,bcr);
+ GMT_bcr_init ((*grd+i), pad, bilinear,1.0,loc_bcr);
#else
- GMT_bcr_init ((*grd+i), pad, bilinear);
+ my_GMT_bcr_init ((*grd+i), pad, bilinear,loc_bcr);
#endif
}
/* Set boundary conditions */
@@ -709,7 +698,7 @@
int nz, /* number of depth layers */
double *value, /* output value */
ggrd_boolean verbose,
- struct GMT_BCR *bcr)
+ struct GMT_BCR *loc_bcr)
#else
ggrd_boolean ggrd_grdtrack_interpolate(double *in, /* lon/lat/z [2/3] in degrees/km */
ggrd_boolean three_d, /* use 3-D inetrpolation or 2-D? */
@@ -719,8 +708,9 @@
int mm, /* nx * ny */
float *z, /* depth layers */
int nz, /* number of depth layers */
- double *value, /* output value */
- ggrd_boolean verbose
+ double *value, /* output value */
+ ggrd_boolean verbose,
+ struct BCR *loc_bcr
)
#endif
{
@@ -770,12 +760,12 @@
*/
#ifdef USE_GMT4
- val1 = GMT_get_bcr_z((grd+i1), in[0], in[1], (f+i1*mm), (edgeinfo+i1),bcr);
- val2 = GMT_get_bcr_z((grd+i2), in[0], in[1], (f+i2*mm), (edgeinfo+i2),bcr);
+ val1 = GMT_get_bcr_z((grd+i1), in[0], in[1], (f+i1*mm), (edgeinfo+i1),loc_bcr);
+ val2 = GMT_get_bcr_z((grd+i2), in[0], in[1], (f+i2*mm), (edgeinfo+i2),loc_bcr);
#else
- ggrd_gt_bcr_init_loc ();
+ ggrd_global_bcr_assign(loc_bcr);
val1 = GMT_get_bcr_z((grd+i1), in[0], in[1], (f+i1*mm), (edgeinfo+i1));
- ggrd_gt_bcr_init_loc ();
+ ggrd_global_bcr_assign(loc_bcr);
val2 = GMT_get_bcr_z((grd+i2), in[0], in[1], (f+i2*mm), (edgeinfo+i2));
#endif
/* fprintf(stderr,"z(%3i/%3i): %11g z: %11g z(%3i/%3i): %11g f1: %11g f2: %11g v1: %11g v2: %11g rms: %11g %11g\n", */
@@ -786,8 +776,9 @@
}else{
/* single layer */
#ifdef USE_GMT4
- *value = GMT_get_bcr_z(grd, in[0], in[1], f, edgeinfo,bcr);
+ *value = GMT_get_bcr_z(grd, in[0], in[1], f, edgeinfo,loc_bcr);
#else
+ ggrd_global_bcr_assign(loc_bcr);
*value = GMT_get_bcr_z(grd, in[0], in[1], f, edgeinfo);
#endif
}
@@ -951,7 +942,8 @@
// not recalculate the weights
//
void ggrd_interpol_time(GGRD_CPREC time,struct ggrd_t *thist,
- int *i1,int *i2,GGRD_CPREC *f1,GGRD_CPREC *f2,
+ int *i1,int *i2,GGRD_CPREC *f1,
+ GGRD_CPREC *f2,
GGRD_CPREC dxlimit)
{
@@ -1131,7 +1123,8 @@
right++;
left = right - 1;
- f2 = (age - v->age_time[left])/(v->age_time[right]-v->age_time[left]);
+ f2 = (age - v->age_time[left])/
+ (v->age_time[right]-v->age_time[left]);
f1 = 1.0-f2;
//fprintf(stderr,"sai: %g %g %g\t%g %g\n",v->age_time[left],age,v->age_time[right],f1,f2);
v->sf_old_age = age;
@@ -1212,16 +1205,28 @@
#ifndef USE_GMT4
/*
-this is aweful, but works
+this is aweful, but works?
*/
-void ggrd_gt_bcr_init_loc(void)
+void ggrd_global_bcr_assign(struct BCR *loc_bcr)
{
- /* Initialize i,j so that they cannot look like they have been used: */
- bcr.i = -10;
- bcr.j = -10;
+ /* copy to global */
+ memcpy((void *)(&bcr),(void *)loc_bcr,sizeof(struct BCR));
}
+
+
+void my_GMT_bcr_init (struct GRD_HEADER *grd, int *pad,
+ int bilinear,struct BCR *loc_bcr)
+{
+
+ GMT_bcr_init(grd,pad,bilinear);
+ /* assign to local bcr */
+ memcpy((void *)loc_bcr,(void *)(&bcr),sizeof(struct BCR));
+}
+
+
+
#endif
/*
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.h
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.h 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.h 2007-09-27 22:43:58 UTC (rev 8044)
@@ -63,7 +63,10 @@
#else
/* GMT 3.4.5 */
ggrd_boolean ggrd_grdtrack_interpolate(double *, ggrd_boolean , struct GRD_HEADER *, float *,
- struct GMT_EDGEINFO *, int, float *, int , double *,ggrd_boolean);
+ struct GMT_EDGEINFO *, int,
+ float *, int ,
+ double *,ggrd_boolean,
+ struct BCR *);
int ggrd_grdtrack_init(double *, double *,double *, double *, /* geographic bounds,
set all to zero to
@@ -81,8 +84,14 @@
float **, /* layers, pass as NULL */
int *, /* number of layers */
ggrd_boolean , /* linear/cubic? */
- ggrd_boolean ,ggrd_boolean);
-void ggrd_gt_bcr_init_loc(void);
+ ggrd_boolean ,ggrd_boolean,
+ struct BCR *);
+
+void ggrd_global_bcr_assign(struct BCR *);
+
+void my_GMT_bcr_init (struct GRD_HEADER *, int *,
+ int ,struct BCR *);
+
#endif
/*
@@ -90,7 +99,6 @@
*/
-void ggrd_gt_bcr_init_loc(void);
void ggrd_gt_interpolate_z(double,float *,int ,
int *, int *, double *, double *,ggrd_boolean,
ggrd_boolean *); /* */
Modified: mc/1D/hc/trunk/ggrd_struc.h
===================================================================
--- mc/1D/hc/trunk/ggrd_struc.h 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/ggrd_struc.h 2007-09-27 22:43:58 UTC (rev 8044)
@@ -4,9 +4,8 @@
*/
-#ifdef USE_GMT4
#include "gmt.h"
-#endif
+
/*
plate tectonic stages interpolation structure
@@ -62,7 +61,10 @@
double west,east,south,north;
#ifdef USE_GMT4
- struct GMT_BCR bcr;
+#include "gmt_bcr.h"
+ struct GMT_BCR loc_bcr[1];
+#else
+ struct BCR loc_bcr[1];
#endif
};
Modified: mc/1D/hc/trunk/ggrd_test.c
===================================================================
--- mc/1D/hc/trunk/ggrd_test.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/ggrd_test.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -18,7 +18,7 @@
/*
read in velocities
*/
- if(ggrd_read_vel_grids(v,1.0,FALSE,TRUE,"/home/rcf-106/twb/data/plates/past/clb/hall/")){
+ if(ggrd_read_vel_grids(v,1.0,FALSE,TRUE,"/home/walter/becker/data/plates/past/clb/hall/")){
fprintf(stderr,"error opening grids\n");
exit(-1);
}
Added: mc/1D/hc/trunk/grdinttester.c
===================================================================
--- mc/1D/hc/trunk/grdinttester.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/grdinttester.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -0,0 +1,64 @@
+/*
+
+
+read in two grid files and compute the correlation of
+points within a certain radius around each location
+
+
+*/
+#include "ggrd_grdtrack_util.h"
+
+int main(int argc, char **argv)
+{
+ struct ggrd_gt g1[1],g2[1];
+ char char_dummy;
+ float x,y,x1,y1,dx,dy,rad,rad_km;
+ int n;
+ double *z1,*z2;
+ z1=(double *)malloc(sizeof(double));
+ z2=(double *)malloc(sizeof(double));
+
+ rad_km = 500;
+ if(argc < 3){
+ fprintf(stderr,"%s: usage:\n\t%s file1.grd file2.grd [radius, %g km]\n\n",
+ argv[0],argv[0],rad_km);
+ exit(-1);
+ }
+ if(argc > 3 )
+ sscanf(argv[3],"%f",&rad_km);
+
+ fprintf(stderr,"%s: assuming grids %s and %s are global, computing for radius %g km\n",
+ argv[0],argv[1],argv[2],rad_km);
+
+ if(ggrd_grdtrack_init_general(FALSE,argv[1],&char_dummy,
+ "-Lx",g1,TRUE,FALSE)){
+ fprintf(stderr,"%s: error reading %s\n",argv[0],argv[1]);
+ exit(-1);
+ }
+ if(ggrd_grdtrack_init_general(FALSE,argv[2],&char_dummy,
+ "-Lx",g2,TRUE,FALSE)){
+ fprintf(stderr,"%s: error reading %s\n",argv[0],argv[2]);
+ exit(-1);
+ }
+
+ /*
+
+ compute total correlation and best-fit slopes
+
+ */
+ dy = 2;
+ for(n=0,y1 = 90-dy/2+1e-7,y = -90+dy/2;y <= y1;y += dy){
+ dx = dy/cos(y/GGRD_PIF); /* adjust for sphericity */
+ for(x=0,x1=360-dx+1e-7;x<=x1;x+=dx,n++){
+ ggrd_grdtrack_interpolate_xy((double)x,(double)y,g1,(z1+n),FALSE);
+ ggrd_grdtrack_interpolate_xy((double)x,(double)y,g2,(z2+n),FALSE);
+
+ fprintf(stderr,"%g %g %g %g %.6e\n",x,y,z1[n],z2[n],fabs(z1[n]-z2[n]));
+ z1=(double *)realloc(z1,(n+2)*sizeof(double));
+ z2=(double *)realloc(z2,(n+2)*sizeof(double));
+ }
+ }
+
+
+
+}
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/hc.h 2007-09-27 22:43:58 UTC (rev 8044)
@@ -162,13 +162,19 @@
HEALPIX or RICK
*/
- int nrad; /* number of output radial layers */
+ int nrad; /* number of output radial
+ layers
+
+ this is determined by the density
+ model
+
+ */
HC_PREC *r; /* non-dim radii */
+ HC_PREC *dvisc; /* viscosity at the density layers */
/*
viscosity structure
*/
- int nvis; /* number of viscosity
- layers */
+ int nvis; /* number of viscosity layers */
HC_PREC *visc; /* values of normalized viscosities */
HC_PREC *rvisc; /* radii of normalized viscosities */
@@ -278,8 +284,10 @@
#define HC_VEL 0 /* compute velocities
v_r v_t v_p
*/
-#define HC_TRACTIONS 1 /* compute tractions
+#define HC_RTRACTIONS 1 /* compute tractions
trr trt trp */
+#define HC_HTRACTIONS 2 /* compute tractions ttt, ttp, tpp */
+
/*
init and assignment modes
@@ -309,6 +317,11 @@
#define HC_TPPROD 3
#define HC_NRNTNP 4
+
+/* for solution */
+#define HC_RAD 0
+#define HC_POL 1
+#define HC_TOR 2
/*
boolean
*/
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2007-09-27 22:43:58 UTC (rev 8044)
@@ -80,7 +80,7 @@
void hc_print_vector_label(double *, int, FILE *, char *);
void hc_print_matrix_label(double *, int, int, FILE *, char *);
void hc_print_vector_row(double *, int, FILE *);
-void hc_compute_solution_scaling_factors(struct hcs *, int, double, double *);
+void hc_compute_solution_scaling_factors(struct hcs *, int, double, double, double *);
void hc_print_poloidal_solution(struct sh_lms *, struct hcs *, int, char *, unsigned short, unsigned short);
void hc_print_toroidal_solution(double *, int, struct hcs *, int, char *, unsigned short);
/* hc_polsol.c */
@@ -93,7 +93,7 @@
void hc_sum(struct hcs *, int, struct sh_lms *, struct sh_lms *, int, unsigned short, struct sh_lms *, unsigned short);
void hc_compute_sol_spatial(struct hcs *, struct sh_lms *, float **, unsigned short);
/* hc_torsol.c */
-void hc_torsol(int, int, int, double *, double **, double **, struct sh_lms *, struct sh_lms *, double *, unsigned short);
+void hc_torsol(struct hcs *, int, int, int, double *, double **, double **, struct sh_lms *, struct sh_lms *, double *, unsigned short);
/* main.c */
/* prem_util.c */
int prem_find_layer_x(double, double, double *, int, int, double *);
Modified: mc/1D/hc/trunk/hc_extract_sh_layer.c
===================================================================
--- mc/1D/hc/trunk/hc_extract_sh_layer.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/hc_extract_sh_layer.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -12,7 +12,7 @@
int main(int argc, char **argv)
{
- int ilayer,nradp2,nsol,i,mode,shps,nset=1,loop,i1,i2;
+ int ilayer,nsol,i,mode,shps,nset=1,loop,i1,i2;
FILE *in;
struct sh_lms *sol=NULL;
struct hcs *model;
@@ -66,27 +66,26 @@
in = hc_open(argv[1],"r","hc_extract_sh_layer");
hc_read_sh_solution(model,&sol,in,binary,verbose);
fclose(in);
- nradp2 = model->nrad+2;
- nsol = nradp2 * 3;
+ nsol = model->nradp2 * 3;
/*
deal with selection
*/
loop = 0;
if(ilayer == -1)
- ilayer = nradp2;
+ ilayer = model->nradp2;
if(ilayer == -2){
- ilayer = nradp2;
+ ilayer = model->nradp2;
loop =1;
}
- if((ilayer<1)||(ilayer > nradp2)){
+ if((ilayer<1)||(ilayer > model->nradp2)){
fprintf(stderr,"%s: ilayer (%i) out of range, use 1 ... %i\n",
- argv[0],ilayer,nradp2);
+ argv[0],ilayer,model->nradp2);
exit(-1);
}
if(loop){
- i1=0;i2=nradp2-1;
+ i1=0;i2=model->nradp2-1;
if(short_format)
- printf("%i\n",nradp2);
+ printf("%i\n",model->nradp2);
}else{
i1=ilayer-1;i2 = i1;
}
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/hc_init.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -68,7 +68,7 @@
assign NULL pointers to allow reallocating
*/
(*hc)->r = (*hc)->visc = (*hc)->rvisc =
- (*hc)->dfact = (*hc)->rden = NULL;
+ (*hc)->dfact = (*hc)->rden = (*hc)->dvisc = NULL;
(*hc)->rpb = (*hc)->fpb= NULL;
(*hc)->dens_anom = NULL; /* expansions */
(*hc)->plm = NULL;
@@ -302,7 +302,8 @@
hc_name_boolean(p->print_pt_sol));
fprintf(stderr,"-px\t\tprint the spatial solution to file (%s)\n",
hc_name_boolean(p->print_spatial));
- fprintf(stderr,"-str\t\tcompute srr,srt,srp tractions [MPa] instead of velocities [cm/yr] (default: vel)\n");
+ fprintf(stderr,"-rtrac\t\tcompute srr,srt,srp tractions [MPa] instead of velocities [cm/yr] (default: vel)\n");
+ fprintf(stderr,"-htrac\t\tcompute stt,stp,spp tractions [MPa] instead of velocities [cm/yr] (default: vel)\n");
fprintf(stderr,"-v\t-vv\t-vvv: verbosity levels (%i)\n",
(int)(p->verbose));
fprintf(stderr,"\n\n");
@@ -340,8 +341,12 @@
}else if(strcmp(argv[i],"-visc")==0){ /* viscosity filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->visc_filename,argv[i],HC_CHAR_LENGTH);
- }else if(strcmp(argv[i],"-str")==0){ /* compute tractions */
- p->solution_mode = HC_TRACTIONS;
+ }else if(strcmp(argv[i],"-rtrac")==0){ /* compute radial
+ tractions */
+ p->solution_mode = HC_RTRACTIONS;
+ }else if(strcmp(argv[i],"-htrac")==0){ /* compute horizontal
+ tractions */
+ p->solution_mode = HC_HTRACTIONS;
}else if(strcmp(argv[i],"-v")==0){ /* verbosities */
p->verbose = 1;
}else if(strcmp(argv[i],"-vv")==0){ /* verbosities */
@@ -665,7 +670,13 @@
*/
hc->nrad = hc->inho;
- hc_vecrealloc(&hc->r,(2+hc->nrad),"hc_assign_density");
+ hc->nradp2 = hc->nrad + 2;
+ hc_vecrealloc(&hc->r,hc->nradp2,"hc_assign_density");
+ /*
+ viscosity at density layers for horizontal stress
+ computations */
+ hc_vecrealloc(&hc->dvisc,hc->nradp2,"hc_assign_density");
+
hc->r[0] = hc->r_cmb; /* CMB */
if(hc->rden[0] <= hc->r[0])
HC_ERROR("hc_assign_density","first density layer has to be above internal CMD limit");
@@ -676,6 +687,18 @@
hc->r[hc->nrad+1] = 1.0; /* surface */
/*
+ assign viscosity at density layers
+
+ */
+ hc->dvisc[0] = hc->visc[0];
+ for(i=1;i < hc->nradp2;i++){
+ for(j=hc->nvis-1;j>=0;j--)
+ if(hc->rvisc[j] < hc->r[i-1])
+ break;
+ hc->dvisc[i] = hc->visc[j];
+ }
+ /*
+
assign the density jump factors
*/
Modified: mc/1D/hc/trunk/hc_input.c
===================================================================
--- mc/1D/hc/trunk/hc_input.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/hc_input.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -39,6 +39,7 @@
for(i=0;i < nsol;i++) /* init expansions on irregular grid */
sh_init_expansion((*sol+i),lmax,hc->sh_type,1,verbose,FALSE);
hc->nrad = nset - 2;
+ hc->nradp2 = hc->nrad + 2;
hc_vecrealloc(&hc->r,nset,"hc_read_sh_solution");
}
/* assign depth */
Modified: mc/1D/hc/trunk/hc_output.c
===================================================================
--- mc/1D/hc/trunk/hc_output.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/hc_output.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -20,7 +20,7 @@
hc_boolean binary,
hc_boolean verbose)
{
- int nradp2,i,os;
+ int i,os;
static int ntype = 3; /* three sets of solutions, r/pol/tor */
HC_PREC fac[3];
if(!hc->spectral_solution_computed)
@@ -28,18 +28,17 @@
/*
number of solution sets of ntype solutions
*/
- nradp2 = hc->nrad+2;
- for(i=os=0;i < nradp2;i++,os += ntype){
+ for(i=os=0;i < hc->nradp2;i++,os += ntype){
/*
scale to cm/yr, or MPa for stress solutions
*/
- hc_compute_solution_scaling_factors(hc,sol_mode,hc->r[i],fac);
+ hc_compute_solution_scaling_factors(hc,sol_mode,hc->r[i],hc->dvisc[i],fac);
/*
write parameters, convert radius to depth in [km]
*/
- sh_print_parameters_to_file((sol+os),ntype,i,nradp2,
+ sh_print_parameters_to_file((sol+os),ntype,i,hc->nradp2,
HC_Z_DEPTH(hc->r[i]),
out,FALSE,binary,verbose);
/*
@@ -58,13 +57,20 @@
sqrt(sh_total_power((sol+os+2))),
fac[0]/11.1194926644559);
break;
- case HC_TRACTIONS:
- fprintf(stderr,"hc_print_spectral_solution: z: %8.3f str: |r|: %11.3e |pol|: %11.3e |tor|: %11.3e (scale: %g MPa)\n",
+ case HC_RTRACTIONS:
+ fprintf(stderr,"hc_print_spectral_solution: z: %8.3f rtrac: |r|: %11.3e |pol|: %11.3e |tor|: %11.3e (scale: %g MPa)\n",
HC_Z_DEPTH(hc->r[i]),sqrt(sh_total_power((sol+os))),
sqrt(sh_total_power((sol+os+1))),
sqrt(sh_total_power((sol+os+2))),
fac[0]/(0.553073278428428/hc->r[i]));
break;
+ case HC_HTRACTIONS:
+ fprintf(stderr,"hc_print_spectral_solution: z: %8.3f htrac: |r|: %11.3e |pol|: %11.3e |tor|: %11.3e (scale: %g MPa)\n",
+ HC_Z_DEPTH(hc->r[i]),sqrt(sh_total_power((sol+os))),
+ sqrt(sh_total_power((sol+os+1))),
+ sqrt(sh_total_power((sol+os+2))),
+ fac[0]/(0.553073278428428/hc->r[i]));
+ break;
default:
fprintf(stderr,"hc_print_spectral_solution: sol mode %i undefined\n",sol_mode);
exit(-1);
@@ -74,7 +80,7 @@
}
if(verbose)
fprintf(stderr,"hc_print_spectral_solution: wrote solution at %i levels\n",
- nradp2);
+ hc->nradp2);
}
/*
@@ -111,7 +117,7 @@
int sol_mode, hc_boolean binary,
hc_boolean verbose)
{
- int nradp2,i,j,k,os[3],los,np,np2,np3;
+ int i,j,k,os[3],los,np,np2,np3;
FILE *file_dummy=NULL,*out,*dout;
float flt_dummy=0,*xy=NULL,value[3];
HC_PREC fac[3];
@@ -119,7 +125,7 @@
if(!hc->spatial_solution_computed)
HC_ERROR("hc_print_spatial_solution","spectral solution not computed");
/* number of solution sets of ntype solutions */
- nradp2 = hc->nrad+2;
+
/* number of lateral points */
np = sol[0].npoints;
np2 = np*2;
@@ -137,14 +143,15 @@
if(verbose >= 2)
fprintf(stderr,"hc_print_spatial_solution: writing depth levels to %s\n",
dfilename);
- for(i=0;i < nradp2;i++){
+ for(i=0;i < hc->nradp2;i++){
/*
compute the scaling factors, those do depend on radius
in the case of the stresses, so leave inside loop!
*/
- hc_compute_solution_scaling_factors(hc,sol_mode,hc->r[i],fac);
+ hc_compute_solution_scaling_factors(hc,sol_mode,hc->r[i],
+ hc->dvisc[i],fac);
/* write depth in [km] to dout file */
fprintf(dout,"%g\n",HC_Z_DEPTH(hc->r[i]));
@@ -196,7 +203,7 @@
fclose(dout);
if(verbose)
fprintf(stderr,"hc_print_spatial_solution: wrote solution at %i levels\n",
- nradp2);
+ hc->nradp2);
free(xy);
}
@@ -208,10 +215,9 @@
void hc_print_depth_layers(struct hcs *hc, FILE *out,
hc_boolean verbose)
{
- int nradp2,i;
+ int i;
/* number of solution sets of ntype solutions */
- nradp2 = hc->nrad + 2;
- for(i=0;i < nradp2;i++)
+ for(i=0;i < hc->nradp2;i++)
fprintf(out,"%g\n",HC_Z_DEPTH(hc->r[i]));
}
@@ -288,6 +294,7 @@
void hc_compute_solution_scaling_factors(struct hcs *hc,
int sol_mode,
HC_PREC radius,
+ HC_PREC viscosity,
HC_PREC *fac)
{
@@ -295,9 +302,13 @@
case HC_VEL:
fac[0]=fac[1]=fac[2] = hc->vel_scale; /* go to cm/yr */
break;
- case HC_TRACTIONS:
+ case HC_RTRACTIONS: /* radial tractions */
fac[0]=fac[1]=fac[2] = hc->stress_scale/radius; /* go to MPa */
break;
+ case HC_HTRACTIONS: /* horizontal tractions, are actually
+ given as strain-rates */
+ fac[0]=fac[1]=fac[2] = 2.0*viscosity*hc->stress_scale/radius; /* go to MPa */
+ break;
default:
HC_ERROR("hc_print_spectral_solution","mode undefined");
break;
Modified: mc/1D/hc/trunk/hc_polsol.c
===================================================================
--- mc/1D/hc/trunk/hc_polsol.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/hc_polsol.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -190,7 +190,7 @@
// SUBROUTINE EVPPOT (L,RATIO,PPOT): OBTAINS PROPAGATOR FOR
// NON-EQUILIBRIUM POTENTIAL AND DERIVATIVE (RATIO IS R(I)/
// R(I+1), FOR PROPAGATION FROM R(I) TO R(I+1) AT L),
- int i,i2,i3,i6,j,l,m,nih,nxtv,ivis,os,pos1,pos2,nradp2,
+ int i,i2,i3,i6,j,l,m,nih,nxtv,ivis,os,pos1,pos2,
prop_s1,prop_s2,nvisp1,nzero,n6,ninho,nl=0,ip1;
int newprp,newpot,jpb,inho2,ibv,indx[3],a_or_b,ilayer,lmax,
nprops_max;
@@ -205,14 +205,22 @@
/*
define a few offset and size pointers
*/
- nradp2 = nrad + 2;
+
+#ifdef DEBUG
+ if(hc->nradp2 != nrad + 2){
+ fprintf(stderr,"hc_polsol: radius number mismatch\n");
+ exit(-1);
+ }
+#endif
+
+
inho2 = inho + 2;
nvisp1 = hc->nvis+1;
lmax = pol_sol[0].lmax ;
/*
max number of propagator levels, choose this generously
*/
- nprops_max = nradp2 * 3;
+ nprops_max = hc->nradp2 * 3;
/*
for prop and ppot: one set of propagators for all layers, there
lmax of those
@@ -223,8 +231,7 @@
/*
check if still same general number of layers
*/
- if((hc->psp.prop_params_init)&&((nradp2 != hc->nradp2)||
- (inho2 != hc->inho2)||
+ if((hc->psp.prop_params_init)&&((inho2 != hc->inho2)||
(nvisp1 != hc->nvisp1))){
HC_ERROR("hc_polsol","layer structure changed from last call");
}
@@ -323,7 +330,7 @@
/*
save in case we want to check if parameters changed later
*/
- hc->inho2 = inho2;hc->nvisp1=nvisp1;hc->nradp2=nradp2;
+ hc->inho2 = inho2;hc->nvisp1=nvisp1;
}
//
// SET RDEN(INHO+1) = 1.1 TO PREVENT TESTING OF THAT VALUE
@@ -339,7 +346,7 @@
hc->nprops = ivis = nih = 0;
hc->rprops[0] = rad[0];
hit = FALSE;
- for(i=1;(i < nradp2)&&(!hit);i++){
+ for(i=1;(i < hc->nradp2)&&(!hit);i++){
//
// INITIALIZE
//
Modified: mc/1D/hc/trunk/hc_solve.c
===================================================================
--- mc/1D/hc/trunk/hc_solve.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/hc_solve.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -123,7 +123,7 @@
hc_vecalloc(&tvec,(hc->nrad+2)*(hc->pvel[1].lmax+1)*2,
"hc_solve");
/* compute kernels, and assign kernel*pvel to tor_sol */
- hc_torsol(hc->nrad,hc->nvis,hc->pvel[1].lmax,hc->r,
+ hc_torsol(hc,hc->nrad,hc->nvis,hc->pvel[1].lmax,hc->r,
&hc->rvisc,&hc->visc,(hc->pvel+1),hc->tor_sol,tvec,
verbose);
if(print_pt_sol)
@@ -140,10 +140,14 @@
if(verbose)
fprintf(stderr,"hc_solve: computing solution for velocities\n");
break;
- case HC_TRACTIONS:
+ case HC_RTRACTIONS:
if(verbose)
- fprintf(stderr,"hc_solve: computing solution for tractions\n");
+ fprintf(stderr,"hc_solve: computing solution for radial tractions\n");
break;
+ case HC_HTRACTIONS:
+ if(verbose)
+ fprintf(stderr,"hc_solve: computing solution for horizontal tractions\n");
+ break;
default:
fprintf(stderr,"hc_solve: error: solution mode %i undefined\n",
solve_mode);
@@ -201,10 +205,10 @@
hc_boolean verbose)
{
int itchoose,irchoose,ipchoose; /* indices for which solutions to use */
- int i,j,i3,i6,nradp2;
+ int i,j,i3,i6;
if(sol[0].lmax > hc->lfac_init)
hc_init_l_factors(hc,sol[0].lmax);
- nradp2 = nrad + 2;
+
/*
pick the right components for the radial, poloidal, and toroidal
@@ -220,7 +224,7 @@
ipchoose = 1; // y2 for poloidal
itchoose = 0; // y9 for toroidal
break;
- case HC_TRACTIONS:
+ case HC_RTRACTIONS:
//
// srr srt srp stress output requested
//
@@ -228,6 +232,10 @@
ipchoose = 3;// y4 for poloidal
itchoose = 1;// y10 for toroidal
break;
+ case HC_HTRACTIONS:
+ fprintf(stderr,"hc_sum: horizontal tractions not implemented yet\n");
+ exit(-1);
+ break;
default:
HC_ERROR("hc_sum","solve mode undefined");
break;
@@ -240,24 +248,24 @@
*/
- for(i=i3=i6=0;i < nradp2;i++,i3+=3,i6+=6){
+ for(i=i3=i6=0;i < hc->nradp2;i++,i3+=3,i6+=6){
/*
radial part
*/
- sh_aexp_equals_bexp_coeff((sol+i3+0),(pol_sol+i6+irchoose));
+ sh_aexp_equals_bexp_coeff((sol+i3+HC_RAD),(pol_sol+i6+irchoose));
/*
poloidal part, need to scale with sqrt(l(l+1))
*/
- sh_aexp_equals_bexp_coeff((sol+i3+1),(pol_sol+i6+ipchoose));
- sh_scale_expansion_l_factor((sol+i3+1),hc->lfac);
+ sh_aexp_equals_bexp_coeff((sol+i3+HC_POL),(pol_sol+i6+ipchoose));
+ sh_scale_expansion_l_factor((sol+i3+HC_POL),hc->lfac);
for(j=0;j<3;j++)
sol[i3+j].spectral_init = TRUE;
if(!free_slip){
/*
toroidal part, need to scale with sqrt(l(l+1))
*/
- sh_aexp_equals_bexp_coeff((sol+i3+2),(tor_sol+i*2+itchoose));
- sh_scale_expansion_l_factor((sol+i3+2),hc->lfac);
+ sh_aexp_equals_bexp_coeff((sol+i3+HC_TOR),(tor_sol+i*2+itchoose));
+ sh_scale_expansion_l_factor((sol+i3+HC_TOR),hc->lfac);
}else{
/* no toroidal part for free-slip */
sh_clear_alm((sol+i3+2));
@@ -278,9 +286,8 @@
void hc_compute_sol_spatial(struct hcs *hc, struct sh_lms *sol_w,
float **sol_x, hc_boolean verbose)
{
- int i,i3,nradp2,np,np2,np3,os;
+ int i,i3,np,np2,np3,os;
static int ntype = 3;
- nradp2 = hc->nrad + 2;
np = sol_w[0].npoints;
np2 = np * 2;
np3 = np2 + np; /*
@@ -288,19 +295,19 @@
expansions for r, pol, tor
*/
/* allocate space for spatial solution*/
- hc_svecrealloc(sol_x,np3*nradp2,"sol_x");
+ hc_svecrealloc(sol_x,np3*hc->nradp2,"sol_x");
/*
compute the plm factors
*/
sh_compute_plm(sol_w,1,&hc->plm,verbose);
- for(i=i3=0;i < nradp2;i++,i3 += ntype){
+ for(i=i3=0;i < hc->nradp2;i++,i3 += ntype){
os = i*np3;
/* radial component */
- sh_compute_spatial((sol_w+i3),0,TRUE,&hc->plm,
+ sh_compute_spatial((sol_w+i3+HC_RAD),0,TRUE,&hc->plm,
(*sol_x+os),verbose);
os += np;
/* poloidal/toroidal component */
- sh_compute_spatial((sol_w+i3+1),1,TRUE,&hc->plm,
+ sh_compute_spatial((sol_w+i3+HC_POL),1,TRUE,&hc->plm,
(*sol_x+os),verbose);
}
hc->spatial_solution_computed = TRUE;
Modified: mc/1D/hc/trunk/hc_torsol.c
===================================================================
--- mc/1D/hc/trunk/hc_torsol.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/hc_torsol.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -48,7 +48,8 @@
//
//
//
-void hc_torsol(int nrad,int nvis,int lmax,HC_PREC *r,
+void hc_torsol(struct hcs *hc,
+ int nrad,int nvis,int lmax,HC_PREC *r,
HC_PREC **rv,HC_PREC **visc, struct sh_lms *pvel_tor,
struct sh_lms *tor_sol,double *tvec,
hc_boolean verbose)
@@ -62,7 +63,7 @@
//
double coef,*vecnor,hold,rlast,rnext,tloc[2],*tvec1,*tvec2;
double exp_fac[2],p[2][2],diflog,el,elp2,elm1,efdiff;
- int l,jvisp1,jvis,i,j,nvisp1,nradp2,lmaxp1,os;
+ int l,jvisp1,jvis,i,j,nvisp1,lmaxp1,os;
hc_boolean qvis;
//
// PASSED PARAMETERS: NRADP2: NUMBER OF OUTPUT RADII,
@@ -88,7 +89,7 @@
work!)
*/
nvisp1 = nvis + 1; /* length of rv and visc */
- nradp2 = nrad + 2; /* radius array */
+ //nradp2 = nrad + 2; /* radius array */
lmaxp1 = lmax+1; /* length of 0:lmax array */
/*
@@ -106,10 +107,14 @@
*/
HC_TVISC(nvis) = HC_TVISC(nvis-1); /* last entry in viscosity array */
#ifdef DEBUG
+ if(hc->nradp2 != nrad + 2){
+ fprintf(stderr,"hc_torsol: radius number mismatch\n");
+ exit(-1);
+ }
/*
test size of expansions
*/
- j = nradp2 * 2;
+ j = hc->nradp2 * 2;
for(i=0;i < j;i++){
if(tor_sol[i].lmax < pvel_tor->lmax){
fprintf(stderr,"hc_torsol: error: toroidal expansion %i has lmax %i, plates have %i\n",
@@ -131,7 +136,7 @@
/* solution factors as f(l,r) */
/* set local pointes */
tvec1 = tvec;
- tvec2 = (tvec + nradp2 * lmaxp1);
+ tvec2 = (tvec + hc->nradp2 * lmaxp1);
//
// (PREVENTS THE REQUESTING OF NON-EXISTANT VALUES)
//
@@ -174,7 +179,7 @@
tvec1[os] = tloc[0];
tvec2[os] = tloc[1];
- for(i=1;i < nradp2;i++){ /* loop through radii */
+ for(i=1;i < hc->nradp2;i++){ /* loop through radii */
os += lmaxp1;
//
// TEST FOR CHANGE IN VISCOSITY IN NEXT LAYER
@@ -224,12 +229,12 @@
// this
//
hc_dvecalloc(&vecnor,lmaxp1,"hc_torsol: vecnor");
- os = (nradp2-1) * lmaxp1;
+ os = (hc->nradp2-1) * lmaxp1;
vecnor[0] = 1.0;
for(l=1;l < lmaxp1;l++)
vecnor[l] = 1.0 / tvec1[os+l];
/* normalize */
- for(i=os=0;i < nradp2;i++,os+=lmaxp1)
+ for(i=os=0;i < hc->nradp2;i++,os+=lmaxp1)
for(l=0;l < lmaxp1;l++){
tvec1[os+l] *= vecnor[l];
tvec2[os+l] *= vecnor[l];
@@ -242,7 +247,7 @@
of l and depth
*/
- for(os=i=j=0;i < nradp2;i++,os+=lmaxp1,j+=2){
+ for(os=i=j=0;i < hc->nradp2;i++,os+=lmaxp1,j+=2){
/*
assign toroidal plate motion fields to solution expansion
*/
Modified: mc/1D/hc/trunk/main.c
===================================================================
--- mc/1D/hc/trunk/main.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/main.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -153,8 +153,10 @@
switch(p->solution_mode){
case HC_VEL:
sprintf(file_prefix,"vel");break;
- case HC_TRACTIONS:
- sprintf(file_prefix,"str");break;
+ case HC_RTRACTIONS:
+ sprintf(file_prefix,"rtrac");break;
+ case HC_HTRACTIONS:
+ sprintf(file_prefix,"htrac");break;
default:
HC_ERROR(argv[0],"solution mode undefined");break;
}
Modified: mc/1D/hc/trunk/make_tar
===================================================================
--- mc/1D/hc/trunk/make_tar 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/make_tar 2007-09-27 22:43:58 UTC (rev 8044)
@@ -1,22 +1,30 @@
#!/bin/bash
-ver=0.1
-date=`date '+%m-%d-%y'`
+#
+# create a tar file
+#
+ver=${1-0.1}
+date=`date '+%m%d%y'`
+tf=$HOME/tmp/hc-$ver.$date.tgz
+
rm -rf $HOME/tmp/hc
mkdir $HOME/tmp/hc
mkdir $HOME/tmp/hc/example_data
mkdir $HOME/tmp/hc/prem/
+mkdir $HOME/tmp/hc/hcplates/
cp *.c *.h Makefile README.TXT INSTALLATION.TXT calc_vel_and_plot \
$HOME/tmp/hc
cp example_data/* $HOME/tmp/hc/example_data/
-cp prem/prem* $HOME/tmp/hc/prem
+cp prem/prem* $HOME/tmp/hc/prem/
+cp hcplates/* $HOME/tmp/hc/hcplates/
cd $HOME/tmp/
-tar --create --gzip --verbose --file $HOME/tmp/hc-$ver.$date.tgz hc/
+tar --create --gzip --verbose --file $tf hc/
rm -r hc/
cd -
+echo $0: output in $tf
Modified: mc/1D/hc/trunk/rick_sh_c.c
===================================================================
--- mc/1D/hc/trunk/rick_sh_c.c 2007-09-27 22:41:54 UTC (rev 8043)
+++ mc/1D/hc/trunk/rick_sh_c.c 2007-09-27 22:43:58 UTC (rev 8044)
@@ -921,6 +921,7 @@
// first call, set up some factors. the arrays were allocated in rick_init
//
for(k=0,i=1;k < rick->nlon;k++,i++){
+ /* plm_srt[n] = sqrt(n+1) */
rick->plm_srt[k] = sqrt((SH_RICK_PREC)(i));
}
// initialize plm factors
@@ -961,9 +962,9 @@
mstop = l - 1;
for(m=1;m <= mstop;m++){
k++;
- rick->plm_fac1[k] = rick->plm_srt[l-m-1] * rick->plm_srt[l+m];
- rick->plm_fac2[k] = rick->plm_srt[l+m-1] * rick->plm_srt[l-m];
- if(m == 1){
+ rick->plm_fac1[k] = rick->plm_srt[l-m-1] * rick->plm_srt[l+m]; /* sqrt((l-m)(l+m+1) */
+ rick->plm_fac2[k] = rick->plm_srt[l+m-1] * rick->plm_srt[l-m]; /* sqrt((l+m)(l-m+1) */
+ if(m == 1){ /* multiply with sqrt(2) */
rick->plm_fac2[k] = rick->plm_fac2[k] * rick->plm_srt[1];
}
}
@@ -1054,6 +1055,7 @@
if(ivec){
//
// derivatives
+ //
// ---derivatives of P(z) wrt theta, where z=cos(theta)
//
dp[1] = -p[2];
@@ -1062,13 +1064,12 @@
for(l=2;l <= lmax;l++){
k++;
// treat m=0 and m=l separately
- dp[k] = -rick->plm_srt[l-1] * rick->plm_srt[l] / rick->plm_srt[1] * p[k+1];
- dp[k+l] = rick->plm_srt[l-1] / rick->plm_srt[1] * p[k+l-1];
+ dp[k] = -rick->plm_srt[l-1] * rick->plm_srt[l] / rick->plm_srt[1] * p[k+1]; /* m = 0 */
+ dp[k+l] = rick->plm_srt[l-1] / rick->plm_srt[1] * p[k+l-1]; /* m = l */
mstop = l-1;
- for(m=1;m <= mstop;m++){
+ for(m=1;m <= mstop;m++){ /* rest */
k++;
- dp[k] = rick->plm_fac2[k] * p[k-1] -
- rick->plm_fac1[k] * p[k+1];
+ dp[k] = rick->plm_fac2[k] * p[k-1] - rick->plm_fac1[k] * p[k+1];
dp[k] *= 0.5;
}
k++;
More information about the cig-commits
mailing list