[cig-commits] r22400 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Fri Jun 21 16:55:36 PDT 2013
Author: becker
Date: 2013-06-21 16:55:36 -0700 (Fri, 21 Jun 2013)
New Revision: 22400
Added:
mc/1D/hc/trunk/calc_dyn_topo
mc/1D/hc/trunk/dtopo_ref.ab
mc/1D/hc/trunk/geoid_ref.ab
mc/1D/hc/trunk/hc.c
mc/1D/hc/trunk/hc_invert_dtopo.c
mc/1D/hc/trunk/hc_visc_scan.c
mc/1D/hc/trunk/visc.dat
Removed:
mc/1D/hc/trunk/hc.tar
mc/1D/hc/trunk/hc_init.con
mc/1D/hc/trunk/main.c
Modified:
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/Makefile.include
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_constants.h
mc/1D/hc/trunk/hc_extract_sh_layer.c
mc/1D/hc/trunk/hc_extract_spatial.c
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/hc_input.c
mc/1D/hc/trunk/hc_misc.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/itg-hc-geoid.ab
mc/1D/hc/trunk/sh_ana.c
mc/1D/hc/trunk/sh_corr.c
mc/1D/hc/trunk/sh_exp.c
mc/1D/hc/trunk/sh_extract_layer.c
mc/1D/hc/trunk/sh_model.c
mc/1D/hc/trunk/sh_power.c
mc/1D/hc/trunk/sh_syn.c
Log:
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/Makefile 2013-06-21 23:55:36 UTC (rev 22400)
@@ -140,7 +140,8 @@
sh_tools: $(BDIR)/sh_syn $(BDIR)/sh_corr $(BDIR)/sh_ana $(BDIR)/sh_power \
$(BDIR)/sh_extract_layer
-hc_tools: $(BDIR)/hc $(BDIR)/hc_extract_sh_layer $(BDIR)/hc_extract_spatial \
+hc_tools: $(BDIR)/hc $(BDIR)/hc_visc_scan $(BDIR)/hc_invert_dtopo \
+ $(BDIR)/hc_extract_sh_layer $(BDIR)/hc_extract_spatial \
$(BDIR)/rotvec2vel $(BDIR)/print_gauss_lat
weird_tools: $(BDIR)/convert_bernhard_dens
@@ -204,13 +205,23 @@
$(CC) $(ODIR)/prem2dsm.o $(PREM_OBJS) -o $(BDIR)/prem2dsm -lm $(LDFLAGS)
-$(BDIR)/hc: $(LIBS) $(INCS) $(ODIR)/main.o $(PREM_OBJS)
- $(CC) $(LIB_FLAGS) $(ODIR)/main.o -o $(BDIR)/hc \
+$(BDIR)/hc: $(LIBS) $(INCS) $(ODIR)/hc.o $(PREM_OBJS)
+ $(CC) $(LIB_FLAGS) $(ODIR)/hc.o -o $(BDIR)/hc \
-lhc -lrick $(HEAL_LIBS_LINKLINE) $(PREM_OBJS) \
$(GGRD_LIBS_LINKLINE) -lm $(LDFLAGS)
-$(BDIR)/hc.dbg: $(LIBS) $(INCS) $(ODIR)/main.dbg.o $(PREM_OBJS)
- $(CC) $(LIB_FLAGS) $(ODIR)/main.dbg.o -o $(BDIR)/hc.dbg \
+$(BDIR)/hc_visc_scan: $(LIBS) $(INCS) $(ODIR)/hc_visc_scan.o $(PREM_OBJS)
+ $(CC) $(LIB_FLAGS) $(ODIR)/hc_visc_scan.o -o $(BDIR)/hc_visc_scan \
+ -lhc -lrick $(HEAL_LIBS_LINKLINE) $(PREM_OBJS) \
+ $(GGRD_LIBS_LINKLINE) -lm $(LDFLAGS)
+
+$(BDIR)/hc_invert_dtopo: $(LIBS) $(INCS) $(ODIR)/hc_invert_dtopo.o $(PREM_OBJS)
+ $(CC) $(LIB_FLAGS) $(ODIR)/hc_invert_dtopo.o -o $(BDIR)/hc_invert_dtopo \
+ -lhc -lrick $(HEAL_LIBS_LINKLINE) $(PREM_OBJS) \
+ $(GGRD_LIBS_LINKLINE) -lm $(LDFLAGS)
+
+$(BDIR)/hc.dbg: $(LIBS) $(INCS) $(ODIR)/hc.dbg.o $(PREM_OBJS)
+ $(CC) $(LIB_FLAGS) $(ODIR)/hc.dbg.o -o $(BDIR)/hc.dbg \
-lhc.dbg -lrick.dbg $(HEAL_LIBS_LINKLINE) $(PREM_OBJS) \
$(GGRD_LIBS_LINKLINE) -lm $(LDFLAGS)
Modified: mc/1D/hc/trunk/Makefile.include
===================================================================
--- mc/1D/hc/trunk/Makefile.include 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/Makefile.include 2013-06-21 23:55:36 UTC (rev 22400)
@@ -12,3 +12,4 @@
GGRD_INC_FLAGS = -I$(GMTHOME)/include -I$(NETCDFHOME)/include
GGRD_LIBS_LINKLINE = -lggrd -lgmt -lpsl -lnetcdf
+LDFLAGS = -lnetcdf
\ No newline at end of file
Added: mc/1D/hc/trunk/calc_dyn_topo
===================================================================
--- mc/1D/hc/trunk/calc_dyn_topo (rev 0)
+++ mc/1D/hc/trunk/calc_dyn_topo 2013-06-21 23:55:36 UTC (rev 22400)
@@ -0,0 +1,38 @@
+#!/bin/bash
+
+density=3377.97
+
+lmax=31
+
+
+#eo="$eo -dens $datadir/tomography/models/smean.$lmax.m.ab -dshs " # density model
+#eo="$eo -dnp"
+#eo="$eo -ds 0.2" # density scaling
+
+#eo="$eo -vf visc.dat " # viscosity
+#eo="$eo -cbckl 63" # solver trick
+
+rm rtrac.sol.bin geoid.ab
+
+hc -fs -rtrac $eo -v
+
+if [ ! -s rtrac.sol.bin ];then
+ echo $0: error, no output produced
+ exit
+
+fi
+
+hc_extract_sh_layer rtrac.sol.bin 2 4 | awk '{print($2)}' > sdepth.dat
+nlay=`lc sdepth.dat`
+
+# convert to m
+fac=`gawk -v rho=$density -v g=10 'BEGIN{print(-1/(rho*g)*1e6)}' `
+echo $nlay $fac
+
+
+hc_extract_sh_layer rtrac.sol.bin $nlay 1 0 > tmp.ab
+cat tmp.ab | \
+ sh_syn 0 0 0 360 -90 90 1 2> /dev/null | \
+ gawk -v fac=$fac '{print($1,$2,$3*fac)}' | \
+ xyz2grd -fg -Rg -I1 -Nnan -Gdtopo.grd
+grd2map dtopo.grd
Property changes on: mc/1D/hc/trunk/calc_dyn_topo
___________________________________________________________________
Added: svn:executable
+ *
Added: mc/1D/hc/trunk/dtopo_ref.ab
===================================================================
--- mc/1D/hc/trunk/dtopo_ref.ab (rev 0)
+++ mc/1D/hc/trunk/dtopo_ref.ab 2013-06-21 23:55:36 UTC (rev 22400)
@@ -0,0 +1,529 @@
+31 0 0 1 1 0
+ -1.4146777e+00 0.0000000e+00
+ -9.3923765e+02 0.0000000e+00
+ 3.4337773e+02 3.3885726e+02
+ -2.0916390e+02 0.0000000e+00
+ 4.4486779e+01 1.9547667e+02
+ -1.9990898e+00 3.4977607e+02
+ 1.2091933e+02 0.0000000e+00
+ -4.6918729e+00 1.0277326e+02
+ 2.5040465e+02 -4.3509441e+02
+ -2.1525507e+02 -1.8998632e+01
+ 2.2957090e+01 0.0000000e+00
+ -4.4357497e+01 7.8824059e+01
+ 2.9046278e+02 -3.1131340e+02
+ -1.2882295e+01 1.1740292e+02
+ 7.8207528e+01 -6.8961751e+02
+ 2.9691390e+02 0.0000000e+00
+ 9.2462863e+01 5.7470770e+01
+ 2.8860174e+02 -7.7037652e+01
+ -2.4650419e+02 3.5156377e+02
+ -6.3518407e+02 2.8516135e+02
+ 1.6741287e+02 1.3721536e+02
+ -3.6836650e+01 0.0000000e+00
+ -1.8812926e+02 -1.1643186e+02
+ 2.2012300e+02 -1.6307218e+02
+ -1.2200384e+02 3.1351889e+02
+ -4.4525125e+01 -1.2626583e+02
+ 1.3003690e+02 5.9695283e+00
+ -1.1867401e+02 2.6219670e+01
+ 4.9340226e+01 0.0000000e+00
+ 8.3324684e+00 1.5026884e+02
+ -1.0359203e+02 6.4761527e-01
+ -9.0335253e+01 1.9507903e+02
+ -6.7907083e+01 -8.9566067e+01
+ 9.9134708e+01 -3.7653929e+01
+ 1.0762806e+02 1.9098077e+02
+ 1.8066112e+01 -1.1321397e+02
+ 1.1893787e+02 0.0000000e+00
+ -7.9713004e+00 -1.7194117e+01
+ 2.7311794e+02 3.4349763e+01
+ 3.0739695e+01 -6.0474202e+01
+ -5.3773544e+00 -4.4525091e+01
+ -9.9919865e+01 -1.1893148e+02
+ 1.4145815e+02 -9.4057280e+01
+ -2.1504219e+01 9.8379379e+01
+ -6.3282068e+01 -1.4992536e+01
+ 5.5978673e+01 0.0000000e+00
+ -1.0058600e+02 -1.2615773e+02
+ 1.3423287e+02 -3.3046076e+01
+ 1.7628348e+02 1.4119629e+01
+ -3.3719773e+01 -9.0401877e+01
+ -3.4164803e+01 3.2271942e+01
+ 4.4627926e+01 2.9247631e+01
+ 1.8548288e+02 2.1357529e+01
+ -2.7829066e+01 6.3623054e+01
+ 5.9923544e+00 1.3794903e+02
+ 1.7065521e+02 0.0000000e+00
+ 4.4277142e+01 6.4964136e+01
+ 8.5889283e+01 -3.5586597e+01
+ -7.0284399e+01 -6.1279387e+01
+ -1.0407171e+02 3.1116815e+01
+ 8.3841608e+01 6.5894248e+01
+ -1.5120026e+02 -1.4511573e+02
+ 6.2857986e+01 -8.2029542e+01
+ 4.9004567e+01 3.2335530e+01
+ -1.1373545e+02 -1.3354153e+01
+ 1.0232285e+02 -1.3385828e+02
+ -1.6358785e+01 0.0000000e+00
+ 3.8087754e+01 -1.7713173e+02
+ 1.9012181e+01 -1.2423236e+02
+ 6.3457463e+01 8.1148954e+01
+ 5.6815283e+01 -4.7638285e+01
+ -2.6600229e+01 -3.9966799e+01
+ -1.2999735e+02 -3.7581893e+01
+ 1.3703544e+01 3.9295600e+01
+ -3.2715102e+01 8.7508390e+01
+ -1.5872161e+02 -1.4291772e+02
+ 1.8164509e+02 1.5866982e+02
+ -4.2755430e+01 -4.9408599e+01
+ 1.2278674e+02 0.0000000e+00
+ 4.3639738e+01 -7.8548898e+00
+ 6.2998561e+01 1.3948076e+02
+ -8.5893009e+01 -4.5095603e+01
+ 5.9722374e+01 -1.0883548e+02
+ 3.3356890e+01 1.0584684e+02
+ 1.4956741e+01 -4.7834492e+01
+ 8.4346090e+01 -1.9388808e+02
+ 4.4277645e+01 7.2348598e+01
+ 1.0057214e+02 4.5899418e+01
+ 8.9413477e+01 -5.7944608e+01
+ -5.9407441e+01 -1.3066683e+02
+ 2.9006785e+01 6.5724116e+01
+ 1.4707956e+02 0.0000000e+00
+ 2.4748917e+01 -9.0197929e+01
+ 4.5992252e+01 -2.1348114e+01
+ 9.1470042e+01 -8.8442562e+01
+ -3.3391209e+01 -5.3201859e+01
+ -4.4362557e+01 -5.8527642e+01
+ -1.0848703e+01 -1.8676098e+01
+ -9.5314190e+00 -3.6282922e+01
+ -2.6278418e+01 6.8591109e+01
+ -5.0416337e+00 -1.0773317e+02
+ 3.1953347e+01 -4.3829456e+01
+ 1.0840956e+02 4.8629976e+01
+ -8.3859807e+01 3.4046545e+00
+ -1.8936116e+01 -7.8891938e+01
+ 2.9725320e+01 0.0000000e+00
+ 2.1310448e+01 5.1916931e+01
+ 3.1157284e+00 5.0983528e+01
+ -9.6765388e+01 -1.1003603e+00
+ 3.2113406e+01 -6.2406912e+01
+ -6.5415220e+01 3.7575741e+01
+ -1.4759141e+02 5.0553992e+01
+ -1.9903839e+01 1.7391778e+01
+ -4.6981521e+01 3.8873433e+01
+ 1.7920460e+01 -5.3208975e+01
+ 2.8586248e+00 -3.3248346e+00
+ -5.3015405e+01 2.8206137e+01
+ -1.3958848e+02 3.4617718e+01
+ -2.4581963e+01 -1.2170622e+02
+ -1.6212277e+02 -2.4870432e+01
+ 4.6645727e+01 0.0000000e+00
+ -3.3260319e+01 -4.5584463e+01
+ -9.1762368e+01 -8.0074323e+01
+ -1.2514841e+02 -1.8671027e+01
+ -1.0000887e+02 -9.1338044e+01
+ 5.8099218e+01 2.0874544e+01
+ -1.1197311e+01 -1.9035437e+01
+ -1.7501702e+01 -7.1283508e+01
+ 1.6465270e+01 -3.7938517e+01
+ -9.9237637e-01 -4.4940090e+01
+ 3.5423968e+01 1.3710316e+01
+ -1.0697114e+02 -1.2486784e+02
+ 5.1531226e+01 7.4523244e+01
+ 6.5945003e+01 -5.7991774e+01
+ -6.0681388e+01 5.7073201e+01
+ 1.4461369e+02 2.7898416e+01
+ -3.4440601e+01 0.0000000e+00
+ -2.7220461e+01 -3.3206123e+01
+ -4.9651509e+01 4.2430798e+01
+ 1.0328937e+02 5.7626612e+01
+ 1.2492710e+02 6.8418102e+01
+ -2.2621133e+01 7.8246180e+01
+ -2.3062923e+00 -2.3178460e+01
+ 4.5776178e+01 7.1281267e+01
+ -1.9724381e+01 4.9193878e+01
+ 6.9062587e+00 4.6797416e+01
+ 1.2452067e+01 8.1727358e+01
+ -1.2687707e+02 -5.6476089e+01
+ 1.0287196e+01 4.6716421e+01
+ -2.4445677e+01 -2.9128640e+01
+ 4.6786752e+00 2.1452640e+01
+ -6.6306548e+01 5.2974628e+01
+ -4.0559157e+01 3.8019281e+01
+ 3.7819187e+01 0.0000000e+00
+ 6.6202259e+01 6.2770731e+01
+ -7.6702089e+01 -3.0103205e+00
+ -4.9215111e+01 3.8627546e+00
+ -3.7727078e+01 7.6585679e+01
+ -1.2967773e+01 9.9467500e-01
+ 6.8057384e+01 -4.1922341e+00
+ -7.1833303e+01 -1.3452505e+00
+ 1.2005191e+01 -7.4086062e+00
+ -7.9267829e+00 1.7010296e+01
+ -6.2629355e+01 6.9999164e+01
+ 7.3982390e+01 -4.6248216e+01
+ 6.0215828e+01 -8.7254276e+00
+ -4.1110356e+01 8.6916646e+01
+ -8.5090238e+00 5.4090375e+01
+ -1.2960539e+01 5.5420087e+01
+ 4.0508345e+01 -2.4653366e+01
+ 3.6514500e+01 -1.4831667e+01
+ -1.8101948e+01 0.0000000e+00
+ 5.5021002e+01 9.3405330e+01
+ 3.1898559e+01 -5.3294625e+01
+ 6.8216624e+01 1.1346464e+01
+ 6.0562518e+01 1.0903702e+01
+ -1.0268670e+02 -3.5844652e+01
+ 3.8198769e+01 8.3815384e+00
+ -2.6466976e+01 2.0125565e+01
+ 7.7190524e+01 -4.3455177e+01
+ 2.2729017e+01 7.2973388e+00
+ 1.0064869e+02 1.4859146e+01
+ -8.2231633e+00 -7.5779791e+00
+ -6.0315499e+01 -1.6185393e+01
+ 5.7425131e+01 2.6984024e+01
+ -9.2934592e+00 -5.6856278e+01
+ 5.6427007e+01 5.8619248e+00
+ 1.9750928e+01 -2.2273828e-01
+ -5.7509606e+01 9.3641718e+01
+ 9.2876420e-01 -6.1030208e+01
+ 8.3007350e+00 0.0000000e+00
+ 5.4286126e+01 -4.0452959e+01
+ 7.6237367e+01 -1.1223128e+01
+ -1.3824817e+01 5.5345907e+01
+ 5.2055395e+01 6.8063874e+01
+ 4.8714093e+01 -5.3327714e+01
+ -5.7550020e+01 4.1401530e+01
+ -7.1338327e+00 2.2233782e+01
+ 3.3147525e+01 -7.2096475e+01
+ -2.7999638e+01 -3.2896594e+01
+ -1.3018104e+01 7.5454435e+01
+ -8.4737978e-01 -3.3772539e+01
+ 4.1733915e+00 -2.8286274e+01
+ -1.2766708e+01 5.0539067e+01
+ -5.4354241e+01 -2.3574925e+00
+ 6.7824936e+01 -3.3280312e+01
+ -4.2956758e+00 -2.3754436e+01
+ 1.9216567e+01 2.8325118e+01
+ 6.5348928e+01 -9.6671105e+01
+ 1.2118859e+01 3.5897806e+01
+ 5.1213782e+01 0.0000000e+00
+ -7.1681396e+00 -6.4847350e+00
+ 1.8330805e+01 6.2463560e+01
+ 5.7362109e+00 -6.6392731e+01
+ 3.5713896e+01 -7.4217490e+01
+ -3.0282235e+01 2.4578286e+01
+ 2.4052388e+01 -1.8144182e+01
+ 1.2059094e+01 2.8637918e+01
+ -2.4408013e+01 -2.8423036e+01
+ -2.9716573e+01 -3.3340653e+01
+ 4.4031704e+01 -1.7707452e+01
+ -3.0400114e+01 4.6574897e+00
+ 2.1107294e+00 -8.4947288e+00
+ -4.1438604e+01 -3.1975509e+01
+ -3.9555917e+01 -3.0220022e+01
+ 1.6299683e+01 1.0710352e+01
+ -5.4933158e+00 5.5989506e+01
+ -1.1456286e+01 2.4180565e+01
+ 4.9939058e+01 5.5354422e+01
+ -2.8393513e+01 -2.9446818e+01
+ 1.8157862e+01 5.9503916e+01
+ 8.4505846e+01 0.0000000e+00
+ 2.2709631e+01 -3.1451552e+01
+ -3.3881068e+01 8.7131280e+01
+ -4.8312195e+00 7.2646044e+00
+ 4.5794057e+01 2.8847098e+01
+ 1.6591975e+00 4.8170408e+00
+ -7.9141435e+00 -3.6550509e+01
+ 5.1970429e+00 5.7001841e+01
+ -1.1027508e+01 -1.8464319e+01
+ -2.8468874e+01 -4.4998602e+01
+ 2.6469289e+01 8.0132511e+01
+ -7.3896024e+01 2.0681795e+01
+ -1.5692202e+01 2.6379713e+01
+ 1.0739479e+01 -5.5642214e+01
+ 5.4143141e+01 8.2276766e+00
+ -7.0328391e+01 4.5513122e+01
+ 6.4510027e+01 -1.0211944e+00
+ 9.2797308e+00 -2.6402397e+01
+ -3.2485541e+01 -8.8555196e+01
+ 7.0198260e+01 -1.5534128e+01
+ -1.7194155e+01 8.0395723e+01
+ -1.7489610e+01 -4.7756818e+01
+ -4.0821954e+01 0.0000000e+00
+ 3.5757193e+01 4.0296801e+01
+ -4.1095250e+01 3.0839164e+01
+ -3.9663004e+01 -2.6269079e+01
+ 7.6422066e+00 1.9677685e+01
+ -2.3286954e+01 -2.0033649e+01
+ 5.2229148e+01 -4.4164400e+01
+ -7.2243719e+01 -9.2712355e+00
+ -1.3212851e+01 -1.6013861e+01
+ 3.0646317e+00 -2.3921735e+01
+ 2.6157144e+01 6.3581861e+01
+ 3.7264611e+01 2.7353963e+01
+ -6.0158306e+01 -2.7880254e+01
+ -1.8144104e+01 -6.3623697e+01
+ 2.7370442e+01 3.4721852e+01
+ -7.4290519e+01 -2.7134280e+01
+ 6.4410794e+01 1.1997069e+01
+ 2.0233515e+01 -2.2409899e+00
+ 1.2160085e+01 -1.5195094e+01
+ 2.9445193e+01 2.6531731e+01
+ -1.1191808e+01 -2.2861418e+00
+ 6.4972289e+01 -4.5063450e+01
+ 5.6096448e+00 3.7656556e+01
+ -6.8202944e+01 0.0000000e+00
+ 9.8800638e+00 3.0932918e+01
+ -7.0388997e+01 -2.0304819e+01
+ 1.8451796e+01 5.9917052e+01
+ -6.0136561e+01 1.9591195e+01
+ -2.0472641e+01 -3.4379353e+01
+ -4.3145131e+01 2.7682809e+01
+ 1.7177457e+01 -8.9062213e+00
+ -1.2382546e+01 -1.0032579e+01
+ -1.6689468e+00 2.8623500e+01
+ 5.6591732e+01 -7.0491306e+00
+ 1.9840863e+01 1.0378577e+01
+ -1.3031395e+01 -1.3297048e+01
+ 2.0683963e+01 1.1928655e+01
+ 1.8898050e+01 2.9013807e+01
+ -2.3525635e+01 8.4098321e+00
+ 4.5883076e+01 -1.2697616e+01
+ 4.5280197e+01 6.9218448e+00
+ 1.0375083e+01 -1.3492909e+01
+ 1.9501073e+01 -2.2990850e+01
+ -2.2463793e-01 1.6013647e+00
+ 4.1225878e+01 -6.9276345e+01
+ -1.0762563e+02 -3.6686108e+01
+ 2.6886836e+01 3.8234077e+00
+ -8.4090306e+00 0.0000000e+00
+ 1.9590186e+01 1.7802099e+01
+ -1.0391804e+01 4.6622033e+01
+ 2.1414064e+01 1.7437170e+01
+ -1.8897553e+01 -3.4525934e+01
+ -1.3092789e+01 3.5264822e+01
+ 5.4716605e+01 2.2942391e+00
+ -1.0328640e+01 4.7268602e+01
+ 9.2786702e+01 -3.3078648e+01
+ 1.0749106e+01 2.4701371e+01
+ 2.2541202e+01 7.4164157e+01
+ 2.4528497e+00 -5.4940545e+01
+ 1.9234027e+01 -1.2630478e+01
+ 5.7350511e+01 -7.5498834e+00
+ -1.6988366e+01 9.8179864e+00
+ -1.8871878e+01 1.0875079e+01
+ -6.1491627e+00 -2.8636499e+01
+ 4.0908319e+01 -3.9334714e+01
+ -2.1236922e+01 -1.3722055e+01
+ 3.1704287e-01 -2.0328483e+00
+ 3.5984358e+01 -2.6121999e+01
+ -6.8318512e+01 1.4931524e+01
+ -6.5744734e+01 6.9248023e+01
+ 8.3297107e+01 -1.1103183e+01
+ -9.0696538e+00 -5.5475913e+00
+ -9.7594232e+00 0.0000000e+00
+ -5.1897769e+00 5.4780580e+01
+ -1.7812300e+01 -1.6816995e+00
+ 2.8115310e+01 -8.1279560e+00
+ 2.3537266e+01 -1.5092234e+01
+ 5.2979857e+01 3.0563523e+00
+ 5.5877532e+01 -7.2442905e+00
+ 2.5185766e+01 3.5416245e+00
+ 5.2386639e+00 4.3412631e+01
+ 4.9901652e+01 1.6741925e+01
+ 4.6622413e+01 2.0609948e+01
+ -1.5667334e+01 -1.4737037e+00
+ 2.1153161e+01 6.2100458e+01
+ -8.0225282e+00 -1.1239188e+01
+ -1.6553406e+01 2.7525883e+01
+ -3.7654055e+01 -1.0752777e+01
+ -5.4821550e+01 -4.0953219e+01
+ 5.9882407e+01 4.9345315e+01
+ -7.4696150e+01 2.6424253e+01
+ -5.2429990e+01 -5.2978129e+01
+ -1.7817541e+01 -1.0765344e+01
+ -4.7471054e+01 -6.0951046e+00
+ 3.3615345e+01 -1.1696034e+01
+ -4.6893316e+00 6.3407073e+00
+ -2.1394777e+01 -3.1354672e+01
+ -4.1912239e+01 -1.0212161e+01
+ 1.3991615e+01 0.0000000e+00
+ -3.1519543e+01 4.5953787e+01
+ 2.5901833e+01 3.4638551e+01
+ 1.1622710e+01 7.4469037e+00
+ 7.2923341e+01 1.4038816e+01
+ -1.9336391e+01 -9.9010173e+00
+ 3.2230605e+01 2.8813758e+01
+ -4.6981794e+01 7.7201416e+01
+ 2.9853811e+01 1.9150565e+01
+ 1.9004238e+01 1.8824528e+01
+ 1.8225677e+01 1.6092542e+01
+ -1.2739864e+01 -1.4563192e+01
+ -6.5015174e+01 1.7675419e+00
+ -1.9399023e+01 -2.3207053e-01
+ -7.8970190e+00 -7.8075538e+00
+ -8.3062052e+01 -2.2982266e+01
+ -1.4793211e+01 9.1795790e+00
+ 2.5105127e+01 1.1565839e+01
+ -5.0106625e+01 3.6395882e+01
+ 1.3626047e+01 -5.2048017e+01
+ -4.5421462e+01 -4.1286729e+01
+ -7.6586574e+00 2.4402085e+01
+ 2.4634729e+01 -2.0026926e+01
+ -5.1336799e+00 -6.2432902e+01
+ 3.1267223e+01 -2.6698478e+01
+ -2.7618256e+01 1.5905559e+01
+ 2.9322787e+01 3.8911336e+01
+ -2.4363495e+01 0.0000000e+00
+ -3.9514974e+01 2.2980328e+01
+ -1.4307020e+01 -5.5303316e+01
+ -6.1767379e+01 -3.2934723e+01
+ -9.2356659e+00 1.5869317e+01
+ -7.2088719e+00 -6.5594074e+01
+ 5.9946707e+01 2.0628049e+01
+ 3.0096951e+01 -3.4382650e+00
+ 2.9428720e+01 -3.1411272e+01
+ -2.7782152e+01 -8.1980990e+01
+ -3.4046817e+00 5.4651352e+00
+ 3.6711514e+01 -2.0672564e+01
+ -5.6751753e+01 3.0962583e+00
+ -2.2718332e+00 2.3064306e+01
+ -1.2812031e+00 -4.4426641e+00
+ -4.7123888e+01 1.5841012e+01
+ 4.3568062e+00 8.0209319e+00
+ 1.0747306e+01 -6.4019147e+00
+ 4.0433172e+01 -1.0179526e+01
+ -1.6010444e+01 1.6707502e+01
+ 3.7731998e+00 9.0297801e-01
+ 1.4111034e+01 2.6346447e+01
+ -2.5927625e+00 3.9680514e+01
+ -1.6143167e+01 -1.9195056e+01
+ -5.0931590e+01 -1.1325205e+01
+ -1.4586306e+00 9.4937438e+01
+ -1.0199240e+01 -2.6123409e+01
+ -3.2287453e+01 2.6493503e+01
+ -6.4931130e+00 0.0000000e+00
+ -1.8460750e+01 -4.9396058e+01
+ -5.9137957e-01 -1.5112294e+01
+ 7.1485951e+00 -3.7195080e+01
+ 2.5823671e+01 3.6045806e+01
+ -4.9534999e+00 -1.9933241e+01
+ 2.4522613e+01 7.4905249e-01
+ 1.5577903e+01 -4.2176794e+01
+ 5.7484115e+01 -3.1521196e+01
+ -2.8612616e+01 -2.0423934e+01
+ -4.5909705e+01 -1.1546959e+01
+ 5.4399777e+01 6.5456005e+00
+ 2.1274828e+01 3.6288556e+01
+ 1.4374500e-01 -8.5800664e+00
+ -2.9692051e+01 -2.2914103e+01
+ 7.5477733e+01 2.1915305e+01
+ -1.2091463e+01 1.0496121e+00
+ 6.3532432e+00 2.4800170e+01
+ 4.0826931e+01 -7.6589933e+00
+ -3.7007927e+01 -3.4640714e+01
+ 1.5625529e+01 1.5524385e+01
+ 5.7154250e+00 1.1177867e+01
+ 4.9067740e+00 7.2748196e+01
+ -6.0429002e+00 8.9142904e+00
+ 1.5259382e+01 1.6655248e+01
+ -4.8284509e+00 -2.7661888e-01
+ 4.3410545e+01 5.3308177e+00
+ -5.8590442e-01 -1.7549403e+01
+ 2.3981500e+01 1.3399154e+01
+ -1.3987577e+01 0.0000000e+00
+ -4.3303033e+00 -7.9708505e+01
+ 6.8856779e+00 -9.9873622e+00
+ 1.7750284e+01 -6.4480815e+01
+ -7.1983747e+01 -5.4615034e+00
+ 3.6366416e+01 -4.9221173e+01
+ -6.6195129e+00 -2.7595268e+01
+ 3.7256703e+01 5.8395431e+00
+ -4.4307833e+01 -4.5228053e+01
+ -1.4688500e+01 2.0546274e+01
+ 7.3563936e+01 1.7462069e+01
+ 1.1615264e+01 -3.1820408e+00
+ 9.5171740e+00 3.4178168e+00
+ -7.5458571e+00 1.0319494e+01
+ -1.5909717e+01 -3.3280155e+00
+ 6.3541176e+00 3.3380907e+01
+ 2.4451728e+01 2.2454565e+01
+ 1.6445748e+01 1.9600643e+00
+ -3.1219959e+01 2.9474588e+01
+ -3.1670568e+01 -1.7364201e+01
+ -5.8923204e+01 6.1033746e+00
+ 1.9023130e+01 6.9452766e+00
+ 3.7330066e+00 -1.4753482e+01
+ -3.2668976e+01 2.9399502e+01
+ 3.2250878e+01 -1.8002414e+01
+ -2.6320271e+01 -2.8821113e+00
+ 3.3165539e+01 1.2750613e+00
+ 6.7400235e+00 -1.9909086e+01
+ 1.4190740e+01 -2.3422263e+01
+ -1.0186071e+01 -3.3257413e+01
+ -7.0928036e-01 0.0000000e+00
+ -2.5657456e+01 -3.5950237e+01
+ -1.9747483e+01 -1.2467827e+01
+ -1.8311908e+01 -2.9260259e+01
+ 1.5906124e+01 -1.4262033e+01
+ -1.2006994e+01 8.1589998e+00
+ 2.8454213e+00 3.4001604e+01
+ -5.0224740e+01 2.0806849e+01
+ 3.4433848e+01 4.7775774e+00
+ -6.4992455e+00 2.6500780e+01
+ 2.9650382e+00 1.5689118e+00
+ 3.0667000e+01 2.2959308e+01
+ -5.3376878e+00 2.9613091e+01
+ -6.4644347e+01 -5.8867484e+00
+ 1.5610753e+00 2.5145409e+01
+ -2.2551589e+01 -1.4103767e+01
+ -1.7447960e+00 -2.5510234e+01
+ 4.8626095e+00 -3.3638567e+01
+ 2.5282414e+01 -1.2179362e+01
+ 1.1852109e+01 -5.0866010e-01
+ -3.5477568e+01 3.0300867e+01
+ 3.2114010e+01 -3.6207815e+01
+ 3.7396053e+01 3.6165693e+00
+ 5.2585866e+01 -9.3292963e+00
+ 3.8514868e+01 -4.2781571e+01
+ -2.1688373e+00 2.9977468e+01
+ -7.5446410e+00 8.3778051e+00
+ 4.8219587e+00 -1.8010763e+01
+ 2.5464649e+01 2.4784509e+01
+ -4.7739706e+01 -2.7107885e+01
+ -1.4799859e+01 5.2765782e+01
+ -1.0636076e+01 0.0000000e+00
+ 3.0668796e+01 1.9299831e+01
+ -1.1132955e+01 2.3288033e+01
+ -2.4078840e+01 1.6407655e+01
+ 4.5125007e+00 -1.5750263e+00
+ -6.5664139e+00 1.0139465e+01
+ -6.5565026e+00 1.2648138e+01
+ -2.8487205e+01 1.7840685e+00
+ -1.7420635e+01 2.4041437e+01
+ -1.2448712e+01 2.8054257e+00
+ -1.8755394e+01 3.4621914e+00
+ 8.3314760e+00 -4.3952759e+01
+ 1.1341652e+00 1.0589179e+01
+ -3.0046551e+01 -1.2255459e+01
+ 3.2717119e+00 -1.4792226e+01
+ -3.4847347e+01 8.6426779e-01
+ 1.5554738e+00 2.5538360e+01
+ -3.1340086e+01 -3.1303788e+01
+ 1.0539030e+01 -5.7391529e+00
+ -2.5569941e+01 -1.2790877e+01
+ -9.5466667e+00 -8.5172875e+00
+ 2.4449944e+01 -1.0842394e+01
+ -3.6580371e+00 1.1394211e+01
+ -6.6129329e+00 -1.1802587e+00
+ -2.9181253e+01 -3.8194032e+01
+ 4.9433036e+01 3.3296006e+01
+ -1.2256459e+01 -1.3520069e+01
+ -2.8510369e+01 1.0547512e+01
+ -5.1983899e+00 -8.8298359e+01
+ -9.2447936e-01 8.6946379e-01
+ 2.1149474e+01 3.8427102e+01
+ -3.4419920e+00 -1.4605041e+01
Added: mc/1D/hc/trunk/geoid_ref.ab
===================================================================
--- mc/1D/hc/trunk/geoid_ref.ab (rev 0)
+++ mc/1D/hc/trunk/geoid_ref.ab 2013-06-21 23:55:36 UTC (rev 22400)
@@ -0,0 +1,529 @@
+31 0 0 1 1 0
+ 0 0
+ 0 0
+ 0 0
+ -1.1552017e+02 0.0000000e+00
+ 6.0024606e-03 -3.3358521e-02
+ 5.5154257e+01 -3.1660072e+01
+ 2.1641332e+01 0.0000000e+00
+ -4.5908573e+01 -5.6117892e+00
+ 2.0457177e+01 -1.3995660e+01
+ -1.6309025e+01 -3.1978316e+01
+ -1.0084792e+01 0.0000000e+00
+ 1.2122473e+01 1.0707316e+01
+ 7.9248115e+00 1.4978616e+01
+ -2.2403185e+01 4.5436135e+00
+ -4.2624132e+00 6.9820287e+00
+ 1.5526296e+00 0.0000000e+00
+ 1.4226425e+00 2.1336924e+00
+ 1.4743427e+01 -7.3109878e+00
+ 1.0216225e+01 4.8601233e+00
+ -6.6773573e+00 1.1261335e+00
+ -3.9524793e+00 1.5134622e+01
+ -3.3904447e+00 0.0000000e+00
+ 1.7165669e+00 -5.9943994e-01
+ 1.0999473e+00 -8.4513437e+00
+ -1.2943092e+00 -2.0240396e-01
+ -1.9449907e+00 -1.0658889e+01
+ 6.0406091e+00 1.2130062e+01
+ 2.1413135e-01 -5.3671935e+00
+ 2.0464702e+00 0.0000000e+00
+ -6.3508423e+00 -2.1507888e+00
+ 7.4704952e+00 2.1026523e+00
+ -5.6628423e+00 4.9090253e+00
+ -6.2175885e+00 -2.8049488e+00
+ -3.7254661e-02 -4.0535459e-01
+ -8.1124005e+00 3.4321448e+00
+ -3.4083299e-02 -5.4505491e-01
+ 1.1186388e+00 0.0000000e+00
+ -5.2366355e-01 -1.3316660e+00
+ 1.8091174e+00 1.4759864e+00
+ 4.3805692e-01 1.9436370e+00
+ -5.5249689e+00 1.5783359e+00
+ 5.8110093e-01 -2.0168827e+00
+ -1.4914593e+00 6.9852572e+00
+ -1.5206739e+00 -1.6927716e+00
+ -2.8041431e+00 2.7256666e+00
+ 6.3348611e-01 0.0000000e+00
+ -3.2140295e+00 -4.8386258e-01
+ 4.8417848e-01 -7.1669830e-01
+ 3.6314306e+00 1.6791448e+00
+ -2.1174864e-01 4.4999766e-01
+ 3.6884447e-01 1.2218277e+00
+ 1.4196297e+00 5.0411585e+00
+ 2.6676050e+00 2.1914030e+00
+ 4.2537420e+00 -6.7951880e-02
+ 1.0752552e+00 -2.1904579e+00
+ 1.2057962e+00 0.0000000e+00
+ -1.8938581e+00 2.9639854e+00
+ -2.1250933e+00 -1.1593157e+00
+ 1.5842986e-01 3.4850895e+00
+ -1.9098938e+00 -1.7867603e+00
+ 1.1144287e+00 1.1443718e+00
+ -8.4979112e-01 -1.8035664e+00
+ -1.8680534e-01 6.8938096e-02
+ 9.1792058e-01 -2.0736418e+00
+ -2.8347539e+00 8.5790235e-01
+ 2.2708476e+00 -5.3946421e-01
+ -1.1478683e+00 0.0000000e+00
+ -3.5300326e-01 6.1326054e-01
+ 4.5476489e-01 -2.2383896e+00
+ 6.9135113e-01 3.3651531e+00
+ -8.5804374e-01 -1.4417656e+00
+ -8.4604554e-01 -1.1212439e+00
+ -3.5368299e-02 7.7492096e-01
+ -1.0524050e-01 2.0309399e+00
+ -1.4248203e-01 5.5495253e-01
+ 7.0255311e-01 -9.5116010e-01
+ -1.1812434e+00 -4.1651145e-01
+ -1.0453473e+00 1.5752579e+00
+ 8.2381883e-01 0.0000000e+00
+ 1.2115660e+00 9.7597022e-01
+ 3.2256633e-01 7.0302605e-01
+ -8.9582990e-01 -5.6665534e-01
+ -1.5313341e+00 8.6782206e-02
+ -6.9813829e-01 -1.7162411e-01
+ 7.0864204e-02 8.8133841e-01
+ 4.3075966e-01 -8.0778107e-01
+ -5.8529564e-01 3.8292711e-01
+ -9.4768857e-01 -5.6440114e-01
+ -1.4017131e-01 6.9954689e-01
+ -2.5695060e-01 1.4437696e-01
+ -5.4800504e-02 -2.5095501e-01
+ 9.4349552e-01 0.0000000e+00
+ 1.1631010e+00 -8.7480084e-01
+ 1.2505957e+00 -1.4175132e+00
+ 4.8740249e-01 -2.2086861e+00
+ -8.2555077e-02 -2.6569509e-01
+ -1.3197453e+00 -1.5199387e+00
+ -7.9235390e-01 -1.4184480e-01
+ -6.8148962e-02 1.6552011e-01
+ -2.2730250e-01 -2.2288545e-01
+ -5.6005351e-01 -1.0375134e+00
+ 9.2944808e-01 -8.3295758e-01
+ 1.0066225e+00 1.0946479e-01
+ -7.0798537e-01 1.9882619e+00
+ 1.3837363e+00 -1.5408683e+00
+ -5.1252387e-01 0.0000000e+00
+ 4.2444413e-01 -6.5252642e-01
+ -8.1211768e-01 -9.1644026e-02
+ -8.2557902e-01 -4.4528378e-01
+ 3.6217430e-02 -5.1239724e-01
+ -6.6267961e-01 3.7960724e-01
+ -4.3111297e-01 5.5543702e-02
+ -8.5080500e-01 8.8939340e-02
+ -7.9002966e-01 -3.4926779e-01
+ -7.2242653e-01 -6.4357335e-01
+ 8.7728291e-01 -2.9246385e-02
+ -3.5378974e-01 8.8269882e-01
+ 1.9135094e-01 -7.0364525e-01
+ -7.2902589e-01 -1.0207761e+00
+ -1.1726633e+00 -1.0889321e-01
+ 4.9564537e-02 0.0000000e+00
+ -2.1319845e-01 -2.3702698e-01
+ -4.6418805e-01 -6.8509660e-01
+ -1.2077371e+00 -3.9936823e-01
+ -9.0828893e-01 1.5405099e-01
+ -2.7680073e-01 -1.7230440e-01
+ 7.4278591e-01 -8.2456307e-01
+ -1.3487732e+00 -1.1472124e-01
+ -7.2549525e-01 5.0132781e-01
+ -3.0071159e-01 -8.5898369e-01
+ 2.3206599e-01 3.3215392e-01
+ 2.9497681e-02 -4.1875925e-01
+ -7.3289324e-01 3.5289650e-01
+ 6.4133104e-01 1.0345233e-01
+ 1.1754122e-01 -5.5156973e-01
+ 4.3059130e-01 1.0625352e-01
+ -1.0650115e-01 0.0000000e+00
+ -5.9204579e-01 -7.5386766e-01
+ -5.5420903e-01 6.3378904e-01
+ 7.6680889e-01 4.8249775e-01
+ 9.2370533e-01 1.0849982e+00
+ 2.7405309e-01 7.7842335e-02
+ 3.1370625e-01 -8.0484238e-01
+ 1.8227859e-01 1.9561596e-01
+ -4.7942974e-01 1.2225567e-01
+ 5.0680541e-01 8.9690427e-01
+ -2.6694299e-01 2.6086036e-01
+ -4.3211705e-01 7.2363896e-02
+ 4.4232154e-01 1.5206076e-01
+ -3.1143752e-01 -2.3682779e-02
+ -4.3739247e-01 -8.7392344e-01
+ 3.2602454e-01 7.4106262e-01
+ -8.6594307e-01 6.6895512e-02
+ 4.3382967e-01 0.0000000e+00
+ 5.7348925e-01 7.1682983e-01
+ -4.5449842e-01 1.5405844e-01
+ -1.4263221e-01 -1.1492811e-01
+ 1.4640258e-01 5.7283233e-01
+ 3.6671902e-01 -1.8155800e-01
+ -2.6526012e-01 -6.6563835e-01
+ -5.6471519e-01 9.9156565e-02
+ 8.8180097e-01 8.2173451e-02
+ -7.8660140e-02 6.2499849e-01
+ -8.5824750e-02 4.1552066e-01
+ 3.6215416e-01 -2.5067718e-01
+ 6.4940586e-01 4.6255842e-01
+ -3.7320606e-01 -4.5526500e-01
+ -3.2275133e-01 2.6169192e-01
+ -1.2520965e-01 -1.1855238e-01
+ -6.8732396e-01 8.1457776e-02
+ 7.8449753e-01 4.4935566e-01
+ 1.3788946e-01 0.0000000e+00
+ -1.6282574e-01 8.8850136e-01
+ 3.3293400e-01 2.4494189e-01
+ 1.1405743e-01 1.3225059e-01
+ 1.2350736e+00 -1.7788308e-02
+ -1.3512834e-01 -5.9101961e-01
+ 3.0698756e-01 -2.9945936e-01
+ -1.5359958e-01 -1.6884867e-01
+ 6.8958528e-01 9.8226325e-02
+ 4.4246866e-01 -8.1691095e-01
+ 1.1791713e-01 -9.5881496e-02
+ 1.5570110e-01 -4.7916300e-02
+ -6.7262171e-01 -3.7444121e-01
+ 1.4139723e-01 7.8997124e-01
+ -1.8754861e-01 -2.9015456e-01
+ 9.1500006e-01 4.5854365e-01
+ 2.2983682e-01 1.4705991e-01
+ -7.8755718e-02 -9.8950184e-02
+ 6.7619320e-02 -2.4552017e-01
+ -7.4683590e-02 0.0000000e+00
+ 2.0284654e-01 -2.6985302e-02
+ 8.0805137e-01 -5.3589083e-02
+ 1.7093820e-01 -2.4936495e-02
+ 3.5742944e-01 -1.8410757e-01
+ -2.3489247e-01 -6.2065387e-01
+ -1.0851947e-01 4.2407984e-01
+ -1.2755008e-01 1.9765385e-01
+ 6.7578987e-01 -2.2592181e-01
+ -7.4343793e-02 -1.6364742e-01
+ -7.6599096e-01 -1.7136150e-01
+ -3.6864520e-01 -2.3417023e-01
+ -5.5163696e-02 2.1380449e-01
+ 1.7050757e-01 6.4491929e-01
+ -1.0725415e-01 -2.9054800e-01
+ 3.9770151e-01 3.1815696e-01
+ -4.9165494e-01 -1.6047799e-01
+ -6.5126577e-01 3.4689572e-01
+ 7.9154864e-01 -2.1908425e-01
+ 6.1142012e-02 -1.1751521e-01
+ 4.8745010e-01 0.0000000e+00
+ -1.2588497e-01 -1.5892615e-01
+ 4.5879685e-01 3.8854289e-01
+ 1.0741226e-01 -8.7952349e-01
+ 9.6545080e-02 -5.1156569e-01
+ 2.2894158e-01 1.8675606e-01
+ 2.7498255e-01 -9.8670087e-02
+ 4.7992468e-01 1.5947336e-02
+ 1.1520104e-01 4.8696888e-02
+ -3.8859706e-01 1.5855756e-01
+ -7.3020511e-01 -1.0829557e-01
+ -3.2671146e-01 4.3300745e-01
+ -1.4622963e-01 4.1002844e-01
+ -6.1979249e-01 -1.5292980e-01
+ 2.6060780e-01 -3.2517624e-01
+ 5.8312466e-01 1.9558166e-02
+ -2.8091050e-01 -7.7652063e-03
+ -1.0182331e-01 3.1048587e-01
+ 3.4721398e-01 -2.0291260e-02
+ 6.8514252e-02 -2.4695664e-01
+ 8.4444088e-02 -2.8703100e-01
+ 1.4130671e-01 0.0000000e+00
+ 3.6534193e-01 -6.4872874e-01
+ -1.2699291e-01 9.4485779e-02
+ -4.3288999e-01 -5.1377364e-01
+ -1.1553034e-01 4.2924595e-01
+ -5.8926312e-02 -2.4806854e-02
+ -3.0635105e-01 3.1005953e-04
+ 1.4946147e-01 -5.7004250e-02
+ -3.8959731e-01 4.8583959e-02
+ -3.5851115e-01 -2.0526216e-01
+ -2.5774692e-01 -2.5785994e-02
+ -1.5587827e-01 8.0442622e-01
+ -7.5206089e-02 3.3405902e-01
+ 4.3780698e-01 -3.1513684e-01
+ 4.5818841e-01 1.6627079e-01
+ -3.9732611e-01 -2.3444144e-01
+ 1.6940902e-01 -1.5043448e-01
+ 1.6882493e-01 1.6680517e-01
+ 5.9850159e-01 -2.5033960e-01
+ 6.1504979e-01 -3.7431128e-01
+ -6.0932811e-01 3.5992488e-01
+ -1.9125652e-01 8.2596874e-02
+ -2.4367263e-01 0.0000000e+00
+ -3.5454617e-01 8.6835663e-02
+ -5.9892003e-01 -2.6332819e-02
+ -2.5648420e-01 -2.2884447e-01
+ -8.1734864e-02 4.2342631e-01
+ -2.1888649e-03 7.2012142e-03
+ 2.2982526e-01 -1.2452710e-01
+ -3.9487150e-01 -1.0733988e-01
+ -5.3130574e-01 8.6246782e-02
+ -1.5109905e-01 -1.8939529e-01
+ 1.2723412e-01 5.0536216e-01
+ 1.0887340e-01 3.9711132e-01
+ 5.4680327e-02 -1.8558948e-01
+ 3.8796623e-01 -4.4118132e-01
+ 2.4577967e-01 1.8826578e-01
+ -5.8257648e-01 -1.0652508e-01
+ 2.5508395e-03 -1.6480390e-01
+ -1.9369833e-01 3.3010575e-01
+ 2.3054934e-01 -3.6639404e-01
+ -3.1783106e-01 7.6021674e-02
+ -3.7923317e-01 4.4276069e-01
+ 5.7239164e-01 -5.4193935e-01
+ -2.2789189e-01 5.2795431e-02
+ -5.0327896e-01 0.0000000e+00
+ -2.0743778e-01 -3.6617248e-01
+ -3.2543847e-01 -1.0368304e-01
+ 5.4438282e-01 3.9557705e-01
+ -5.3966238e-01 1.8749005e-01
+ 1.9969446e-02 -3.7231073e-03
+ -2.6968245e-01 3.6470797e-01
+ 1.5572359e-01 4.3274816e-02
+ 1.6642223e-01 4.1230161e-03
+ -4.1109021e-02 2.9053651e-01
+ 3.8123994e-01 -5.1316546e-02
+ -2.0871338e-01 -3.2318002e-01
+ 3.7229772e-01 -2.7151345e-01
+ 2.6253907e-01 1.1768620e-01
+ 1.5729637e-01 -3.3093795e-02
+ -4.2836754e-01 8.1195635e-02
+ 1.3357679e-01 2.5044397e-01
+ 1.2352769e-01 2.9242938e-01
+ 1.9316889e-01 -3.3501492e-01
+ 1.2380095e-01 -2.4094826e-01
+ 1.8021701e-01 -1.2331446e-01
+ -3.5518384e-01 -2.6547013e-01
+ -4.0639237e-01 1.0703449e-01
+ -6.9082191e-02 2.7259751e-01
+ -1.6488603e-03 0.0000000e+00
+ 6.1735029e-02 3.4653391e-02
+ 2.1704776e-02 3.4263929e-01
+ 1.0872352e-01 2.2011490e-01
+ 1.3501932e-01 1.1506453e-01
+ 1.6093573e-01 4.8444486e-01
+ 8.1894565e-02 3.0759670e-02
+ 1.3709820e-01 -1.1116684e-01
+ 3.4835629e-01 -8.1579220e-02
+ 2.6171283e-01 4.0326011e-01
+ 2.4807415e-01 4.8061696e-01
+ -3.4332078e-01 -4.0254767e-01
+ 2.5970124e-01 -1.3969006e-01
+ 6.8328263e-02 -6.6354973e-02
+ -4.5122488e-01 -4.6269784e-02
+ -1.4368405e-01 3.6077946e-01
+ 1.9103192e-01 6.1505021e-02
+ 2.7002266e-01 1.4570465e-01
+ -8.3518301e-03 -2.3300428e-01
+ 1.0052413e-01 1.9614516e-01
+ -1.2789522e-01 1.9539943e-01
+ -1.2523768e-01 -3.1924289e-01
+ 9.0345359e-02 -8.8993370e-02
+ 1.3997243e-01 1.9943073e-01
+ 2.8620732e-01 -8.3369222e-02
+ 7.2487019e-02 0.0000000e+00
+ -1.4405643e-01 2.0703756e-01
+ 5.1364093e-01 2.1133332e-01
+ 2.1889104e-01 3.2346946e-01
+ 2.2977746e-01 6.7128928e-03
+ 2.4964012e-01 9.3280846e-03
+ 3.6795863e-01 3.7020049e-03
+ -2.1689136e-01 1.4973488e-01
+ 4.7090553e-02 4.1509150e-03
+ 6.8672872e-01 -5.1896898e-01
+ 2.0338075e-01 -9.9856633e-02
+ -4.2723522e-02 -2.3333240e-01
+ -1.7546362e-01 2.6162404e-01
+ -1.8006089e-01 2.6342849e-01
+ -4.5044011e-01 1.4835151e-01
+ 1.0052074e-01 1.6757867e-01
+ 2.9519066e-02 -2.9171202e-01
+ 3.5071014e-01 8.3758835e-02
+ 3.4762331e-02 -3.3651686e-01
+ -1.7556118e-01 -2.2172512e-01
+ -1.8053796e-01 -2.1776832e-02
+ -2.4177561e-01 -1.7544787e-01
+ -3.1968262e-01 8.7304604e-02
+ -1.9411642e-01 2.8576308e-01
+ 9.5900932e-02 -1.8985625e-01
+ -2.3725711e-01 -1.1186895e-01
+ 1.3322612e-01 0.0000000e+00
+ 3.2404337e-03 1.5033015e-01
+ -3.5860698e-02 2.6009940e-01
+ -3.4152502e-01 -9.1461626e-02
+ 4.3132465e-01 -4.5732916e-01
+ -2.9743337e-01 -1.7686635e-01
+ 2.1713395e-01 -2.3775328e-01
+ 3.5887511e-02 -1.0055884e-01
+ 7.8474825e-02 4.0269808e-02
+ 2.9607875e-01 -1.6104883e-02
+ -3.3844187e-01 -1.2800046e-01
+ 1.0462576e-01 -3.7881584e-02
+ -3.8464626e-01 5.2710939e-02
+ 3.9681901e-03 -3.4742573e-02
+ 1.7336716e-01 1.8234260e-01
+ 3.0204173e-01 -1.8760916e-01
+ 3.5852012e-02 -1.4469108e-01
+ 2.7139090e-01 -1.7898693e-01
+ -2.9251436e-01 1.1071905e-01
+ 5.1195180e-02 -8.5217152e-02
+ 1.4089344e-01 -2.7084618e-01
+ 2.0713480e-01 -4.3901445e-02
+ 2.3465076e-01 1.5973969e-01
+ -2.3148592e-02 -2.4246980e-01
+ 2.0416785e-01 3.3423231e-01
+ -8.5500225e-02 1.3130519e-02
+ 1.4879662e-02 4.6955367e-02
+ 8.1703180e-02 0.0000000e+00
+ -6.4695470e-02 -3.5887714e-02
+ 1.4924033e-02 -1.3935984e-02
+ -7.9831823e-02 -2.3701413e-01
+ -3.6159109e-02 1.8219454e-01
+ -3.3031056e-01 -3.3893427e-01
+ 1.1233064e-01 1.4498057e-01
+ 2.5168875e-01 1.5612661e-01
+ -1.8172447e-01 -1.9812508e-01
+ -6.9329238e-02 -2.5480816e-01
+ -3.0553243e-01 2.8667997e-02
+ -5.0571882e-02 2.1342699e-01
+ -2.4790052e-01 4.6800800e-02
+ 1.0336709e-01 1.0566485e-01
+ 3.4678273e-01 2.6730095e-01
+ 4.8209406e-02 -3.3262054e-02
+ 5.8154802e-02 7.1672295e-02
+ -1.0004239e-01 -1.1412671e-02
+ -4.6840109e-02 1.7811092e-01
+ 6.1699462e-03 5.4266122e-02
+ -2.8746028e-02 8.2760268e-02
+ -1.1131281e-01 1.5041323e-01
+ -1.3161903e-01 6.5361572e-02
+ 1.2030895e-01 2.6226935e-01
+ 8.6943264e-03 -5.0648282e-02
+ -2.7347079e-01 -1.3588011e-01
+ -1.5430779e-01 -5.3959255e-02
+ -1.8231483e-01 -2.4276404e-02
+ -2.2095634e-01 0.0000000e+00
+ 1.1932987e-01 -1.9712991e-01
+ -3.3935290e-01 -1.8665900e-01
+ -6.5635383e-02 -2.3025254e-01
+ -3.1952073e-02 2.5514473e-01
+ -2.7072772e-01 1.0677352e-01
+ -1.0426007e-01 2.9227563e-01
+ 5.8502652e-02 -1.2779416e-01
+ -6.9396199e-02 -3.9721662e-02
+ -2.0986135e-01 2.8705412e-01
+ -2.4028495e-01 1.8014485e-01
+ 8.0135754e-02 4.1571444e-02
+ -1.9611118e-02 2.3377602e-01
+ -4.8072465e-02 -1.3341153e-01
+ -1.8363542e-01 -2.9081433e-01
+ 2.7698353e-01 4.9423924e-02
+ -1.0489157e-01 -3.0025982e-01
+ -3.1192199e-01 1.1267751e-01
+ 1.2900241e-01 -7.5426452e-02
+ -1.4183103e-01 -5.4506433e-01
+ -3.4394343e-02 1.4693694e-01
+ -1.5323999e-01 -1.4958407e-01
+ -5.0405486e-02 -1.6528435e-01
+ -1.5575468e-01 -6.0322659e-02
+ 2.5008488e-01 -3.1390796e-01
+ -1.6653051e-01 3.9457516e-01
+ 2.7288167e-01 8.6527239e-02
+ 1.8162860e-01 -2.7302327e-02
+ 1.5406588e-01 1.5108586e-01
+ -1.1932299e-01 0.0000000e+00
+ -9.1078522e-02 3.0590469e-01
+ -5.3821973e-02 -5.5209791e-02
+ -8.6221952e-02 2.6175883e-01
+ -5.7405056e-01 5.4260089e-02
+ 1.8049595e-01 -1.3467397e-01
+ 2.6341004e-01 2.3471841e-01
+ 8.1165527e-02 7.1223772e-02
+ -4.4267166e-01 2.1657905e-01
+ -1.7441155e-03 -6.0891220e-02
+ 3.1109021e-01 1.2371573e-03
+ 1.4082972e-01 -1.5526721e-01
+ -1.7635212e-02 -1.0980670e-01
+ 3.2742204e-02 5.9624600e-02
+ -6.9596516e-02 -1.0677727e-01
+ 2.1443968e-01 1.5940691e-01
+ 3.4891352e-03 -3.3197113e-01
+ 5.1612088e-02 1.0652919e-01
+ -8.0288746e-02 -6.3806019e-02
+ 1.4044482e-01 -1.5413007e-01
+ -1.8363936e-01 1.3081697e-01
+ 2.2520087e-01 1.1620847e-01
+ 2.6476532e-01 -2.8001993e-02
+ 4.3219082e-02 -5.0147129e-02
+ 7.0727572e-03 -5.0630769e-02
+ -1.3097939e-01 -1.9766550e-01
+ 1.7784148e-01 -1.5408588e-01
+ 1.7482924e-01 1.5994048e-02
+ 2.1557115e-01 -1.2973634e-01
+ -2.8950633e-01 1.1785615e-01
+ 1.4145228e-01 0.0000000e+00
+ -6.1719402e-03 -2.4578660e-02
+ -2.2524324e-01 -2.6530928e-02
+ -8.8599099e-02 3.1056733e-01
+ -3.7662674e-02 -8.4962799e-02
+ 5.9524435e-02 1.2828934e-01
+ -3.5985441e-03 1.8913083e-02
+ -5.4099685e-02 -8.6208276e-02
+ 1.0525864e-01 7.3945161e-02
+ 1.8137023e-01 2.6474035e-01
+ 1.1053889e-01 -1.1206217e-01
+ 2.1592062e-01 -2.3487375e-01
+ 3.8659356e-01 -2.4273926e-01
+ -3.1187313e-01 -6.7359088e-02
+ 1.0932508e-01 1.8364109e-01
+ -1.4820920e-02 2.4038158e-02
+ -2.3889150e-01 9.3674429e-02
+ 1.5567563e-01 1.3289161e-01
+ -2.5238407e-01 -1.8016021e-01
+ 3.0822838e-01 -6.2169727e-02
+ -1.1807753e-01 3.0161424e-01
+ 2.6366804e-01 1.2535790e-01
+ -1.1854330e-01 -1.9689787e-01
+ -1.2192421e-01 2.4814268e-01
+ -6.1896829e-02 -5.8285696e-02
+ -6.9572481e-02 3.5543989e-01
+ 3.0650678e-02 2.8377261e-01
+ 1.7476068e-01 -2.8382568e-01
+ -1.3169714e-01 -1.8298210e-01
+ -8.9047092e-02 -4.2433656e-02
+ 5.8530542e-02 1.9159661e-01
+ 1.5283284e-01 0.0000000e+00
+ -1.6100197e-01 4.1776383e-01
+ 1.4157277e-01 1.1569515e-01
+ 1.4912327e-01 1.7054802e-01
+ 2.5754950e-01 -1.1703479e-01
+ 2.5959565e-01 -1.2337952e-01
+ -1.9241755e-03 -6.6799069e-03
+ 3.4074349e-02 3.0342366e-02
+ -7.0836315e-02 4.8380438e-02
+ -2.0471655e-02 -1.3224747e-01
+ 9.8796186e-02 -1.1061616e-01
+ -4.4909054e-02 -4.9073517e-01
+ 1.2965930e-02 9.0343083e-02
+ -2.2060312e-01 -1.0551363e-01
+ -1.8213150e-01 7.3123342e-02
+ -9.4007717e-02 5.3443604e-02
+ -1.5554344e-01 1.3891478e-01
+ 6.6568288e-02 -1.5496309e-01
+ -1.1959706e-02 -5.4871489e-02
+ -4.8384732e-02 -1.0783429e-01
+ -5.5372950e-02 1.2962901e-01
+ 2.2509866e-01 -1.6801958e-01
+ -2.2362351e-01 -2.5387974e-01
+ -1.8231359e-01 -9.7324607e-02
+ -6.5384398e-02 -8.4011488e-02
+ 3.7616019e-01 3.9852540e-02
+ -2.9051514e-01 4.0287411e-02
+ 2.5715213e-02 -2.4951215e-01
+ 2.3090840e-01 6.4346974e-02
+ 4.5750285e-02 5.2584868e-02
+ -1.8661596e-02 -1.8393582e-01
+ 1.9954284e-01 3.5717752e-02
Added: mc/1D/hc/trunk/hc.c
===================================================================
--- mc/1D/hc/trunk/hc.c (rev 0)
+++ mc/1D/hc/trunk/hc.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -0,0 +1,293 @@
+#include "hc.h"
+/*
+
+
+implementation of Hager & O'Connell (1981) method of solving mantle
+circulation given internal density anomalies, only radially varying
+viscosity, and either free-slip or plate velocity boundary condition
+at surface
+
+based on Hager & O'Connell (1981), Hager & Clayton (1989), and
+Steinberger (2000)
+
+the original code is due to Brad Hager, Rick O'Connell, and was
+largely modified by Bernhard Steinberger
+
+this version by Thorsten Becker (twb at usc.edu)
+
+
+>>> SOME COMMENTS FROM THE ORIGINAL CODE <<<
+
+C * It uses the following Numerical Recipes (Press et al.) routines:
+C four1, realft, gauleg, rk4, rkdumb, ludcmp, lubksb;
+C !!!!!!!!!!!!!!!!!!! rkqc, odeint !!!!!!!!!! take out !!!!!!
+C and the following routines by R.J. O'Connell:
+C shc2d, shd2c, ab2cs, cs2ab, plmbar, plvel2sh, pltgrid, pltvel,
+C vshd2c, plmbar1.
+C Further subroutines are: kinsub, evalpa,
+C torsol (all based on "kinflow" by Hager & O'Connell),
+C densub and evppot (based on "denflow" by Hager & O'Connell),
+C sumsub (based on "sumkd" by Hager & O'Connell, but
+C substantially speeded up),
+C convert, derivs and shc2dd (based on R.J. O'Connell's shc2d).
+C
+C bugs found:
+C * The combination of (1) high viscosity lithosphere
+C (2) compressible flow (3) kinematic (plate driven) flow
+C doesn't work properly. The problem presumably only occurs
+C at degree 1 (I didn't make this sure) but this is sufficient
+C to screw up everything. It will usually work to reduce the
+C contrast between lithospheric and asthenospheric viscosity.
+C Then make sure that (1) two viscosity structures give similar
+C results for incompressible and (2) incompressible and compressible
+C reduced viscosity look similar (e.g. anomalous mass flux vs. depth)
+
+<<< END OLD COMMENTS
+
+*/
+
+int main(int argc, char **argv)
+{
+ struct hcs *model; /* main structure, make sure to initialize with
+ zeroes */
+ struct sh_lms *sol_spectral=NULL, *geoid = NULL; /* solution expansions */
+ struct sh_lms *pvel=NULL; /* local plate velocity expansion */
+ int nsol,lmax,i;
+ FILE *out;
+ struct hc_parameters p[1]; /* parameters */
+ char filename[HC_CHAR_LENGTH],file_prefix[10];
+ HC_PREC *sol_spatial = NULL; /* spatial solution,
+ e.g. velocities */
+ HC_PREC corr[2]; /* correlations */
+ static hc_boolean geoid_binary = FALSE; /* type of geoid output */
+ static HC_CPREC unitya[1] = {1.0};
+ /*
+
+
+ (1)
+
+ initialize the model structure, this is needed to initialize some
+ of the default values before callign the parameter handling
+ routine this call also involves initializing the hc parameter
+ structure
+
+ */
+ hc_struc_init(&model);
+ /*
+
+ (2)
+ init parameters to default values
+
+ */
+ hc_init_parameters(p);
+ /*
+ handle command line arguments
+ */
+ hc_handle_command_line(argc,argv,1,p);
+ /*
+
+ begin main program part
+
+ */
+#ifdef __TIMESTAMP__
+ if(p->verbose)
+ fprintf(stderr,"%s: starting version compiled on %s\n",argv[0],__TIMESTAMP__);
+#else
+ if(p->verbose)
+ fprintf(stderr,"%s: starting main program\n",argv[0]);
+#endif
+ /*
+
+ (3)
+
+ initialize all variables
+
+ - choose the internal spherical harmonics convention
+ - assign constants
+ - assign phase boundaries, if any
+ - read in viscosity structure
+ - assign density anomalies
+ - read in plate velocities
+
+ */
+ hc_init_main(model,SH_RICK,p);
+ nsol = (model->nradp2) * 3; /*
+ number of solutions (r,pol,tor) * (nlayer+2)
+
+ total number of layers is nlayer +2,
+
+ because CMB and surface are added
+ to intermediate layers which are
+ determined by the spacing of the
+ density model
+
+ */
+ if(p->free_slip) /* maximum degree is determined by the
+ density expansion */
+ lmax = model->dens_anom[0].lmax;
+ else /* max degree is determined by the
+ plate velocities */
+ lmax = model->pvel.p[0].lmax; /* shouldn't be larger than that*/
+ /*
+ make sure we have room for the plate velocities
+ */
+ sh_allocate_and_init(&pvel,2,lmax,model->sh_type,1,p->verbose,FALSE);
+
+ /* init done */
+ /*
+
+
+
+ SOLUTION PART
+
+
+ */
+ /*
+ make room for the spectral solution on irregular grid
+ */
+ sh_allocate_and_init(&sol_spectral,nsol,lmax,model->sh_type,HC_VECTOR,
+ p->verbose,FALSE);
+ if(p->compute_geoid == 1)
+ /* make room for geoid solution at surface */
+ sh_allocate_and_init(&geoid,1,model->dens_anom[0].lmax,
+ model->sh_type,HC_SCALAR,p->verbose,FALSE);
+ else if(p->compute_geoid == 2) /* all layers */
+ sh_allocate_and_init(&geoid,model->nradp2,model->dens_anom[0].lmax,
+ model->sh_type,HC_SCALAR,p->verbose,FALSE);
+
+ /*
+ solve poloidal and toroidal part and sum
+ */
+ if(!p->free_slip)
+ hc_select_pvel(p->pvel_time,&model->pvel,pvel,p->verbose);
+ hc_solve(model,p->free_slip,p->solution_mode,sol_spectral,
+ TRUE,TRUE,TRUE,p->print_pt_sol,p->compute_geoid,
+ pvel,model->dens_anom,geoid,
+ p->verbose);
+ /*
+
+ OUTPUT PART
+
+ */
+ /*
+
+ output of spherical harmonics solution
+
+ */
+ switch(p->solution_mode){
+ case HC_VEL:
+ sprintf(file_prefix,"vel");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;
+ }
+ if(p->sol_binary_out)
+ sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_BINARY);
+ else
+ sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_ASCII);
+ if(p->verbose)
+ fprintf(stderr,"%s: writing spherical harmonics solution to %s\n",
+ argv[0],filename);
+ out = ggrd_open(filename,"w","main");
+ hc_print_spectral_solution(model,sol_spectral,out,
+ p->solution_mode,
+ p->sol_binary_out,p->verbose);
+ fclose(out);
+ /* */
+ if(p->print_density_field){
+ /*
+ print the density field
+ */
+ sprintf(file_prefix,"dscaled");
+ if(p->sol_binary_out)
+ sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_BINARY);
+ else
+ sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_ASCII);
+ if(p->verbose)
+ fprintf(stderr,"%s: writing scaled density anomaly field to %s\n",
+ argv[0],filename);
+
+ out = ggrd_open(filename,"w","main");
+ hc_print_dens_anom(model,out,p->sol_binary_out,p->verbose);
+ fclose(out);
+ }
+
+
+ /* compute the geoid? */
+ if(p->compute_geoid){
+ if(p->compute_geoid_correlations){
+ if(p->compute_geoid == 2){ /* check if all geoids were computed */
+ fprintf(stderr,"%s: ERROR: can only compute correlation for surface geoid, geoid = %i\n",
+ argv[0],p->compute_geoid);
+ exit(-1);
+ }
+ if(p->verbose)
+ fprintf(stderr,"%s: correlation for geoid with %s\n",argv[0],p->ref_geoid_file);
+ hc_compute_correlation(geoid,p->ref_geoid,corr,1,p->verbose);
+ fprintf(stdout,"%10.7f %10.7f\n",(double)corr[0],(double)corr[1]);
+ }else{
+ /*
+ print geoid solution
+ */
+ sprintf(filename,"%s",HC_GEOID_FILE);
+ if(p->verbose)
+ fprintf(stderr,"%s: writing geoid to %s, %s\n",
+ argv[0],filename,(p->compute_geoid == 1)?("at surface"):("all layers"));
+ out = ggrd_open(filename,"w","main");
+ if(p->compute_geoid == 1) /* surface layer */
+ hc_print_sh_scalar_field(geoid,out,FALSE,geoid_binary,p->verbose);
+ else{ /* all layers */
+ for(i=0;i < model->nradp2;i++){
+ sh_print_parameters_to_stream((geoid+i),1,i,model->nradp2,
+ HC_Z_DEPTH(model->r[i]),out,FALSE,geoid_binary,p->verbose);
+ sh_print_coefficients_to_stream((geoid+i),1,out,unitya,geoid_binary,p->verbose);
+ }
+ }
+ fclose(out);
+ }
+ }
+ if(p->print_spatial){
+ /*
+ we wish to use the spatial solution
+
+ expand velocities to spatial base, compute spatial
+ representation
+
+ */
+ hc_compute_sol_spatial(model,sol_spectral,&sol_spatial,
+ p->verbose);
+ /*
+
+ output of spatial solution
+
+ */
+ sprintf(filename,"%s.%s",file_prefix,HC_SPATIAL_SOLOUT_FILE);
+ /* print lon lat z v_r v_theta v_phi */
+ hc_print_spatial_solution(model,sol_spectral,sol_spatial,
+ filename,HC_LAYER_OUT_FILE,
+ p->solution_mode,p->sol_binary_out,
+ p->verbose);
+ }
+ /*
+
+ free memory
+
+ */
+ sh_free_expansion(sol_spectral,nsol);
+ /* local copies of plate velocities */
+ sh_free_expansion(pvel,2);
+ /* */
+ if(p->compute_geoid == 1) /* only one layer */
+ sh_free_expansion(geoid,1);
+ else if(p->compute_geoid == 2) /* all layers */
+ sh_free_expansion(geoid,model->nradp2);
+ free(sol_spatial);
+ if(p->verbose)
+ fprintf(stderr,"%s: done\n",argv[0]);
+ hc_struc_free(&model);
+ return 0;
+}
+
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc.h 2013-06-21 23:55:36 UTC (rev 22400)
@@ -226,7 +226,7 @@
*/
HC_PREC *rdf,*sdf;
int ndf;
- struct sh_lms *ref_geoid;
+ struct sh_lms *ref_geoid,*ref_dtopo;
HC_PREC rlayer[3]; /* for four layer approaches (first
layer radius is CMB radius) */
@@ -237,6 +237,7 @@
char prem_model_filename[HC_CHAR_LENGTH]; /* PREM model filename */
char dens_scaling_filename[HC_CHAR_LENGTH]; /* */
char ref_geoid_file[HC_CHAR_LENGTH]; /* reference geoid */
+ char ref_dtopo_file[HC_CHAR_LENGTH]; /* reference dynamic topography */
};
/* plate velocity structure */
@@ -281,7 +282,7 @@
int inho; /* number of density layers */
HC_PREC *dfact; /* density factors, derived from layer thickness */
HC_PREC *rden; /* radii of density layers [normalized] */
-
+ HC_PREC rho_top_kg; /* top density in kg/m^3 */
struct sh_lms *dens_anom; /*
expansions of density
anomalies has to be [inho] (if
@@ -421,6 +422,7 @@
#define HC_SOLVER_MODE_DEFAULT 0
#define HC_SOLVER_MODE_VISC_SCAN 1
+#define HC_SOLVER_MODE_DYNTOPO_INVERT 2
/*
Deleted: mc/1D/hc/trunk/hc.tar
===================================================================
(Binary files differ)
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2013-06-21 23:55:36 UTC (rev 22400)
@@ -34,6 +34,7 @@
void ggrd_get_velocities(double *, double *, double *, int, struct ggrd_master *, double, double);
void ggrd_weights(double, double *, int, int, double [(5 +1)][(1 +1)]);
/* grdinttester.c */
+/* hc.c */
/* hc_extract_sh_layer.c */
/* hc_extract_spatial.c */
/* hc_init.c */
@@ -42,7 +43,7 @@
void hc_init_polsol_struct(struct hc_ps *);
void hc_init_main(struct hcs *, int, struct hc_parameters *);
void hc_init_constants(struct hcs *, long double, char *, unsigned short);
-void hc_handle_command_line(int, char **, struct hc_parameters *);
+void hc_handle_command_line(int, char **, int, struct hc_parameters *);
void hc_assign_viscosity(struct hcs *, int, long double [4], struct hc_parameters *);
void hc_assign_density(struct hcs *, unsigned short, int, char *, int, unsigned short, unsigned short, unsigned short, unsigned short, unsigned short, unsigned short, int, long double *, long double *, unsigned short);
long double hc_find_dens_scale(long double, long double, unsigned short, long double *, long double *, int);
@@ -53,10 +54,11 @@
void hc_get_blank_expansions(struct sh_lms **, int, int, char *);
void hc_struc_free(struct hcs **);
void hc_assign_dd_scaling(int, long double [4], struct hc_parameters *, long double);
-void hc_read_geoid(struct hc_parameters *);
+void hc_read_scalar_shexp(char *, struct sh_lms **, char *, struct hc_parameters *);
void hc_select_pvel(long double, struct pvels *, struct sh_lms *, unsigned short);
/* hc_input.c */
int hc_read_sh_solution(struct hcs *, struct sh_lms **, FILE *, unsigned short, unsigned short);
+/* hc_invert_dtopo.c */
/* hc_matrix.c */
void hc_ludcmp_3x3(long double [3][3], int, int *);
void hc_lubksb_3x3(long double [3][3], int, int *, long double *);
@@ -111,7 +113,7 @@
void hc_flipit(void *, void *, size_t);
void hc_print_dens_anom(struct hcs *, FILE *, unsigned short, unsigned short);
/* hc_polsol.c */
-void hc_polsol(struct hcs *, int, long double *, int, long double *, unsigned short, struct sh_lms *, unsigned short, int, long double *, long double *, unsigned short, struct sh_lms *, struct sh_lms *, unsigned short, struct sh_lms *, unsigned short, unsigned short);
+void hc_polsol(struct hcs *, int, long double *, int, long double *, unsigned short, struct sh_lms *, unsigned short, int, long double *, long double *, unsigned short, struct sh_lms *, struct sh_lms *, unsigned short, struct sh_lms *, unsigned short, unsigned short, unsigned short);
/* hc_propagator.c */
void hc_evalpa(int, long double, long double, long double, long double *);
void hc_evppot(int, long double, long double *);
@@ -119,9 +121,10 @@
void hc_solve(struct hcs *, unsigned short, int, struct sh_lms *, unsigned short, unsigned short, unsigned short, unsigned short, unsigned short, struct sh_lms *, struct sh_lms *, struct sh_lms *, unsigned short);
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 *, long double **, unsigned short);
+void hc_compute_dynamic_topography(struct hcs *, struct sh_lms *, struct sh_lms **, unsigned short, unsigned short);
/* hc_torsol.c */
void hc_torsol(struct hcs *, int, int, int, long double *, long double **, long double **, struct sh_lms *, struct sh_lms *, long double *, unsigned short);
-/* main.c */
+/* hc_visc_scan.c */
/* prem2dsm.c */
/* prem_util.c */
int prem_find_layer_x(double, double, double *, int, int, double *);
@@ -178,12 +181,14 @@
void sh_compute_power_per_degree(struct sh_lms *, long double *);
long double sh_correlation(struct sh_lms *, struct sh_lms *, int);
long double sh_correlation_per_degree(struct sh_lms *, struct sh_lms *, int, int);
-void sh_print_parameters_to_file(struct sh_lms *, int, int, int, long double, FILE *, unsigned short, unsigned short, unsigned short);
-unsigned short sh_read_parameters_from_file(int *, int *, int *, int *, int *, long double *, int *, FILE *, unsigned short, unsigned short, unsigned short);
-void sh_print_coefficients_to_file(struct sh_lms *, int, FILE *, long double *, unsigned short, unsigned short);
-void sh_read_coefficients_from_file(struct sh_lms *, int, int, FILE *, unsigned short, long double *, unsigned short);
+void sh_single_par_and_exp_to_file(struct sh_lms *, char *, unsigned short, unsigned short);
+void sh_single_par_and_exp_to_stream(struct sh_lms *, FILE *, unsigned short, unsigned short);
+void sh_print_parameters_to_stream(struct sh_lms *, int, int, int, long double, FILE *, unsigned short, unsigned short, unsigned short);
+unsigned short sh_read_parameters_from_stream(int *, int *, int *, int *, int *, long double *, int *, FILE *, unsigned short, unsigned short, unsigned short);
+void sh_print_coefficients_to_stream(struct sh_lms *, int, FILE *, long double *, unsigned short, unsigned short);
+void sh_read_coefficients_from_stream(struct sh_lms *, int, int, FILE *, unsigned short, long double *, unsigned short);
void sh_print_nonzero_coeff(struct sh_lms *, FILE *);
-void sh_read_spatial_data_from_file(struct sh_lms *, FILE *, unsigned short, int, long double *, long double *);
+void sh_read_spatial_data_from_stream(struct sh_lms *, FILE *, unsigned short, int, long double *, long double *);
void sh_read_spatial_data_from_grd(struct sh_lms *, struct ggrd_gt *, unsigned short, int, long double *, long double *);
void sh_read_spatial_data(struct sh_lms *, FILE *, struct ggrd_gt *, unsigned short, unsigned short, int, long double *, long double *);
void sh_compute_spatial_basis(struct sh_lms *, FILE *, unsigned short, long double, long double **, int, unsigned short);
@@ -193,10 +198,10 @@
void sh_compute_spatial_irreg(struct sh_lms *, int, long double *, long double *, int, long double *, unsigned short);
void sh_exp_type_error(char *, struct sh_lms *);
void sh_print_plm(long double *, int, int, int, FILE *);
-void sh_print_spatial_data_to_file(struct sh_lms *, int, long double *, unsigned short, long double, FILE *);
+void sh_print_spatial_data_to_stream(struct sh_lms *, int, long double *, unsigned short, long double, FILE *);
void sh_get_coordinates(struct sh_lms *, int, long double *, long double *);
-void sh_print_reg_spatial_data_to_file(struct sh_lms *, int, long double *, unsigned short, long double, long double *, int, long double *, int, FILE *);
-void sh_print_irreg_spatial_data_to_file(struct sh_lms *, int, long double *, unsigned short, long double, long double *, long double *, int, FILE *);
+void sh_print_reg_spatial_data_to_stream(struct sh_lms *, int, long double *, unsigned short, long double, long double *, int, long double *, int, FILE *);
+void sh_print_irreg_spatial_data_to_stream(struct sh_lms *, int, long double *, unsigned short, long double, long double *, long double *, int, FILE *);
void sh_compute_plm(struct sh_lms *, int, long double **, unsigned short);
void sh_compute_plm_reg(struct sh_lms *, int, long double **, unsigned short, long double *, int);
void sh_get_coeff(struct sh_lms *, int, int, int, unsigned short, long double *);
Modified: mc/1D/hc/trunk/hc_constants.h
===================================================================
--- mc/1D/hc/trunk/hc_constants.h 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc_constants.h 2013-06-21 23:55:36 UTC (rev 22400)
@@ -51,3 +51,6 @@
use 0.01, if input/output is to be in cm/yr
*/
#define HC_VEL_IO_SCALE 0.01
+
+
+
Modified: mc/1D/hc/trunk/hc_extract_sh_layer.c
===================================================================
--- mc/1D/hc/trunk/hc_extract_sh_layer.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc_extract_sh_layer.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -114,9 +114,9 @@
/* SH header */
if(short_format && loop)
fprintf(stdout,"%g\n",(double)HC_Z_DEPTH(model->r[ilayer]));
- sh_print_parameters_to_file((sol+ilayer*shps_read),shps,
- ilayer,nset,(HC_PREC)(HC_Z_DEPTH(model->r[ilayer])),
- stdout,short_format,FALSE,verbose);
+ sh_print_parameters_to_stream((sol+ilayer*shps_read),shps,
+ ilayer,nset,(HC_PREC)(HC_Z_DEPTH(model->r[ilayer])),
+ stdout,short_format,FALSE,verbose);
}
switch(mode){
case 1:
@@ -124,21 +124,21 @@
if(verbose)
fprintf(stderr,"%s: printing x_r SHE at layer %i (depth: %g)\n",
argv[0],ilayer,(double)HC_Z_DEPTH(model->r[ilayer]));
- sh_print_coefficients_to_file((sol+ilayer*shps_read),shps,stdout,fac,FALSE,verbose);
+ sh_print_coefficients_to_stream((sol+ilayer*shps_read),shps,stdout,fac,FALSE,verbose);
break;
case 2:
/* */
if(verbose)
fprintf(stderr,"%s: printing x_pol x_tor SHE at layer %i (depth: %g)\n",
argv[0],ilayer,(double)HC_Z_DEPTH(model->r[ilayer]));
- sh_print_coefficients_to_file((sol+ilayer*shps_read+1),shps,stdout,fac,FALSE,verbose);
+ sh_print_coefficients_to_stream((sol+ilayer*shps_read+1),shps,stdout,fac,FALSE,verbose);
break;
case 3:
/* mode == 3 */
if(verbose)
fprintf(stderr,"%s: printing x_r x_pol x_tor SHE at layer %i (depth: %g)\n",
argv[0],ilayer,(double)HC_Z_DEPTH(model->r[ilayer]));
- sh_print_coefficients_to_file((sol+ilayer*shps_read),shps,stdout,fac,FALSE,verbose);
+ sh_print_coefficients_to_stream((sol+ilayer*shps_read),shps,stdout,fac,FALSE,verbose);
break;
case 4:
fprintf(stdout,"%5i %11g\n",ilayer,(double)HC_Z_DEPTH(model->r[ilayer]));
@@ -148,14 +148,14 @@
if(verbose)
fprintf(stderr,"%s: printing x_pol SHE at layer %i (depth: %g)\n",
argv[0],ilayer,(double)HC_Z_DEPTH(model->r[ilayer]));
- sh_print_coefficients_to_file((sol+ilayer*shps_read+1),shps,stdout,fac,FALSE,verbose);
+ sh_print_coefficients_to_stream((sol+ilayer*shps_read+1),shps,stdout,fac,FALSE,verbose);
break;
case 6:
/* */
if(verbose)
fprintf(stderr,"%s: printing x_tor SHE at layer %i (depth: %g)\n",
argv[0],ilayer,(double)HC_Z_DEPTH(model->r[ilayer]));
- sh_print_coefficients_to_file((sol+ilayer*shps_read+2),shps,stdout,fac,FALSE,verbose);
+ sh_print_coefficients_to_stream((sol+ilayer*shps_read+2),shps,stdout,fac,FALSE,verbose);
break;
default:
Modified: mc/1D/hc/trunk/hc_extract_spatial.c
===================================================================
--- mc/1D/hc/trunk/hc_extract_spatial.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc_extract_spatial.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -149,21 +149,21 @@
(double)zlabel);
ivec=FALSE;sh_compute_spatial((vsol+ilayer*shps_read),ivec,TRUE,&plm,data,verbose);
- sh_print_spatial_data_to_file((vsol+ilayer*shps_read),shps,data,TRUE,zlabel,stdout);
+ sh_print_spatial_data_to_stream((vsol+ilayer*shps_read),shps,data,TRUE,zlabel,stdout);
break;
case 2:
/* */
if(verbose)
fprintf(stderr,"%s: printing v_theta v_phi SHE at layer %i (depth: %g)\n",argv[0],ilayer,(double)zlabel);
ivec=TRUE;sh_compute_spatial((vsol+ilayer*shps_read+1),ivec,TRUE,&plm,data,verbose);
- sh_print_spatial_data_to_file((vsol+ilayer*shps_read+1),shps,data,TRUE,zlabel,stdout);
+ sh_print_spatial_data_to_stream((vsol+ilayer*shps_read+1),shps,data,TRUE,zlabel,stdout);
break;
case 3:
if(verbose)
fprintf(stderr,"%s: printing v_r v_theta v_phi SHE at layer %i (depth: %g)\n",argv[0],ilayer,(double)zlabel);
ivec=FALSE;sh_compute_spatial((vsol+ilayer*shps_read), ivec,TRUE,&plm,data,verbose); /* radial */
ivec=TRUE; sh_compute_spatial((vsol+ilayer*shps_read+1),ivec,TRUE,&plm,(data+npoints),verbose); /* theta,phi */
- sh_print_spatial_data_to_file((vsol+ilayer*shps_read),shps,data,TRUE,zlabel,stdout);
+ sh_print_spatial_data_to_stream((vsol+ilayer*shps_read),shps,data,TRUE,zlabel,stdout);
break;
case 4:
fprintf(stdout,"%5i %11g\n",ilayer,(double)HC_Z_DEPTH(model->r[ilayer]));
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc_init.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -356,12 +356,29 @@
visc_filename[] needs to be [HC_CHAR_LENGTH]
*/
-void hc_handle_command_line(int argc, char **argv,
+void hc_handle_command_line(int argc, char **argv,int start_from_i,
struct hc_parameters *p)
{
int i;
- for(i=1;i < argc;i++){
- if(strcmp(argv[i],"-help")==0 || strcmp(argv[i],"-h")==0 || strcmp(argv[i],"-?")==0){// help
+ hc_boolean used_parameter;
+ for(i=start_from_i;i < argc;i++){
+ used_parameter = FALSE; /* */
+ if((p->solver_mode == HC_SOLVER_MODE_VISC_SCAN && argc < 2) || /* need
+ one
+ or
+ two
+ arguments
+ for
+ some
+ of
+ the
+ other
+ modes */
+ (p->solver_mode == HC_SOLVER_MODE_DYNTOPO_INVERT && argc < 3) ||
+ (strcmp(argv[i],"-help")==0) ||
+ (strcmp(argv[i],"--help")==0) ||
+ (strcmp(argv[i],"-h")==0) ||
+ (strcmp(argv[i],"-?")==0)){// help
/*
help page
*/
@@ -372,16 +389,30 @@
fprintf(stderr,"This particular implementation illustrates one possible way to combine the HC solver routines.\n");
fprintf(stderr,"Based on code by Brad Hager, Richard O'Connell, and Bernhard Steinberger.\n");
fprintf(stderr,"This version by Thorsten Becker and Craig O'Neill\n\n");
-
- fprintf(stderr,"usage example:\n\n");
- fprintf(stderr,"bin/hc -vvv\n\n");
- fprintf(stderr,"Compute mantle flow solution using the default input files:\n");
- fprintf(stderr," viscosity profile visc.dat\n");
- fprintf(stderr," density profile dens.sh.dat\n");
- fprintf(stderr," earth model prem/prem.dat\n");
- fprintf(stderr,"and provide lots of output. Default setting is quiet operation.\n\n");
+ switch(p->solver_mode){
+ case HC_SOLVER_MODE_VISC_SCAN:
+ fprintf(stderr,"usage example:\n\n");
+ fprintf(stderr,"bin/hc_visc_scan geoid.ab\n\n");
+ fprintf(stderr,"Scan through viscosity values and compute correlation with the geoid in geoid.ab\n\n");
+ break;
+ case HC_SOLVER_MODE_DYNTOPO_INVERT:
+ fprintf(stderr,"usage example:\n\n");
+ fprintf(stderr,"bin/hc_invert_dtopo geoid.ab dtopo.ab\n\n");
+ fprintf(stderr,"Invert for density anomalies based on geoid in geoid.ab and res topo wrt to air in dtopo.ab\n\n");
+ break;
+ default:
+ fprintf(stderr,"usage example:\n\n");
+ fprintf(stderr,"bin/hc -vvv\n\n");
+ fprintf(stderr,"Compute mantle flow solution using the default input files:\n");
+ fprintf(stderr," viscosity profile visc.dat\n");
+ fprintf(stderr," density profile dens.sh.dat\n");
+ fprintf(stderr," earth model prem/prem.dat\n");
+ fprintf(stderr,"and provide lots of output. Default setting is quiet operation.\n\n");
+ break;
+ }
fprintf(stderr,"See README.TXT in the installation directory for example for how to plot output, and\n");
- fprintf(stderr,"http://geosys.usc.edu/projects/seatree/ for a graphical user interface.\n\n");
+ fprintf(stderr,"http://geosys.usc.edu/projects/seatree/ for a graphical user interface.\n");
+ fprintf(stderr,"http://geodynamics.usc.edu/~becker/ugesce.html for a VirtualBox install.\n\n");
fprintf(stderr,"density anomaly options:\n");
@@ -409,7 +440,9 @@
fprintf(stderr,"-vf\tname\tviscosity structure filename (%s), use pvisc.py to edit\n",
p->visc_filename);
fprintf(stderr,"\t\tThis file is in non_dim_radius viscosity[Pas] format\n");
- fprintf(stderr,"\t\tif name is \"visc_scan\", will loop through a four layer viscosity scan\n\n");
+ if(p->solver_mode == HC_SOLVER_MODE_VISC_SCAN)
+ fprintf(stderr,"\t\tWARNING: Will loop through a four layer viscosity scan\n\n");
+
fprintf(stderr,"boundary condition options:\n");
fprintf(stderr,"-fs\t\tperform free slip computation (%s)\n",hc_name_boolean(p->free_slip));
fprintf(stderr,"-ns\t\tperform no slip computation (%s)\n",hc_name_boolean(p->no_slip));
@@ -426,18 +459,21 @@
fprintf(stderr,"solution procedure and I/O options:\n");
fprintf(stderr,"-cbckl\tval\twill modify CMB boundary condition for all l > val with solver kludge (%i)\n",
p->solver_kludge_l);
- fprintf(stderr,"-ng\t\tdo not compute and print the geoid (%i)\n",
- p->compute_geoid);
- fprintf(stderr,"-ag\t\tcompute geoid at all layer depths, as opposed to the surface only\n");
- fprintf(stderr,"-rg\tname\tcompute correlation of surface geoid with that in file \"name\",\n");
- fprintf(stderr,"\t\tthis will not print out the geoid file, but only correlations (%s)\n",
- hc_name_boolean(p->compute_geoid_correlations));
- fprintf(stderr,"-pptsol\t\tprint pol[6] and tor[2] solution vectors (%s)\n",
- 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,"-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");
+ if(p->solver_mode == HC_SOLVER_MODE_DEFAULT){
+ /* these only apply to the regular mode */
+ fprintf(stderr,"-ng\t\tdo not compute and print the geoid (%i)\n",
+ p->compute_geoid);
+ fprintf(stderr,"-ag\t\tcompute geoid at all layer depths, as opposed to the surface only\n");
+ fprintf(stderr,"-rg\tname\tcompute correlation of surface geoid with that in file \"name\",\n");
+ fprintf(stderr,"\t\tthis will not print out the geoid file, but only correlations (%s)\n",
+ hc_name_boolean(p->compute_geoid_correlations));
+ fprintf(stderr,"-pptsol\t\tprint pol[6] and tor[2] solution vectors (%s)\n",
+ 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,"-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");
@@ -445,82 +481,108 @@
}else if(strcmp(argv[i],"-ds")==0){ /* density anomaly scaling factor */
hc_advance_argument(&i,argc,argv);
sscanf(argv[i],HC_FLT_FORMAT,&p->dens_anom_scale);
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-cbckl")==0){ /* solver kludge */
hc_advance_argument(&i,argc,argv);
sscanf(argv[i],"%i",&p->solver_kludge_l);
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-vtime")==0){ /* */
hc_advance_argument(&i,argc,argv);
sscanf(argv[i],HC_FLT_FORMAT,&p->pvel_time);
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-fs")==0){ /* free slip flag */
p->free_slip = TRUE;p->no_slip = FALSE;
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-ns")==0){ /* no slip flag */
p->no_slip = TRUE;p->free_slip = FALSE;
- }else if(strcmp(argv[i],"-px")==0){ /* print spatial solution? */
- hc_toggle_boolean(&p->print_spatial);
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-dshs")==0){ /* use short format for densities ? */
hc_toggle_boolean(&p->read_short_dens_sh);
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-vshs")==0){ /* use short format for velocities? */
hc_toggle_boolean(&p->read_short_pvel_sh);
- }else if(strcmp(argv[i],"-pptsol")==0){ /* print
- poloidal/toroidal
- solution
- parameters */
- hc_toggle_boolean(&p->print_pt_sol);
- }else if(strcmp(argv[i],"-ng")==0){ /* do not compute geoid */
- p->compute_geoid = 0;
- }else if(strcmp(argv[i],"-ag")==0){ /* compute geoid at all layers */
- p->compute_geoid = 2;
- }else if(strcmp(argv[i],"-rg")==0){ /* compute geoid correlations */
- hc_toggle_boolean(&p->compute_geoid_correlations);
- p->compute_geoid = TRUE;
- hc_advance_argument(&i,argc,argv); /* filename */
- strncpy(p->ref_geoid_file,argv[i],HC_CHAR_LENGTH);
- hc_read_geoid(p);
- }else if(strcmp(argv[i],"-vf")==0){ /* viscosity filename */
- hc_advance_argument(&i,argc,argv);
- strncpy(p->visc_filename,argv[i],HC_CHAR_LENGTH);
- if(strcmp(p->visc_filename,"visc_scan")==0)
- { /* inversion mode */
- p->solver_mode = HC_SOLVER_MODE_VISC_SCAN;
- p->visc_init_mode = HC_INIT_E_FOUR_LAYERS;
- }
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-prem")==0){ /* PREM filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->prem_model_filename,argv[i],HC_CHAR_LENGTH);
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-dnp")==0){
hc_toggle_boolean(&p->scale_dens_anom_with_prem);
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-pvel")==0){ /* velocity filename, this will switch off free slip */
hc_advance_argument(&i,argc,argv);
strncpy(p->pvel_filename,argv[i],HC_CHAR_LENGTH);
p->platebc = TRUE;p->no_slip = TRUE;p->free_slip = FALSE;
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-dens")==0){ /* density filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->dens_filename,argv[i],HC_CHAR_LENGTH);
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-dsf")==0){ /* density scaling filename */
p->dd_dens_scale = HC_DD_READ_FROM_FILE;
hc_advance_argument(&i,argc,argv);
strncpy(p->dens_scaling_filename,argv[i],HC_CHAR_LENGTH);
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-dsp")==0){
p->dd_dens_scale = HC_DD_POLYNOMIAL;
hc_advance_argument(&i,argc,argv);
- }else if(strcmp(argv[i],"-visc")==0){ /* viscosity filename */
+ used_parameter = TRUE;
+ }else if(strcmp(argv[i],"-vf")==0){ /* viscosity filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->visc_filename,argv[i],HC_CHAR_LENGTH);
- }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;
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-vdir")==0){ /* velocities in dir */
p->pvel_mode = HC_INIT_P_FROM_DIRECTORY;
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-v")==0){ /* verbosities */
p->verbose = 1;
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-vv")==0){ /* verbosities */
p->verbose = 2;
+ used_parameter = TRUE;
}else if(strcmp(argv[i],"-vvv")==0){
p->verbose = 3;
- }else{
+ used_parameter = TRUE;
+ }
+
+ if(p->solver_mode == HC_SOLVER_MODE_DEFAULT){
+ /*
+ those subsequent options only apply for default solver
+ mode
+ */
+ if(strcmp(argv[i],"-px")==0){ /* print spatial solution? */
+ hc_toggle_boolean(&p->print_spatial);
+ used_parameter = TRUE;
+ }else if(strcmp(argv[i],"-pptsol")==0){ /* print
+ poloidal/toroidal
+ solution
+ parameters */
+ hc_toggle_boolean(&p->print_pt_sol);
+ used_parameter = TRUE;
+ }else if(strcmp(argv[i],"-ng")==0){ /* do not compute geoid */
+ p->compute_geoid = 0;
+ used_parameter = TRUE;
+ }else if(strcmp(argv[i],"-ag")==0){ /* compute geoid at all layers */
+ p->compute_geoid = 2;
+ used_parameter = TRUE;
+ }else if(strcmp(argv[i],"-rg")==0){ /* compute geoid correlations */
+ p->compute_geoid = TRUE;
+ hc_toggle_boolean(&p->compute_geoid_correlations);
+ hc_advance_argument(&i,argc,argv); /* filename */
+ strncpy(p->ref_geoid_file,argv[i],HC_CHAR_LENGTH);
+ hc_read_scalar_shexp(p->ref_geoid_file,&(p->ref_geoid),"reference_geoid",p);
+ used_parameter = TRUE;
+ }else if(strcmp(argv[i],"-rtrac")==0){ /* compute radial
+ tractions */
+ p->solution_mode = HC_RTRACTIONS;
+ used_parameter = TRUE;
+ }else if(strcmp(argv[i],"-htrac")==0){ /* compute horizontal
+ tractions */
+ p->solution_mode = HC_HTRACTIONS;
+ used_parameter = TRUE;
+ }
+ } /* end default operation mode branch */
+ if(!used_parameter){
fprintf(stderr,"%s: can not use parameter %s, use -h for help page\n",
argv[0],argv[i]);
exit(-1);
@@ -599,7 +661,8 @@
hc_get_flt_frmt_string(fstring,2,FALSE);
rold = hc->r_cmb;
/* start read loop */
- while(fscanf(in,HC_TWO_FLT_FORMAT,(hc->rvisc+hc->nvis),(hc->visc+hc->nvis))==2){
+ while(fscanf(in,HC_TWO_FLT_FORMAT,
+ (hc->rvisc+hc->nvis),(hc->visc+hc->nvis))==2){
if(hc->visc[hc->nvis] < 1e15)
fprintf(stderr,"hc_assign_viscosity: WARNING: expecting viscosities in Pas, read %g at layer %i\n",
(double)hc->visc[hc->nvis],hc->nvis);
@@ -715,7 +778,8 @@
if(!hc->orig_danom_saved)
HC_ERROR("hc_assign_density","trying to rescale original density anomaly model, but it was not saved");
for(i=0;i<hc->inho;i++){
- sh_aexp_equals_bexp_coeff((hc->dens_anom+i),(hc->dens_anom_orig+i));
+ sh_aexp_equals_bexp_coeff((hc->dens_anom+i),
+ (hc->dens_anom_orig+i));
}
break;
case HC_INIT_D_FROM_FILE:
@@ -763,9 +827,9 @@
ivec = 0;shps = 1;type = HC_DEFAULT_INTERNAL_FORMAT;
ilayer++;
}else{
- read_on = sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer, &nset,
- &zlabel,&ivec,in,FALSE,density_in_binary,
- verbose);
+ read_on = sh_read_parameters_from_stream(&type,&lmax,&shps,&ilayer, &nset,
+ &zlabel,&ivec,in,FALSE,density_in_binary,
+ verbose);
}
if(read_on){
if((verbose)&&(!reported)){
@@ -781,7 +845,7 @@
}
reported = TRUE;
if(verbose >= 2)
- fprintf(stderr,"hc_assign_density: non_dim radius %% factor PREM \\rho/mean_rho layer # depth[km]\n");
+ fprintf(stderr,"hc_assign_density: non_dim radius %% factor PREM \\rho/mean_rho layer # depth[km] rho[kg/m^3]\n");
}
/*
@@ -823,9 +887,9 @@
/* scaling factor without depth dependence */
dens_scale[0] = HC_DENSITY_SCALING * (HC_PREC)rho0;
if(verbose >= 2){
- fprintf(stderr,"hc_assign_density: r: %11g anom scales: %11g x %11g = %11g\t%5i out of %i, z: %11g\n",
+ fprintf(stderr,"hc_assign_density: r: %11g anom scales: %11g x %11g = %11g\t%5i out of %i, z: %11g %6.1f\n",
(double)hc->rden[hc->inho],
- HC_DENSITY_SCALING,rho0/ (double)hc->avg_den_mantle,(double)dens_scale[0],hc->inho+1,nset,(double)zlabel);
+ HC_DENSITY_SCALING,rho0/ (double)hc->avg_den_mantle,(double)dens_scale[0],hc->inho+1,nset,(double)zlabel,rho0*1000);
}
if(hc->inho){
/*
@@ -860,12 +924,15 @@
will assume input is in physical convention
*/
- sh_read_coefficients_from_file((hc->dens_anom+hc->inho),1,lmax,in,density_in_binary,
- dens_scale,verbose);
+ sh_read_coefficients_from_stream((hc->dens_anom+hc->inho),1,lmax,in,density_in_binary,
+ dens_scale,verbose);
hc->inho++;
} /* end actualy read on */
} /* end while */
-
+ /*
+ assign actual top layer density
+ */
+ hc->rho_top_kg = rho0 * 1000;
if(hc->inho != nset)
HC_ERROR("hc_assign_density","file mode: mismatch in number of layers");
fclose(in);
@@ -968,12 +1035,14 @@
hc_dvecrealloc(&hc->dfact,hc->nrad,"hc_assign_density");
for(i=0;i<hc->nrad;i++){
hc->dfact[i] = 1.0/hc->rden[i] *(dbot[i] - dtop[i]);
+
}
if(verbose)
for(i=0;i < hc->nrad;i++)
- fprintf(stderr,"hc_assign_density: dens %3i: r: %8.6f df: %8.6f |rho|: %8.4f\n",
+ fprintf(stderr,"hc_assign_density: dens %3i: r: %8.6f df: %8.6f |rho|: %8.4f dtop: %8.3f\n",
i+1,(double)hc->rden[i],(double)hc->dfact[i],
- (double)sqrt(sh_total_power((hc->dens_anom+i))));
+ (double)sqrt(sh_total_power((hc->dens_anom+i))),
+ (double)dtop[i]);
free(dbot);free(dtop);
} /* end layer structure part */
@@ -1125,9 +1194,9 @@
ivec = 1;shps = 2;type = HC_DEFAULT_INTERNAL_FORMAT;ilayer=0;zlabel=0;nset=1;
fscanf(in,"%i",&lmax);
}else{
- if(!sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer, &nset,
- &zlabel,&ivec,in,FALSE,
- pvel_in_binary,verbose)){
+ if(!sh_read_parameters_from_stream(&type,&lmax,&shps,&ilayer, &nset,
+ &zlabel,&ivec,in,FALSE,
+ pvel_in_binary,verbose)){
fprintf(stderr,"hc_init_single_plate_exp: read error file %s, expecting long header format\n",
filename);
exit(-1);
@@ -1152,7 +1221,7 @@
read in expansions, convert to internal format from
physical
*/
- sh_read_coefficients_from_file(pvel,shps,-1,in,pvel_in_binary,vfac,verbose);
+ sh_read_coefficients_from_stream(pvel,shps,-1,in,pvel_in_binary,vfac,verbose);
fclose(in);
/*
scale by 1/sqrt(l(l+1))
@@ -1329,36 +1398,39 @@
}
/* read in and assign reference geoid */
-void hc_read_geoid(struct hc_parameters *p)
+void hc_read_scalar_shexp(char *filename,
+ struct sh_lms **sh_exp,
+ char *name,
+ struct hc_parameters *p)
{
int type,lmax,shps,ilayer,nset,ivec;
HC_PREC zlabel;
FILE *in;
HC_PREC fac[3]={1,1,1};
- in = fopen(p->ref_geoid_file,"r");
+ in = fopen(filename,"r");
if(!in){
- fprintf(stderr,"hc_read_geoid: error, could not open ref geoid \"%s\", expecting long scalar SH format\n",
- p->ref_geoid_file);
+ fprintf(stderr,"hc_read_scalar_shexp: error, could not open %s \"%s\", expecting long scalar SH format\n",
+ name,filename);
exit(-1);
}
/* read in the file */
- sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer,&nset,
- &zlabel,&ivec,in,FALSE,
- FALSE,p->verbose);
+ sh_read_parameters_from_stream(&type,&lmax,&shps,&ilayer,&nset,
+ &zlabel,&ivec,in,FALSE,
+ FALSE,p->verbose);
if((ivec != 0)||(shps!=1)){
- fprintf(stderr,"hc_read_geoid: error, expecting scalar, long format SH in %s\n",
- p->ref_geoid_file);
+ fprintf(stderr,"hc_read_scalar_shexp: error, expecting scalar, long format SH in %s\n",
+ filename);
exit(-1);
}
- sh_allocate_and_init(&p->ref_geoid,shps,lmax,type,ivec,p->verbose,0);
- sh_read_coefficients_from_file(p->ref_geoid,shps,-1,in,FALSE,fac,p->verbose);
+ sh_allocate_and_init(sh_exp,shps,lmax,type,ivec,p->verbose,0);
+ sh_read_coefficients_from_stream(*sh_exp,shps,-1,in,FALSE,fac,p->verbose);
fclose(in);
if(p->verbose)
- fprintf(stderr,"hc_read_geoid: read reference geoid from %s, L=%i\n",
- p->ref_geoid_file,p->ref_geoid->lmax);
+ fprintf(stderr,"hc_read_scalar_shexp: read %s from %s, L=%i\n",
+ name,filename,(*sh_exp)->lmax);
}
Deleted: mc/1D/hc/trunk/hc_init.con
===================================================================
--- mc/1D/hc/trunk/hc_init.con 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc_init.con 2013-06-21 23:55:36 UTC (rev 22400)
@@ -1,692 +0,0 @@
-#include "hc.h"
-#include "string.h"
-/*
- general routines dealing with the hager & Connell implementation
-
-
- $Id: hc_init.c,v 1.13 2006/01/22 01:11:34 becker Exp $
-
-*/
-
-/* first, call this routine with a blank **hc */
-void hc_struc_init(struct hcs **hc)
-{
-
- (*hc) = (struct hcs *)calloc(1,sizeof(struct hcs *));
- if(!(*hc))
- HC_MEMERROR("hc_struc_init: (*hc)");
- /*
- assign NULL pointers to allow reallocating
- CO'N: this is not working
- */
- (*hc)->r = (*hc)->visc = (*hc)->rvisc =
- (*hc)->dfact = (*hc)->rden = NULL;
- (*hc)->rpb = (*hc)->fpb= NULL;
- (*hc)->dens_anom = (*hc)->geoid = NULL; /* expansions */
- (*hc)->plm = NULL;
- (*hc)->prem_init = FALSE;
- (*hc)->print_pt_sol = FALSE;
- fprintf(stderr,"struc_init: Did we initialize rvisc ok? Value: %d Address: %d \n",(*hc)->rvisc,&(*hc)->rvisc);
- /*
-
- filenames
-
- */
- /* *hc->visc_filename = malloc(HC_CHAR_LENGTH);
- *hc->dens_filename = (char *)malloc(HC_CHAR_LENGTH);
- *hc->pvel_filename = (char *)malloc(HC_CHAR_LENGTH); */
- strncpy((*hc)->visc_filename,HC_VISC_FILE,HC_CHAR_LENGTH);
- strncpy((*hc)->dens_filename,HC_DENS_SH_FILE,HC_CHAR_LENGTH);
- strncpy((*hc)->pvel_filename,HC_PVEL_FILE,HC_CHAR_LENGTH);
- fprintf(stderr,"Done here in hc_struc_init: visc filename %s and *visc %s \n",(*hc)->visc_filename,(*hc)->visc_filename);
-}
-/*
-
-initialize all variables
-
-sh_type: type of expansion storage/spherical haronics scheme to use
-
-compressible: flag for ddensity factors and polsol operation
-
-vel_bc_zero: if true, will set all surface velocities to zero
-
-
-free-slip: free_slip TRUE
-no-slip : free_slip FALSE vel_bc_zero: TRUE
-platevel : free_slip FALSE vel_bc_zero: FALSE
-
-
-*/
-void hc_init(struct hcs *hc,int sh_type,
- hc_boolean compressible,
- hc_boolean free_slip,
- hc_boolean vel_bc_zero,
- HC_PREC dens_anom_scale,
- hc_boolean verbose)
-{
- int dummy=0;
- /* mechanical boundary condition */
- hc->free_slip = free_slip;
- /*
- set the default expansion type, input expansions will be
- converted
- */
- hc->sh_type = sh_type;
- /*
- start by reading in physical constants and PREM model, if compressible
- */
- hc_init_constants(hc,dens_anom_scale,PREM_MODEL_FILE,verbose);
- /*
- initialize viscosity structure from file
- */
- fprintf(stderr,"Before assign_viscosity: rvisc ok? %d %d \n",&hc->rvisc,hc->rvisc);
- fprintf(stderr,"Before assign_viscosity: visc ok? %d %d \n",&hc->visc,hc->visc);
- hc_assign_viscosity(hc,HC_INIT_FROM_FILE,hc->visc_filename,verbose);
- fprintf(stderr,"past assign visc \n");
- if(vel_bc_zero){
- /* no slip (zero velocity) surface boundary conditions */
- if(free_slip)
- HC_ERROR("hc_init","vel_bc_zero and free_slip doesn't make sense");
- /* read in the densities */
- hc_assign_density(hc,compressible,HC_INIT_FROM_FILE,
- hc->dens_filename,-1,FALSE,FALSE,verbose);
- /* assign all zeroes up to the lmax of the density expansion */
- hc_assign_plate_velocities(hc,HC_INIT_FROM_FILE,
- hc->pvel_filename,vel_bc_zero,
- hc->dens_anom[0].lmax,FALSE,verbose);
- }else{
- /* presribed surface velocities */
- if(!free_slip){
- /* read in velocities, which will determine the solution lmax */
- hc_assign_plate_velocities(hc,HC_INIT_FROM_FILE,
- hc->pvel_filename,vel_bc_zero,
- dummy,FALSE,verbose);
- hc_assign_density(hc,compressible,HC_INIT_FROM_FILE,
- hc->dens_filename,hc->pvel[0].lmax,FALSE,FALSE,verbose);
- }else{
- if(verbose)
- fprintf(stderr,"hc_init: initializing for free-slip\n");
- /* read in the density fields */
- hc_assign_density(hc,compressible,HC_INIT_FROM_FILE,
- hc->dens_filename,-1,FALSE,FALSE,verbose);
- }
- }
- /*
- phase boundaries, if any
- */
- hc_init_phase_boundaries(hc,0,verbose);
- /* */
- hc->save_solution = TRUE; /* (don')t save the propagator
- matrices in hc_polsol and the
- poloidal/toroidal solutions
- */
- hc->initialized = TRUE;
-}
-/*
-
- some of those numbers might be a bit funny, but leave them like
- this for now for backward compatibility.
-
-*/
-void hc_init_constants(struct hcs *hc, HC_PREC dens_anom_scale,
- char *prem_filename,hc_boolean verbose)
-{
- static hc_boolean init=FALSE;
- if(init)
- HC_ERROR("hc_init_constants","why would you call this routine twice?")
- if(!hc->prem_init){
- /* PREM constants */
- prem_read_model(prem_filename,hc->prem,verbose);
- hc->prem_init = TRUE;
- }
- /*
- density scale
- */
- hc->dens_scale[0] = dens_anom_scale;
- /*
- constants
- */
- hc->timesc = HC_TIMESCALE_YR; /* timescale [yr]*/
- hc->visnor = 1e21; /* normalizing viscosity [Pas]*/
- hc->gacc = 10.0; /* gravitational acceleration [m/s2]*/
- hc->g = 6.672e-11; /* gravitational constant [Nm2/kg2]*/
- hc->re = HC_RE_KM*1e3; /* raadius of Earth [m]*/
- hc->secyr = 3.1556926e7; /* seconds/year */
- hc->avg_den_mantle = 4.4488e3;/* average density in mantle [kg/m^3] */
- hc->avg_den_core = 9.90344e3; /* same for core */
- /* velocity scale if input is in [cm/yr],
- works out to be ~0.11 */
- hc->vel_scale = hc->re*PIOVERONEEIGHTY/hc->timesc;
- init = TRUE;
-}
-
-/*
-
- handle command line parameters
-
- visc_filename[] needs to be [HC_CHAR_LENGTH]
-
- */
-void hc_handle_command_line(int argc, char **argv,
- HC_PREC *dens_anom_scale,
- hc_boolean *free_slip,
- char *visc_filename,
- hc_boolean *print_pt_sol,
- hc_boolean *verbose)
-{
- int i;
- for(i=1;i < argc;i++){
- if(strcmp(argv[i],"-h")==0 || strcmp(argv[i],"-?")==0){// help
- /*
- help paghe
- */
- fprintf(stderr,"%s - perform Hager & O'Connell flow computation\n\n",
- argv[0]);
- fprintf(stderr,"options:\n\n");
- fprintf(stderr,"-ds\t\tdensity scaling (%g)\n",
- *dens_anom_scale);
- fprintf(stderr,"-fs\t\tperform free slip computation, else no slip or plates (%s)\n",
- hc_name_boolean(*free_slip));
- fprintf(stderr,"-pptsol\t\tprint pol[6] and tor[2] solution vectors\n");
- fprintf(stderr,"-vf\tname\tset viscosity structure filename to name (%s)\n",
- visc_filename);
- fprintf(stderr,"-v\t-vv\t-vvv: verbosity levels (%i)\n",
- (int)(*verbose));
- fprintf(stderr,"\n\n");
- exit(-1);
- }else if(strcmp(argv[i],"-ds")==0){ /* density anomaly scaling factor */
- hc_advance_argument(&i,argc,argv);
- sscanf(argv[i],HC_FLT_FORMAT,dens_anom_scale);
- }else if(strcmp(argv[i],"-fs")==0){ /* free slip flag */
- hc_toggle_boolean(free_slip);
- }else if(strcmp(argv[i],"-pptsol")==0){ /* print
- poloidal/toroidal
- solution
- parameters */
- hc_toggle_boolean(print_pt_sol);
- }else if(strcmp(argv[i],"-vf")==0){ /* viscosity filename */
- hc_advance_argument(&i,argc,argv);
- fprintf(stderr,"Arg is %s \n",argv[i]);
- strncpy(visc_filename,argv[i],HC_CHAR_LENGTH);
- fprintf(stderr,"Viscosity filename is %s \n",visc_filename);
- }else if(strcmp(argv[i],"-v")==0){ /* verbosities */
- *verbose = 1;
- }else if(strcmp(argv[i],"-vv")==0){ /* verbosities */
- *verbose = 2;
- }else if(strcmp(argv[i],"-vvv")==0){
- *verbose = 3;
- }else{
- fprintf(stderr,"%s: can not use parameter %s, use -h for help page\n",
- argv[0],argv[i]);
- exit(-1);
- }
- }
-}
-/*
-
-assign viscosity structure
-
-mode == 0
-
-read in a viscosity structure with layers of constant viscosity in
-format
-
-r visc
-
-where r is non-dim radius and visc non-dim viscosity, r has to be
-ascending
-
-*/
-void hc_assign_viscosity(struct hcs *hc,int mode,char filename[HC_CHAR_LENGTH],
- hc_boolean verbose)
-{
- FILE *in;
- char fstring[100];
- HC_PREC mean;
- static hc_boolean init=FALSE;
- switch(mode){
- case HC_INIT_FROM_FILE:
- /*
-
- init from file part
-
- */
- if(init)
- HC_ERROR("hc_assign_viscosity","viscosity already read from file, really read again?");
- /*
- read viscosity structure from file
-
- format:
-
- r[non-dim] visc[non-dim]
-
- from bottom to top
- */
- in = hc_open(filename,"r","hc_assign_viscosity");
- fprintf(stderr,"Did we open file ok? %s \n",filename);
- fprintf(stderr,"Whats rvisc before realloc? Ad: %d Val: %d \n",&hc->rvisc,hc->rvisc);
- hc_vecrealloc(&hc->rvisc,1,"hc_assign_viscosity");
- fprintf(stderr,"Whats rvisc after realloc? Ad: %d Val: %g \n",&hc->rvisc,hc->rvisc[0]);
- fprintf(stderr,"Got pas rvisc, trying visc: %d %d \n",&hc->visc,hc->visc);
- hc_vecrealloc(&hc->visc,1,"hc_assign_viscosity");
- fprintf(stderr,"Whats visc after realloc? Ad: %d Val: %g \n",&hc->visc,hc->visc[0]);
- fprintf(stderr,"Past the reallocs \n");
- hc->nvis = 0;
- fprintf(stderr,"Whats up with nvis? %i \n",hc->nvis);
- mean = 0.0;
- /* read sscanf string */
- hc_get_flt_frmt_string(fstring,2,FALSE);
- fprintf(stderr,"Get past get_flt?i fstring %s rvisc: %g nvis: %i visc: %g \n",fstring, hc->rvisc[0], hc->nvis, hc->visc[0]);
- /* start read loop */
- while(fscanf(in,"%lf %lf",
- (hc->rvisc+hc->nvis),(hc->visc+hc->nvis))==2){
- fprintf(stderr,"Made it this far? mean: %g [nvis] %i visc %g %g \n",mean,hc->nvis,hc->visc[0],hc->visc[0]);
- mean += hc->visc[hc->nvis];
- fprintf(stderr,"Made it this far? mean: %g hc->visc: %g \n",mean,hc->visc[hc->nvis]);
- if(hc->nvis){
- if(hc->rvisc[hc->nvis] < hc->rvisc[hc->nvis-1]){
- fprintf(stderr,"hc_assign_viscosity: error: radius has to be ascing, entry %i (%g) smaller than last (%g)\n",
- hc->nvis+1,hc->rvisc[hc->nvis],hc->rvisc[hc->nvis-1]);
- exit(-1);
- }
- }
- hc->nvis++;
- hc_vecrealloc(&hc->rvisc,hc->nvis+1,"hc_assign_viscosity");
- hc_vecrealloc(&hc->visc,hc->nvis+1,"hc_assign_viscosity");
- }
- fprintf(stderr,"Past while loop \n");
- fclose(in);
- mean /= hc->nvis;
- if(verbose){
- fprintf(stderr,"hc_assign_viscosity: read %i layers of non-dimensionalized viscosities from %s\n",
- hc->nvis,filename);
- fprintf(stderr,"hc_assign_viscosity: rough estimate of mean viscosity %g Pas\n",
- mean * hc->visnor);
- }
- fprintf(stderr,"End of init visc \n");
- break;
- default:
- HC_ERROR("hc_assign_viscosity","mode undefined");
- break;
- }
- init = TRUE;
-}
-/*
-
-assign/initialize the density anomalies and density factors
-
-if mode==0: expects spherical harmonics of density anomalies [%] with
- respect to the 1-D reference model (PREM) given in SH
- format on decreasing depth levels in [km]
-
-
- spherical harmonics are real, fully normalized as in
- Dahlen & Tromp p. 859
-
-
-this routine assigns the inho density radii, and the total (nrad=inho)+2
-radii
-
-furthermore, the dfact factors are assigned as well
-
-set density_in_binary to TRUE, if expansion given in binary
-
-nominal_lmax: -1: the max order of the density expansion will either
- determine the lmax of the solution (free-slip, or vel_bc_zero) or
- will have to be the same as the plate expansion lmax (!free_slip && !vel_bc_zero)
- else: will zero out all entries > nominal_lmax
-
-*/
-void hc_assign_density(struct hcs *hc,
- hc_boolean compressible,int mode,
- char *filename,int nominal_lmax,
- hc_boolean layer_structure_changed,
- hc_boolean density_in_binary,
- hc_boolean verbose)
-{
- FILE *in;
- int type,lmax,shps,ilayer,nset,ivec,i,j;
- HC_PREC *dtop,*dbot,zlabel,dens_scale[1],rho0;
- hc_boolean reported = FALSE;
- static HC_PREC local_dens_fac = .01; /* this factor will be multiplied with
- the hc->dens_fac factor to arrive at
- fractional anomalies from input. set to
- 0.01 for percent input, for instance
- */
- static hc_boolean init=FALSE;
- hc->compressible = compressible;
-
- if(init) /* clear old expansions, if
- already initialized */
- sh_free_expansion(hc->dens_anom,hc->inho);
- /* get PREM model, if not initialized */
- if(!hc->prem_init)
- HC_ERROR("hc_assign_density","assign 1-D reference model (PREM) first");
- switch(mode){
- case HC_INIT_FROM_FILE:
- if(init)
- HC_ERROR("hc_assign_density","really read dens anomalies again from file?");
-
- /*
-
- read in density anomalies in spherical harmonics format for
- different layers from file.
-
- this assumes that we are reading in anomalies in percent
-
- */
- in = hc_open(filename,"r","hc_assign_density");
- if(verbose)
- fprintf(stderr,"hc_assign_density: reading density anomalies in [%%] from %s, scaling by %g\n",
- filename,hc->dens_scale[0]);
- hc->inho = 0; /* counter for density layers */
- hc->dens_anom = (struct sh_lms *)
- realloc(hc->dens_anom,sizeof(struct sh_lms));
- if(!hc->dens_anom)
- HC_MEMERROR("hc_assign_density: dens anom");
- /*
- read all layes as spherical harmonics assuming real Dahlen & Tromp
- (physical) normalization
-
- */
- while(sh_read_parameters(&type,&lmax,&shps,&ilayer, &nset,
- &zlabel,&ivec,in,FALSE,density_in_binary,
- verbose)){
- if((verbose)&&(!reported)){
- if(nominal_lmax > lmax)
- fprintf(stderr,"hc_assign_density: density lmax: %3i filling up to nominal lmax: %3i with zeroes\n",
- lmax,nominal_lmax);
- if(nominal_lmax != -1){
- fprintf(stderr,"hc_assign_density: density lmax: %3i limiting to lmax: %3i\n",
- lmax,nominal_lmax);
- }else{
- fprintf(stderr,"hc_assign_density: density lmax: %3i determines solution lmax\n",
- lmax);
- }
- reported = TRUE;
- }
- /*
- do tests
- */
- if((shps != 1)||(ivec))
- HC_ERROR("hc_assign_density","vector field read in but only scalar expansion expected");
- /* test and assign depth levels */
- hc->rden=(HC_PREC *)
- realloc(hc->rden,(1+hc->inho)*sizeof(HC_PREC));
- if(!hc->rden)
- HC_MEMERROR("hc_assign_density: rden");
- /*
- assign depth, this assumes that we are reading in depths [km]
- */
- hc->rden[hc->inho] = HC_ND_RADIUS(zlabel);
- /*
-
- get reference density at this level
-
- */
- prem_get_rho(&rho0,hc->rden[hc->inho],hc->prem);
- /*
- general density (add additional depth dependence here)
- */
- dens_scale[0] = hc->dens_scale[0] * local_dens_fac * rho0;
- if(verbose >= 2)
- fprintf(stderr,"hc_assign_density: r: %11g anom scales: %11g x %11g x %11g = %11g\n",
- hc->rden[hc->inho],hc->dens_scale[0],
- local_dens_fac,rho0,dens_scale[0]);
- if(hc->inho){
- /*
- check by comparison with previous expansion
- */
- if(nominal_lmax == -1)
- if(lmax != hc->dens_anom[0].lmax)
- HC_ERROR("hc_assign_density","lmax changed in file");
- if(hc->rden[hc->inho] <= hc->rden[hc->inho-1])
- HC_ERROR("hc_assign_density","depth should decrease, radius increase (give z[km])");
- }
- /*
- make room for new expansion
- */
- hc->dens_anom = (struct sh_lms *)
- realloc(hc->dens_anom,(1+hc->inho)*sizeof(struct sh_lms));
- if(!hc->dens_anom)
- HC_MEMERROR("hc_assign_density");
- /*
- initialize expansion
- */
- sh_init_expansion((hc->dens_anom+hc->inho),(nominal_lmax > lmax) ? (nominal_lmax):(lmax),
- hc->sh_type,1,verbose);
- /*
-
- read parameters and scale (put possible depth dependence of
- scaling here)
-
- will assume input is in physical convention
-
- */
- sh_read_coefficients((hc->dens_anom+hc->inho),1,lmax,
- in,density_in_binary,dens_scale,
- verbose);
- hc->inho++;
- }
- if(hc->inho != nset)
- HC_ERROR("hc_assign_density","file mode: mismatch in number of layers");
- fclose(in);
- break;
- default:
- HC_ERROR("hc_assign_density","mode undefined");
- break;
- }
- if((!init)||(layer_structure_changed)){
- /*
-
- assign the other radii, nrad + 2
-
- */
- hc->nrad = hc->inho;
- hc_vecrealloc(&hc->r,(2+hc->nrad),"hc_assign_density");
- hc->r[0] = HC_RCMB_ND; /* CMB */
- if(hc->rden[0] <= hc->r[0])
- HC_ERROR("hc_assign_density","first density layer has to be above internal CMD limit");
- for(i=0;i<hc->nrad;i++) /* density layers */
- hc->r[i+1] = hc->rden[i];
- if(hc->rden[hc->nrad-1] >= 1.0)
- HC_ERROR("hc_assign_density","uppermost density layer has to be below surface");
- hc->r[hc->nrad+1] = 1.0; /* surface */
- /*
-
- assign the density jump factors
-
- */
- /*
- since we have spherical harmonics at several layers, we assign
- the layer thickness by picking the intermediate depths
- */
- hc_vecalloc(&dbot,hc->nrad,"hc_assign_density");
- hc_vecalloc(&dtop,hc->nrad,"hc_assign_density");
- // top boundaries
- j = hc->nrad-1;
- for(i=0;i < j;i++)
- dtop[i] = 1.0 - (hc->rden[i+1] + hc->rden[i])/2.0;
- dtop[j] = 0.0; // top boundary
- // bottom boundaries
- dbot[0] = 1.0 - HC_RCMB_ND; // bottom boundary, ie. CMB
- for(i=1;i < hc->nrad;i++)
- dbot[i] = dtop[i-1];
- /*
- density layer thickness factors
- */
- hc_dvecrealloc(&hc->dfact,hc->nrad,"hc_assign_density");
- for(i=0;i<hc->nrad;i++){
- hc->dfact[i] = 1.0/hc->rden[i] *(dbot[i] - dtop[i]);
- }
- if(verbose)
- for(i=0;i < hc->nrad;i++)
- fprintf(stderr,"hc_assign_density: dens %3i: r: %8.6f df: %8.6f |rho|: %8.4f\n",
- i+1,hc->rden[i],hc->dfact[i],
- sqrt(sh_total_power((hc->dens_anom+i))));
- free(dbot);free(dtop);
- } /* end layer structure part */
- init = TRUE;
-}
-/*
-
-assign phase boundary jumps
-input:
-
-npb: number of phase boundaries
-
-....
-
-
-
-*/
-void hc_init_phase_boundaries(struct hcs *hc, int npb,
- hc_boolean verbose)
-{
-
- hc->npb = npb; /* no phase boundaries for now */
- if(hc->npb){
- HC_ERROR("hc_init_phase_boundaries","phase boundaries not implemented yet");
- hc_vecrealloc(&hc->rpb,hc->npb,"hc_init_phase_boundaries");
- hc_vecrealloc(&hc->fpb,hc->npb,"hc_init_phase_boundaries");
- }
-
-}
-
-/*
-
-read in plate velocities,
-
-vel_bc_zero: if true, will set all surface velocities to zero
-
-lmax will only be referenced if all velocities are supposed to be set to zero
-
-we expect velocities to be in cm/yr, convert to m/yr
-
-*/
-
-void hc_assign_plate_velocities(struct hcs *hc,int mode,
- char *filename,
- hc_boolean vel_bc_zero,int lmax,
- hc_boolean pvel_in_binary,
- hc_boolean verbose)
-{
- static hc_boolean init = FALSE;
- int type,shps,ilayer,nset,ivec;
- HC_PREC zlabel,vfac[2],t10[2],t11[2];
- FILE *in;
- /* scale to go from cm/yr to internal scale */
- vfac[0] = vfac[1] = hc->vel_scale;
- if(init)
- HC_ERROR("hc_assign_plate_velocities","what to do if called twice?");
- if(!vel_bc_zero){
- /*
-
- velocities are NOT all zero
-
-
- */
- switch(mode){
- case HC_INIT_FROM_FILE:
- /*
- read velocities in pol/tor expansion format from file
- in units of HC_VELOCITY_FILE_FACTOR per year
- */
- if(verbose)
- fprintf(stderr,"hc_assign_plate_velocities: expecting [cm/yr] pol/tor from %s\n",
- filename);
- in = hc_open(filename,"r","hc_assign_plate_velocities");
- if(!sh_read_parameters(&type,&lmax,&shps,&ilayer, &nset,
- &zlabel,&ivec,in,FALSE,
- pvel_in_binary,verbose)){
- fprintf(stderr,"hc_assign_plate_velocities: read error file %s\n",
- filename);
- exit(-1);
- } /* check if we read in two sets of expansions */
- if(shps != 2){
- fprintf(stderr,"hc_assign_plate_velocities: two sets expected but found shps: %i in file %s\n",
- shps,filename);
- exit(-1);
- }
- if((nset > 1)||(fabs(zlabel) > 0.01)){
- fprintf(stderr,"hc_assign_plate_velocities: error: expected one layer at surface, but nset: %i z: %g\n",
- nset, zlabel);
- exit(-1);
- }
- /*
- initialize expansion
- */
- sh_init_expansion((hc->pvel+0),lmax,hc->sh_type,1,verbose);
- sh_init_expansion((hc->pvel+1),lmax,hc->sh_type,1,verbose);
- /*
- read in expansions, convert to internal format from
- physical
- */
- sh_read_coefficients(hc->pvel,shps,-1,in,pvel_in_binary,
- vfac,verbose);
- fclose(in);
- /*
- scale by 1/sqrt(l(l+1))
- */
- if(hc->pvel[0].lmax > hc->lfac_init)
- hc_init_l_factors(hc,hc->pvel[0].lmax);
- sh_scale_expansion_l_factor((hc->pvel+0),hc->ilfac);
- sh_scale_expansion_l_factor((hc->pvel+1),hc->ilfac);
- /*
- check for net rotation
- */
- sh_get_coeff((hc->pvel+1),1,0,0,TRUE,t10);
- sh_get_coeff((hc->pvel+1),1,0,2,TRUE,t11);
- if(fabs(t10[0])+fabs(t11[0])+fabs(t11[1]) > 1.0e-7)
- fprintf(stderr,"\nhc_assign_plate_velocities: WARNING: toroidal A(1,0): %g A(1,1): %g B(1,1): %g\n\n",
- t10[0],t11[0],t11[1]);
- if(verbose)
- fprintf(stderr,"hc_assign_plate_velocities: read velocities, lmax %i: |pol|: %11g |tor|: %11g\n",
- lmax,sqrt(sh_total_power((hc->pvel+0))),sqrt(sh_total_power((hc->pvel+1))));
- break;
- default:
- HC_ERROR("hc_assign_plate_velocities","op mode undefined");
- }
- }else{
- /*
- initialize with zeroes
- */
- if(init){
- sh_clear_alm(hc->pvel);
- sh_clear_alm((hc->pvel+1));
- }else{
- sh_init_expansion(hc->pvel,lmax,hc->sh_type,
- 1,verbose);
- sh_init_expansion((hc->pvel+1),lmax,hc->sh_type,
- 1,verbose);
- }
- if(verbose)
- fprintf(stderr,"hc_assign_plate_velocities: using no-slip surface BC, lmax %i\n",
- lmax);
- }
- init = TRUE;
-}
-
-/*
-
-initialize an array with sqrt(l(l+1)) factors
-from l=0 .. lmax+1
-
-pass lfac initialized (say, as NULL)
-
-*/
-void hc_init_l_factors(struct hcs *hc, int lmax)
-{
- int lmaxp1,l;
- lmaxp1 = lmax + 1;
- hc_vecrealloc(&hc->lfac,lmaxp1,"hc_init_l_factors");
- hc_vecrealloc(&hc->ilfac,lmaxp1,"hc_init_l_factors");
- /* maybe optimize later */
- hc->lfac[0] = 0.0;
- hc->ilfac[0] = 1.0; /* shouldn't matter */
- for(l=1;l < lmaxp1;l++){
- hc->lfac[l] = sqrt((HC_PREC)l * ((HC_PREC)l + 1.0));
- hc->ilfac[l] = 1.0/hc->lfac[l];
- }
- hc->lfac_init = lmax;
-}
Modified: mc/1D/hc/trunk/hc_input.c
===================================================================
--- mc/1D/hc/trunk/hc_input.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc_input.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -23,9 +23,9 @@
*/
n = os = 0;
- while(sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer,
- &nset,&zlabel,&ivec,in,
- FALSE,binary,verbose)){
+ while(sh_read_parameters_from_stream(&type,&lmax,&shps,&ilayer,
+ &nset,&zlabel,&ivec,in,
+ FALSE,binary,verbose)){
hc->sh_type = type;
if(ilayer != n){
fprintf(stderr,"hc_read_sh_solution: error: ilayer %i n %i\n",
@@ -50,8 +50,8 @@
/*
read coefficients
*/
- sh_read_coefficients_from_file((*sol+os),shps,-1,in,
- binary,unity,verbose);
+ sh_read_coefficients_from_stream((*sol+os),shps,-1,in,
+ binary,unity,verbose);
if(verbose){
if(shps == 3)
fprintf(stderr,"hc_read_sh_solution: z: %8.3f |r|: %12.5e |pol|: %12.5e |tor|: %12.5e\n",
Added: mc/1D/hc/trunk/hc_invert_dtopo.c
===================================================================
--- mc/1D/hc/trunk/hc_invert_dtopo.c (rev 0)
+++ mc/1D/hc/trunk/hc_invert_dtopo.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -0,0 +1,128 @@
+#include "hc.h"
+/*
+
+ implementation of Hager & O'Connell (1981) method of solving mantle
+ circulation given internal density anomalies, only radially varying
+ viscosity, and either free-slip or plate velocity boundary
+ condition at surface. based on Hager & O'Connell (1981), Hager &
+ Clayton (1989), and Steinberger (2000). the original code is due to
+ Brad Hager, Rick O'Connell, and was largely modified by Bernhard
+ Steinberger. this version by Thorsten Becker (twb at usc.edu) for
+ additional comments, see hc.c
+
+ invert for compositional anomalies given geoid anomalies [m] and
+ residual topography wrt. to air. [m]
+
+*/
+
+int main(int argc, char **argv)
+{
+ struct hcs *model; /* main structure, make sure to initialize with
+ zeroes */
+ struct sh_lms *sol_spectral=NULL, *geoid = NULL, *dtopo = NULL; /* solution expansions */
+ struct sh_lms *pvel=NULL; /* local plate velocity expansion */
+ int nsol,lmax,solved;
+ struct hc_parameters p[1]; /* parameters */
+ HC_PREC gcorr[3],dcorr[3]; /* correlations */
+ hc_struc_init(&model);
+ hc_init_parameters(p);
+
+ /*
+
+ special options for this computation
+
+ */
+ p->solver_mode = HC_SOLVER_MODE_DYNTOPO_INVERT;
+ p->compute_geoid = 1;
+ p->compute_geoid_correlations = TRUE;
+ p->solution_mode = HC_RTRACTIONS; /* make sure to compute tractions */
+ /* */
+ p->verbose = 1;
+ if(argc > 2){
+ /* read in the reference geoid */
+ strcpy(p->ref_geoid_file,argv[1]);
+ hc_read_scalar_shexp(p->ref_geoid_file,&(p->ref_geoid),
+ "reference geoid",p);
+ /* read in the reference topography */
+ strcpy(p->ref_dtopo_file,argv[2]);
+ hc_read_scalar_shexp(p->ref_dtopo_file,&(p->ref_dtopo),
+ "reference dynamic topography",p);
+ }else{
+ fprintf(stderr,"%s: ERROR: need geoid.ab dtopo.ab file as arguments\n",argv[0]);
+ fprintf(stderr,"%s: usage:\n\n%s geoid.ab dtopo.ab\n\n",argv[0],argv[0]);
+ fprintf(stderr,"%s: for help, use\n\n%s geoid.ab dtopo.ab -h\n\n",argv[0],argv[0]);
+ exit(-1);
+ }
+
+ /*
+ handle other command line arguments
+ */
+ hc_handle_command_line(argc,argv,3,p);
+ /*
+
+ begin main program part
+
+ */
+ hc_init_main(model,SH_RICK,p);
+ nsol = (model->nradp2) * 3;
+ if(p->free_slip) /* maximum degree is determined by the
+ density expansion */
+ lmax = model->dens_anom[0].lmax;
+ else /* max degree is determined by the
+ plate velocities */
+ lmax = model->pvel.p[0].lmax; /* shouldn't be larger than that*/
+
+ sh_allocate_and_init(&pvel,2,lmax,model->sh_type,1,p->verbose,FALSE);
+ sh_allocate_and_init(&sol_spectral,nsol,lmax,model->sh_type,HC_VECTOR,
+ p->verbose,FALSE);
+ sh_allocate_and_init(&geoid,1,model->dens_anom[0].lmax,
+ model->sh_type,HC_SCALAR,p->verbose,FALSE);
+ if(!p->free_slip)
+ hc_select_pvel(p->pvel_time,&model->pvel,pvel,p->verbose);
+
+ /*
+
+
+ */
+ solved=0;
+ {
+ /* compute solution */
+ hc_solve(model,p->free_slip,p->solution_mode,sol_spectral,
+ (solved)?(FALSE):(TRUE), /* density changed? */
+ (solved)?(FALSE):(TRUE), /* plate velocity changed? */
+ TRUE, /* viscosity changed */
+ FALSE,p->compute_geoid,
+ pvel,model->dens_anom,geoid,
+ p->verbose);
+ /* extract the top tractions */
+ hc_compute_dynamic_topography(model,sol_spectral,&dtopo,FALSE,p->verbose);
+
+ sh_single_par_and_exp_to_file(dtopo,"dtopo.ab",TRUE,p->verbose);
+
+ /* geoid correlation */
+ hc_compute_correlation(geoid,p->ref_geoid,(gcorr),0,p->verbose); /* full correlation */
+ hc_compute_correlation(geoid,p->ref_geoid,(gcorr+1),1,p->verbose); /* up to 20 and 4...9 */
+ fprintf(stdout,"geoid full: %10.7f L=20: %10.7f \n",(double)gcorr[0],(double)gcorr[1]);
+
+ hc_compute_correlation(dtopo,p->ref_dtopo,(dcorr),0,p->verbose); /* full correlation */
+ hc_compute_correlation(dtopo,p->ref_dtopo,(dcorr+1),1,p->verbose); /* up to 20 and 4..9 */
+ fprintf(stdout,"dtopo full: %10.7f L=20: %10.7f \n",(double)dcorr[0],(double)dcorr[1]);
+
+ solved++;
+ }
+ /*
+
+ free memory
+
+ */
+ sh_free_expansion(sol_spectral,nsol);
+ /* local copies of plate velocities */
+ sh_free_expansion(pvel,2);
+ /* */
+ sh_free_expansion(geoid,1);
+ if(p->verbose)
+ fprintf(stderr,"%s: done\n",argv[0]);
+ hc_struc_free(&model);
+ return 0;
+}
+
Modified: mc/1D/hc/trunk/hc_misc.c
===================================================================
--- mc/1D/hc/trunk/hc_misc.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc_misc.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -213,18 +213,26 @@
*i += 1;
}
+/*
+ compute the correlation between two scalar fields
+ mode 0 : full correlation
+ 1 : up to 20 and between 4...9
+
+ */
void hc_compute_correlation(struct sh_lms *g1,struct sh_lms *g2,
HC_PREC *c,int mode,hc_boolean verbose)
{
int lmaxg;
lmaxg = MIN(g1->lmax,g1->lmax);
- lmaxg = MIN(20,lmaxg);
+
switch(mode){
- case 0:
+ case 0: /* 1...LMAX */
if(verbose)
fprintf(stderr,"hc_compute_correlation: computing 1...%i\n",lmaxg);
c[0] = sh_correlation_per_degree(g1,g2,1,lmaxg);
+ break;
case 1: /* 1...20 and 4..9 correlations */
+ lmaxg = MIN(20,lmaxg);
if(verbose)
fprintf(stderr,"hc_compute_correlation: computing 1...%i and 4..9 correlations\n",lmaxg);
c[0] = sh_correlation_per_degree(g1,g2,1,lmaxg);
Modified: mc/1D/hc/trunk/hc_output.c
===================================================================
--- mc/1D/hc/trunk/hc_output.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc_output.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -39,15 +39,15 @@
/*
write parameters, convert radius to depth in [km]
*/
- sh_print_parameters_to_file((sol+os),ntype,i,hc->nradp2,
- HC_Z_DEPTH(hc->r[i]),
- out,FALSE,binary,verbose);
+ sh_print_parameters_to_stream((sol+os),ntype,i,hc->nradp2,
+ HC_Z_DEPTH(hc->r[i]),
+ out,FALSE,binary,verbose);
/*
write the set of coefficients in D&T convention
*/
- sh_print_coefficients_to_file((sol+os),ntype,out,fac,
+ sh_print_coefficients_to_stream((sol+os),ntype,out,fac,
binary,verbose);
if(verbose >= 2){
switch(sol_mode){
@@ -99,9 +99,9 @@
hc_boolean verbose)
{
HC_CPREC fac[1] = {1.0};
- sh_print_parameters_to_file(sh,1,0,1,0.0,out,
- short_format,binary,verbose); /* parameters in long format */
- sh_print_coefficients_to_file(sh,1,out,fac,binary,verbose); /* coefficients */
+ sh_print_parameters_to_stream(sh,1,0,1,0.0,out,
+ short_format,binary,verbose); /* parameters in long format */
+ sh_print_coefficients_to_stream(sh,1,out,fac,binary,verbose); /* coefficients */
}
@@ -781,8 +781,8 @@
sh_c_is_a_plus_b_coeff((exp+2),(exp+0),(exp+1)); /* c = a+b */
/* print to file */
- sh_print_parameters_to_file((exp+2),1,i,hc->nradp2,HC_Z_DEPTH(hc->r[i]),out,FALSE,binary,verbose);
- sh_print_coefficients_to_file((exp+2),1,out,fac,binary,verbose);
+ sh_print_parameters_to_stream((exp+2),1,i,hc->nradp2,HC_Z_DEPTH(hc->r[i]),out,FALSE,binary,verbose);
+ sh_print_coefficients_to_stream((exp+2),1,out,fac,binary,verbose);
if(verbose>2)fprintf(stderr,"hc_print_dens_anom: z: %8.3f (f1: %6.3f f2: %6.3f) %3i/%3i pow: %10.3e %10.3e %10.3e\n",
(double)HC_Z_DEPTH(hc->r[i]),
(double)f1,(double)f2,i+1,hc->nradp2,
Modified: mc/1D/hc/trunk/hc_polsol.c
===================================================================
--- mc/1D/hc/trunk/hc_polsol.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc_polsol.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -77,7 +77,10 @@
anomalies changes between
call, but nothing else
*/
- hc_boolean verbose)
+ hc_boolean verbose, /* output options */
+ hc_boolean calc_kernel_only /* only compute the
+ kernels */
+ )
{
// ****************************************************************
@@ -193,7 +196,8 @@
int i,i2,i3,i6,j,l,m,nih,nxtv,ivis,os,pos1,pos2,gi,g1,g2,gic,
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,jsol;
+ nprops_max,jsol,mmax;
+ int klayer = 1;
double *xprem;
HC_HIGH_PREC *b,du1,du2,el,rnext,drho,dadd;
HC_PREC rbound_kludge;
@@ -560,8 +564,8 @@
begin m loop
*/
-
- for(m=0;m <= l;m++){
+ mmax = (calc_kernel_only)?(0):(l);
+ for(m=0;m <= mmax;m++){
/*
@@ -574,7 +578,7 @@
//
a_or_b = 0; /* start with A coefficient */
do{ /* do loop for A/B evaluation */
- if(l <= dens_anom[0].lmax){
+ if((!calc_kernel_only) && (l <= dens_anom[0].lmax)){
/*
obtain the coefficients from the density field expansions
@@ -594,6 +598,9 @@
b[i] = 0.0;
}
b[inho] = 0.0;
+ if(calc_kernel_only)
+ b[klayer] = 1.0;
+
//
// U(C) = [0,VC,SC,0], U(A) = [0,0,SA,SX]
// POT(A) = [U5(A),-(L+1)*U5(A)]T, POT(C) = [U5(C),L*U5(C)]T
Modified: mc/1D/hc/trunk/hc_solve.c
===================================================================
--- mc/1D/hc/trunk/hc_solve.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/hc_solve.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -97,7 +97,7 @@
hc->npb,hc->rpb,hc->fpb,free_slip,
(pvel+0),hc->pol_sol,
compute_geoid,geoid,hc->save_solution,
- verbose);
+ verbose,FALSE);
if(print_pt_sol) /* print poloidal solution without the
scaling factors */
hc_print_poloidal_solution(hc->pol_sol,hc,31, /* print only up
@@ -315,3 +315,45 @@
}
hc->spatial_solution_computed = TRUE;
}
+
+/*
+ calculate dynamic topgoraphy given radial tractions in [MPa]
+
+ pass dtopo as NULL initialized
+
+
+*/
+void hc_compute_dynamic_topography(struct hcs *hc,struct sh_lms *spectral_sol,
+ struct sh_lms **dtopo,
+ hc_boolean scale_from_MPa_to_m,
+ hc_boolean verbose)
+{
+ HC_PREC scale;
+ const int shps = 3; /* radial component of stress */
+ int nlayer;
+ nlayer = hc->nradp2-1; /* top layer */
+ /* original solution is non-dim */
+ scale = hc->stress_scale/hc->r[nlayer]; /* go to MPa */
+ if(scale_from_MPa_to_m){
+ /*
+ output will be in [m]
+ */
+ scale *= -1./(hc->rho_top_kg*(HC_GACC/100))*1e6;
+ }
+
+ if(verbose){
+ if(scale_from_MPa_to_m)
+ fprintf(stderr,"hc_compute_dynamic_topography: density %g to scale stress [MPa] to [m] with %g (g: %g) layer %i\n",(double)hc->rho_top_kg,(double)scale,(double)HC_GACC/100,hc->nradp2);
+ else
+ fprintf(stderr,"hc_compute_dynamic_topography: leaving in MPa, layer %i\n",
+ hc->nradp2);
+ }
+ /* create a new expansion */
+ sh_allocate_and_init(dtopo,1,
+ spectral_sol[nlayer*shps].lmax,
+ spectral_sol[nlayer*shps].type,
+ FALSE, verbose,FALSE);
+ /* assign */
+ sh_copy_lms((spectral_sol+nlayer*shps),*dtopo);
+ sh_scale_expansion(*dtopo,scale); /* scale */
+}
Added: mc/1D/hc/trunk/hc_visc_scan.c
===================================================================
--- mc/1D/hc/trunk/hc_visc_scan.c (rev 0)
+++ mc/1D/hc/trunk/hc_visc_scan.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -0,0 +1,208 @@
+#include "hc.h"
+/*
+
+ implementation of Hager & O'Connell (1981) method of solving mantle
+ circulation given internal density anomalies, only radially varying
+ viscosity, and either free-slip or plate velocity boundary
+ condition at surface. based on Hager & O'Connell (1981), Hager &
+ Clayton (1989), and Steinberger (2000). the original code is due to
+ Brad Hager, Rick O'Connell, and was largely modified by Bernhard
+ Steinberger. this version by Thorsten Becker (twb at usc.edu) for
+ additional comments, see hc.c
+
+ scan through viscosities and compute correlation with the geoid
+
+*/
+
+int main(int argc, char **argv)
+{
+ struct hcs *model; /* main structure, make sure to initialize with
+ zeroes */
+ struct sh_lms *sol_spectral=NULL, *geoid = NULL; /* solution expansions */
+ struct sh_lms *pvel=NULL; /* local plate velocity expansion */
+ int nsol,lmax,solved;
+ struct hc_parameters p[1]; /* parameters */
+ HC_PREC corr[2]; /* correlations */
+ HC_PREC vl[4][3],v[4],dv; /* for viscosity scans */
+ /*
+
+
+ (1)
+
+ initialize the model structure, this is needed to initialize some
+ of the default values before callign the parameter handling
+ routine this call also involves initializing the hc parameter
+ structure
+
+ */
+ hc_struc_init(&model);
+ /*
+
+ (2)
+ init parameters to default values
+
+ */
+ hc_init_parameters(p);
+ /*
+
+ special options for this computation
+
+ */
+ p->solver_mode = HC_SOLVER_MODE_VISC_SCAN;
+ p->visc_init_mode = HC_INIT_E_FOUR_LAYERS;
+ p->compute_geoid = 1;
+ p->compute_geoid_correlations = TRUE;
+
+ if(argc > 1){
+ /* read in the reference geoid */
+ strcpy(p->ref_geoid_file,argv[1]);
+ hc_read_scalar_shexp(p->ref_geoid_file,&(p->ref_geoid),"reference geoid",p);
+ }else{
+ fprintf(stderr,"%s: ERROR: need geoid.ab file as an argument\n",argv[0]);
+ fprintf(stderr,"%s: usage:\n\n%s geoid.ab\n\n",argv[0],argv[0]);
+ fprintf(stderr,"%s: for help, use:\n\n%s geoid.ab -h\n\n",argv[0],argv[0]);
+ exit(-1);
+ }
+
+ /*
+ handle other command line arguments
+ */
+ hc_handle_command_line(argc,argv,2,p);
+ /*
+
+ begin main program part
+
+ */
+#ifdef __TIMESTAMP__
+ if(p->verbose)
+ fprintf(stderr,"%s: starting version compiled on %s\n",
+ argv[0],__TIMESTAMP__);
+#else
+ if(p->verbose)
+ fprintf(stderr,"%s: starting main program\n",argv[0]);
+#endif
+ /*
+
+ (3)
+
+ initialize all variables
+
+ - choose the internal spherical harmonics convention
+ - assign constants
+ - assign phase boundaries, if any
+ - read in viscosity structure
+ - assign density anomalies
+ - read in plate velocities
+
+ */
+ hc_init_main(model,SH_RICK,p);
+ nsol = (model->nradp2) * 3; /*
+ number of solutions (r,pol,tor) * (nlayer+2)
+
+ total number of layers is nlayer +2,
+
+ because CMB and surface are added
+ to intermediate layers which are
+ determined by the spacing of the
+ density model
+
+ */
+ if(p->free_slip) /* maximum degree is determined by the
+ density expansion */
+ lmax = model->dens_anom[0].lmax;
+ else /* max degree is determined by the
+ plate velocities */
+ lmax = model->pvel.p[0].lmax; /* shouldn't be larger than that*/
+ /*
+ make sure we have room for the plate velocities
+ */
+ sh_allocate_and_init(&pvel,2,lmax,model->sh_type,1,p->verbose,FALSE);
+
+ /* init done */
+ /*
+
+
+
+ SOLUTION PART
+
+
+ */
+ /*
+ make room for the spectral solution on irregular grid
+ */
+ sh_allocate_and_init(&sol_spectral,nsol,lmax,model->sh_type,HC_VECTOR,
+ p->verbose,FALSE);
+ /* make room for geoid solution at surface */
+ sh_allocate_and_init(&geoid,1,model->dens_anom[0].lmax,
+ model->sh_type,HC_SCALAR,p->verbose,FALSE);
+
+ /* parameter space log bounds */
+
+ dv = .1; /* spacing */
+
+ vl[0][0]= -3;vl[0][1]=3+1e-5;vl[0][2]=dv; /* 0..100 layer log bounds and spacing */
+ vl[1][0]= -3;vl[1][1]=3+1e-5;vl[1][2]=dv; /* 100..410 */
+ if(p->free_slip){
+ vl[2][0]= 0;vl[2][1]=0+1e-5;vl[2][2]=dv; /* for free slip,
+ only relative
+ viscosisites
+ matter for
+ correlation */
+ }else{
+ vl[2][0]= -3;vl[2][1]=3+1e-5;vl[2][2]=dv; /* need to actually
+ loop 410 .660 */
+ }
+ vl[3][0]= -3;vl[3][1]=3+1e-5;vl[3][2]=dv; /* 660 ... 2871 */
+
+
+ /* */
+ /* select plate velocity */
+ if(!p->free_slip)
+ hc_select_pvel(p->pvel_time,&model->pvel,pvel,p->verbose);
+
+ solved=0;
+ for(v[0]=vl[0][0];v[0] <= vl[0][1];v[0] += vl[0][2])
+ for(v[1]=vl[1][0];v[1] <= vl[1][1];v[1] += vl[1][2])
+ for(v[2]=vl[2][0];v[2] <= vl[2][1];v[2] += vl[2][2])
+ for(v[3]=vl[3][0];v[3] <= vl[3][1];v[3] += vl[3][2]){
+ /* layer viscosity structure */
+ p->elayer[0] = pow(10,v[3]); /* bottom */
+ p->elayer[1] = pow(10,v[2]); /* 660..410 */
+ p->elayer[2] = pow(10,v[1]); /* 410..100 */
+ p->elayer[3] = pow(10,v[0]); /* 100..0 */
+ hc_assign_viscosity(model,HC_INIT_E_FOUR_LAYERS,p->elayer,p);
+ /* print viscosities of 0...100, 100...410, 410 ... 660 and 660...2871 layer */
+ fprintf(stdout,"%14.7e %14.7e %14.7e %14.7e\t",
+ (double)p->elayer[3],(double)p->elayer[2],
+ (double)p->elayer[1],(double)p->elayer[0]);
+ /* compute solution */
+ hc_solve(model,p->free_slip,p->solution_mode,sol_spectral,
+ (solved)?(FALSE):(TRUE), /* density changed? */
+ (solved)?(FALSE):(TRUE), /* plate velocity changed? */
+ TRUE, /* viscosity changed */
+ FALSE,p->compute_geoid,
+ pvel,model->dens_anom,geoid,
+ p->verbose);
+ /* only output are the geoid correlations, for now */
+ hc_compute_correlation(geoid,p->ref_geoid,corr,1,p->verbose);
+ fprintf(stdout,"%10.7f %10.7f ",
+ (double)corr[0],(double)corr[1]);
+ fprintf(stdout,"\n");
+ solved++;
+ }
+ /*
+
+ free memory
+
+ */
+ sh_free_expansion(sol_spectral,nsol);
+ /* local copies of plate velocities */
+ sh_free_expansion(pvel,2);
+ /* */
+ sh_free_expansion(geoid,1);
+ if(p->verbose)
+ fprintf(stderr,"%s: done\n",argv[0]);
+ hc_struc_free(&model);
+ return 0;
+}
+
Modified: mc/1D/hc/trunk/itg-hc-geoid.ab
===================================================================
--- mc/1D/hc/trunk/itg-hc-geoid.ab 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/itg-hc-geoid.ab 2013-06-21 23:55:36 UTC (rev 22400)
@@ -2,14 +2,14 @@
0 0
0 0
0 0
- -1.0046197e+02 0.0000000e+00
+ -1.1552017e+02 0.0000000e+00
6.0024606e-03 -3.3358521e-02
5.5154257e+01 -3.1660072e+01
2.1641332e+01 0.0000000e+00
-4.5908573e+01 -5.6117892e+00
2.0457177e+01 -1.3995660e+01
-1.6309025e+01 -3.1978316e+01
- -1.0341030e+01 0.0000000e+00
+ -1.0084792e+01 0.0000000e+00
1.2122473e+01 1.0707316e+01
7.9248115e+00 1.4978616e+01
-2.2403185e+01 4.5436135e+00
Deleted: mc/1D/hc/trunk/main.c
===================================================================
--- mc/1D/hc/trunk/main.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/main.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -1,359 +0,0 @@
-#include "hc.h"
-/*
-
-
-implementation of Hager & O'Connell (1981) method of solving mantle
-circulation given internal density anomalies, only radially varying
-viscosity, and either free-slip or plate velocity boundary condition
-at surface
-
-based on Hager & O'Connell (1981), Hager & Clayton (1989), and
-Steinberger (2000)
-
-the original code is due to Brad Hager, Rick O'Connell, and was
-largely modified by Bernhard Steinberger
-
-this version by Thorsten Becker (twb at usc.edu)
-
-$Id: main.c,v 1.13 2006/01/22 01:11:34 becker Exp becker $
-
->>> SOME COMMENTS FROM THE ORIGINAL CODE <<<
-
-C * It uses the following Numerical Recipes (Press et al.) routines:
-C four1, realft, gauleg, rk4, rkdumb, ludcmp, lubksb;
-C !!!!!!!!!!!!!!!!!!! rkqc, odeint !!!!!!!!!! take out !!!!!!
-C and the following routines by R.J. O'Connell:
-C shc2d, shd2c, ab2cs, cs2ab, plmbar, plvel2sh, pltgrid, pltvel,
-C vshd2c, plmbar1.
-C Further subroutines are: kinsub, evalpa,
-C torsol (all based on "kinflow" by Hager & O'Connell),
-C densub and evppot (based on "denflow" by Hager & O'Connell),
-C sumsub (based on "sumkd" by Hager & O'Connell, but
-C substantially speeded up),
-C convert, derivs and shc2dd (based on R.J. O'Connell's shc2d).
-C
-C bugs found:
-C * The combination of (1) high viscosity lithosphere
-C (2) compressible flow (3) kinematic (plate driven) flow
-C doesn't work properly. The problem presumably only occurs
-C at degree 1 (I didn't make this sure) but this is sufficient
-C to screw up everything. It will usually work to reduce the
-C contrast between lithospheric and asthenospheric viscosity.
-C Then make sure that (1) two viscosity structures give similar
-C results for incompressible and (2) incompressible and compressible
-C reduced viscosity look similar (e.g. anomalous mass flux vs. depth)
-
-<<< END OLD COMMENTS
-
-*/
-
-int main(int argc, char **argv)
-{
- struct hcs *model; /* main structure, make sure to initialize with
- zeroes */
- struct sh_lms *sol_spectral=NULL, *geoid = NULL; /* solution expansions */
- struct sh_lms *pvel=NULL; /* local plate velocity expansion */
- int nsol,lmax,solved,i;
- FILE *out;
- struct hc_parameters p[1]; /* parameters */
- char filename[HC_CHAR_LENGTH],file_prefix[10];
- HC_PREC *sol_spatial = NULL; /* spatial solution,
- e.g. velocities */
- HC_PREC corr[2]; /* correlations */
- HC_PREC vl[4][3],v[4],dv; /* for viscosity scans */
- static hc_boolean geoid_binary = FALSE; /* type of geoid output */
- static HC_CPREC unitya[1] = {1.0};
- /*
-
-
- (1)
-
- initialize the model structure, this is needed to initialize some
- of the default values before callign the parameter handling
- routine this call also involves initializing the hc parameter
- structure
-
- */
- hc_struc_init(&model);
- /*
-
- (2)
- init parameters to default values
-
- */
- hc_init_parameters(p);
- /*
- handle command line arguments
- */
- hc_handle_command_line(argc,argv,p);
- /*
-
- begin main program part
-
- */
-#ifdef __TIMESTAMP__
- if(p->verbose)
- fprintf(stderr,"%s: starting version compiled on %s\n",argv[0],__TIMESTAMP__);
-#else
- if(p->verbose)
- fprintf(stderr,"%s: starting main program\n",argv[0]);
-#endif
- /*
-
- (3)
-
- initialize all variables
-
- - choose the internal spherical harmonics convention
- - assign constants
- - assign phase boundaries, if any
- - read in viscosity structure
- - assign density anomalies
- - read in plate velocities
-
- */
- hc_init_main(model,SH_RICK,p);
- nsol = (model->nradp2) * 3; /*
- number of solutions (r,pol,tor) * (nlayer+2)
-
- total number of layers is nlayer +2,
-
- because CMB and surface are added
- to intermediate layers which are
- determined by the spacing of the
- density model
-
- */
- if(p->free_slip) /* maximum degree is determined by the
- density expansion */
- lmax = model->dens_anom[0].lmax;
- else /* max degree is determined by the
- plate velocities */
- lmax = model->pvel.p[0].lmax; /* shouldn't be larger than that*/
- /*
- make sure we have room for the plate velocities
- */
- sh_allocate_and_init(&pvel,2,lmax,model->sh_type,1,p->verbose,FALSE);
-
- /* init done */
- /*
-
-
-
- SOLUTION PART
-
-
- */
- /*
- make room for the spectral solution on irregular grid
- */
- sh_allocate_and_init(&sol_spectral,nsol,lmax,model->sh_type,HC_VECTOR,
- p->verbose,FALSE);
- if(p->compute_geoid == 1)
- /* make room for geoid solution at surface */
- sh_allocate_and_init(&geoid,1,model->dens_anom[0].lmax,
- model->sh_type,HC_SCALAR,p->verbose,FALSE);
- else if(p->compute_geoid == 2) /* all layers */
- sh_allocate_and_init(&geoid,model->nradp2,model->dens_anom[0].lmax,
- model->sh_type,HC_SCALAR,p->verbose,FALSE);
-
- switch(p->solver_mode){
- case HC_SOLVER_MODE_DEFAULT:
- /*
- solve poloidal and toroidal part and sum
- */
- if(!p->free_slip)
- hc_select_pvel(p->pvel_time,&model->pvel,pvel,p->verbose);
- hc_solve(model,p->free_slip,p->solution_mode,sol_spectral,
- TRUE,TRUE,TRUE,p->print_pt_sol,p->compute_geoid,
- pvel,model->dens_anom,geoid,
- p->verbose);
- /*
-
- OUTPUT PART
-
- */
- /*
-
- output of spherical harmonics solution
-
- */
- switch(p->solution_mode){
- case HC_VEL:
- sprintf(file_prefix,"vel");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;
- }
- if(p->sol_binary_out)
- sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_BINARY);
- else
- sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_ASCII);
- if(p->verbose)
- fprintf(stderr,"%s: writing spherical harmonics solution to %s\n",
- argv[0],filename);
- out = ggrd_open(filename,"w","main");
- hc_print_spectral_solution(model,sol_spectral,out,
- p->solution_mode,
- p->sol_binary_out,p->verbose);
- fclose(out);
- /* */
- if(p->print_density_field){
- /*
- print the density field
- */
- sprintf(file_prefix,"dscaled");
- if(p->sol_binary_out)
- sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_BINARY);
- else
- sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_ASCII);
- if(p->verbose)
- fprintf(stderr,"%s: writing scaled density anomaly field to %s\n",
- argv[0],filename);
-
- out = ggrd_open(filename,"w","main");
- hc_print_dens_anom(model,out,p->sol_binary_out,p->verbose);
- fclose(out);
- }
-
-
- /* compute the geoid? */
- if(p->compute_geoid){
- if(p->compute_geoid_correlations){
- if(p->compute_geoid == 2){ /* check if all geoids were computed */
- fprintf(stderr,"%s: can only compute correlation for surface geoid, geoid = %i\n",
- argv[0],p->compute_geoid);
- exit(-1);
- }
- if(p->verbose)
- fprintf(stderr,"%s: correlation for geoid with %s\n",argv[0],
- p->ref_geoid_file);
- hc_compute_correlation(geoid,p->ref_geoid,corr,1,p->verbose);
- fprintf(stdout,"%10.7f %10.7f\n",(double)corr[0],(double)corr[1]);
- }else{
- /*
- print geoid solution
- */
- sprintf(filename,"%s",HC_GEOID_FILE);
- if(p->verbose)
- fprintf(stderr,"%s: writing geoid to %s, %s\n",
- argv[0],filename,(p->compute_geoid == 1)?("at surface"):("all layers"));
- out = ggrd_open(filename,"w","main");
- if(p->compute_geoid == 1) /* surface layer */
- hc_print_sh_scalar_field(geoid,out,FALSE,geoid_binary,p->verbose);
- else{ /* all layers */
- for(i=0;i < model->nradp2;i++){
- sh_print_parameters_to_file((geoid+i),1,i,model->nradp2,
- HC_Z_DEPTH(model->r[i]),out,FALSE,geoid_binary,p->verbose);
- sh_print_coefficients_to_file((geoid+i),1,out,unitya,geoid_binary,p->verbose);
- }
- }
- fclose(out);
- }
- }
- if(p->print_spatial){
- /*
- we wish to use the spatial solution
-
- expand velocities to spatial base, compute spatial
- representation
-
- */
- hc_compute_sol_spatial(model,sol_spectral,&sol_spatial,
- p->verbose);
- /*
- output of spatial solution
- */
- sprintf(filename,"%s.%s",file_prefix,HC_SPATIAL_SOLOUT_FILE);
- /* print lon lat z v_r v_theta v_phi */
- hc_print_spatial_solution(model,sol_spectral,sol_spatial,
- filename,HC_LAYER_OUT_FILE,
- p->solution_mode,p->sol_binary_out,
- p->verbose);
- }
- break; /* end default mode */
- case HC_SOLVER_MODE_VISC_SCAN: /* scan through viscosities for a
- four layer problem */
-
-
- /* parameter space log bounds */
-
- dv = .1; /* spacing */
-
- vl[0][0]= -3;vl[0][1]=3+1e-5;vl[0][2]=dv; /* 0..100 layer log bounds and spacing */
- vl[1][0]= -3;vl[1][1]=3+1e-5;vl[1][2]=dv; /* 100..410 */
- if(p->free_slip){
- vl[2][0]= 0;vl[2][1]=0+1e-5;vl[2][2]=dv; /* for free slip,
- only relative
- viscosisites
- matter for
- correlation */
- }else{
- vl[2][0]= -3;vl[2][1]=3+1e-5;vl[2][2]=dv; /* need to actually
- loop 410 .660 */
- }
- vl[3][0]= -3;vl[3][1]=3+1e-5;vl[3][2]=dv; /* 660 ... 2871 */
-
-
- /* */
- /* select plate velocity */
- if(!p->free_slip)
- hc_select_pvel(p->pvel_time,&model->pvel,pvel,p->verbose);
-
- solved=0;
- for(v[0]=vl[0][0];v[0] <= vl[0][1];v[0] += vl[0][2])
- for(v[1]=vl[1][0];v[1] <= vl[1][1];v[1] += vl[1][2])
- for(v[2]=vl[2][0];v[2] <= vl[2][1];v[2] += vl[2][2])
- for(v[3]=vl[3][0];v[3] <= vl[3][1];v[3] += vl[3][2]){
- /* layer viscosity structure */
- p->elayer[0] = pow(10,v[3]); /* bottom */
- p->elayer[1] = pow(10,v[2]); /* 660..410 */
- p->elayer[2] = pow(10,v[1]); /* 410..100 */
- p->elayer[3] = pow(10,v[0]); /* 100..0 */
- hc_assign_viscosity(model,HC_INIT_E_FOUR_LAYERS,p->elayer,p);
- /* print viscosities of 0...100, 100...410, 410 ... 660 and 660...2871 layer */
- fprintf(stdout,"%14.7e %14.7e %14.7e %14.7e\t",
- (double)p->elayer[3],(double)p->elayer[2],
- (double)p->elayer[1],(double)p->elayer[0]);
- /* compute solution */
- hc_solve(model,p->free_slip,p->solution_mode,sol_spectral,
- (solved)?(FALSE):(TRUE), /* density changed? */
- (solved)?(FALSE):(TRUE), /* plate velocity changed? */
- TRUE, /* viscosity changed */
- FALSE,p->compute_geoid,
- pvel,model->dens_anom,geoid,
- p->verbose);
- /* only output are the geoid correlations, for now */
- hc_compute_correlation(geoid,p->ref_geoid,corr,1,p->verbose);
- fprintf(stdout,"%10.7f %10.7f ",(double)corr[0],(double)corr[1]);
- fprintf(stdout,"\n");
- solved++;
- }
- break;
- default:
- HC_ERROR("hc","solver mode undefined");
- }
-
- /*
-
- free memory
-
- */
- sh_free_expansion(sol_spectral,nsol);
- /* local copies of plate velocities */
- sh_free_expansion(pvel,2);
- /* */
- if(p->compute_geoid == 1) /* only one layer */
- sh_free_expansion(geoid,1);
- else if(p->compute_geoid == 2) /* all layers */
- sh_free_expansion(geoid,model->nradp2);
- free(sol_spatial);
- if(p->verbose)
- fprintf(stderr,"%s: done\n",argv[0]);
- hc_struc_free(&model);
- return 0;
-}
-
Modified: mc/1D/hc/trunk/sh_ana.c
===================================================================
--- mc/1D/hc/trunk/sh_ana.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/sh_ana.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -164,7 +164,7 @@
sh_read_spatial_data_from_grd(exp,ggrd,use_3d,shps,data,&zlabel);
}else{
/* read in data from stdin */
- sh_read_spatial_data_from_file(exp,stdin,use_3d,shps,data,&zlabel);
+ sh_read_spatial_data_from_stream(exp,stdin,use_3d,shps,data,&zlabel);
}
/*
perform spherical harmonic expansion
@@ -172,12 +172,12 @@
sh_compute_spectral(data,ivec,FALSE,&dummy,
exp,verbose);
/* print parameters of expansion */
- sh_print_parameters_to_file(exp,shps,ilayer,nset,zlabel,
- stdout,short_format,binary,
- verbose);
+ sh_print_parameters_to_stream(exp,shps,ilayer,nset,zlabel,
+ stdout,short_format,binary,
+ verbose);
/* print coefficients */
- sh_print_coefficients_to_file(exp,shps,stdout,fac,binary,
- verbose);
+ sh_print_coefficients_to_stream(exp,shps,stdout,fac,binary,
+ verbose);
}
fprintf(stderr,"%s: printing to stdout, done\n",argv[0]);
free(exp);free(data);
Modified: mc/1D/hc/trunk/sh_corr.c
===================================================================
--- mc/1D/hc/trunk/sh_corr.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/sh_corr.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -81,9 +81,9 @@
argv[0],(short_format)?("(short format)"):("(long format)"),argv[0]);
/* allocate two expansions */
i=0;
- while(sh_read_parameters_from_file(&type,(lmax+i),&shps,&ilayer,&nset,
- &zlabel,&ivec,stdin,short_format,
- binary,verbose)){
+ while(sh_read_parameters_from_stream(&type,(lmax+i),&shps,&ilayer,&nset,
+ &zlabel,&ivec,stdin,short_format,
+ binary,verbose)){
if(ivec){
fprintf(stderr,"%s: error, vector fields not implemented yet\n",argv[0]);
exit(-1);
@@ -94,10 +94,10 @@
/* input and init */
if(i==0){
sh_allocate_and_init(&exp1,shps,lmax[i],type,ivec,verbose,0);
- sh_read_coefficients_from_file(exp1,shps,-1,stdin,binary,fac,verbose);
+ sh_read_coefficients_from_stream(exp1,shps,-1,stdin,binary,fac,verbose);
}else{
sh_allocate_and_init(&exp2,shps,lmax[i],type,ivec,verbose,0);
- sh_read_coefficients_from_file(exp2,shps,-1,stdin,binary,fac,verbose);
+ sh_read_coefficients_from_stream(exp2,shps,-1,stdin,binary,fac,verbose);
}
if(verbose)
fprintf(stderr,"ok.\n");
Modified: mc/1D/hc/trunk/sh_exp.c
===================================================================
--- mc/1D/hc/trunk/sh_exp.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/sh_exp.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -308,10 +308,33 @@
tmp = sqrt(sum[1]*sum[2]);
return sum[0]/tmp;
}
+void sh_single_par_and_exp_to_file(struct sh_lms *exp, char *name,
+ hc_boolean binary,hc_boolean verbose)
+{
+ FILE *out;
+ out = fopen(name,"w");
+ if(!out){
+ fprintf(stderr,"sh_single_par_and_exp_to_file: ERROR: problem openeing %s\n",name);
+ exit(-1);
+ }
+ sh_single_par_and_exp_to_stream(exp,out,binary,verbose);
+ fclose(out);
+ if(verbose)
+ fprintf(stderr,"sh_single_par_and_exp_to_file: written to %s\n",name);
+}
+void sh_single_par_and_exp_to_stream(struct sh_lms *exp, FILE *out,
+ hc_boolean binary,hc_boolean verbose)
+{
+ HC_PREC fac[1]={1.0};
+ const hc_boolean short_format = FALSE;
+ sh_print_parameters_to_stream(exp,1,0,1,0,out,short_format,binary,verbose);
+ sh_print_coefficients_to_stream(exp,1,out,fac,binary,verbose);
+}
+
/*
print one line with all parameters needed to identify a spherical
@@ -330,10 +353,10 @@
lmax
*/
-void sh_print_parameters_to_file(struct sh_lms *exp, int shps,
- int ilayer, int nset,HC_CPREC zlabel,
- FILE *out, hc_boolean short_format,
- hc_boolean binary,hc_boolean verbose)
+void sh_print_parameters_to_stream(struct sh_lms *exp, int shps,
+ int ilayer, int nset,HC_CPREC zlabel,
+ FILE *out, hc_boolean short_format,
+ hc_boolean binary,hc_boolean verbose)
{
HC_PREC fz;
/*
@@ -412,15 +435,15 @@
*/
-hc_boolean sh_read_parameters_from_file(int *type, int *lmax,
- int *shps,
- int *ilayer, int *nset,
- HC_CPREC *zlabel,
- int *ivec,
- FILE *in,
- hc_boolean short_format,
- hc_boolean binary,
- hc_boolean verbose)
+hc_boolean sh_read_parameters_from_stream(int *type, int *lmax,
+ int *shps,
+ int *ilayer, int *nset,
+ HC_CPREC *zlabel,
+ int *ivec,
+ FILE *in,
+ hc_boolean short_format,
+ hc_boolean binary,
+ hc_boolean verbose)
{
int input1[2],input2[3];
HC_PREC fz;
@@ -536,11 +559,11 @@
fac[3] scales the coefficients
*/
-void sh_print_coefficients_to_file(struct sh_lms *exp,
- int shps, FILE *out,
- HC_CPREC *fac,
- hc_boolean binary,
- hc_boolean verbose)
+void sh_print_coefficients_to_stream(struct sh_lms *exp,
+ int shps, FILE *out,
+ HC_CPREC *fac,
+ hc_boolean binary,
+ hc_boolean verbose)
{
int j,l,m;
HC_PREC value[2];
@@ -587,7 +610,7 @@
fprintf(out,"\n");
} /* end m loop */
} /* end l loop */
- fprintf(out,"\n");
+ //fprintf(out,"\n");
}
}
/*
@@ -605,9 +628,9 @@
*/
-void sh_read_coefficients_from_file(struct sh_lms *exp, int shps, int lmax,
- FILE *in, hc_boolean binary, HC_CPREC *fac,
- hc_boolean verbose)
+void sh_read_coefficients_from_stream(struct sh_lms *exp, int shps, int lmax,
+ FILE *in, hc_boolean binary, HC_CPREC *fac,
+ hc_boolean verbose)
{
int j,k,l,m,lmax_loc;
HC_CPREC value[2]={0,0};
@@ -621,12 +644,12 @@
*/
for(j=1;j < shps;j++){ /* check the lmax */
if(exp[j].lmax != exp[0].lmax){
- fprintf(stderr,"sh_read_coefficients_from_file: error: lmax(%i):%i != lmax(0):%i\n",
+ fprintf(stderr,"sh_read_coefficients_from_stream: error: lmax(%i):%i != lmax(0):%i\n",
j+1,exp[j].lmax,exp[0].lmax);
exit(-1);
}
if(exp[j].type != exp[0].type ){
- fprintf(stderr,"sh_read_coefficients_from_file: error: type(%i):%i != type(0):%i\n",
+ fprintf(stderr,"sh_read_coefficients_from_stream: error: type(%i):%i != type(0):%i\n",
j+1,exp[j].type,exp[0].type);
exit(-1);
}
@@ -636,7 +659,7 @@
for(m=0;m <= l;m++)
for(j=0;j < shps;j++){
if(hc_read_float(fvalue,2,in)!=2){
- fprintf(stderr,"sh_read_coefficients_from_file: read error: set %i l %i m %i\n",
+ fprintf(stderr,"sh_read_coefficients_from_stream: read error: set %i l %i m %i\n",
j+1,l,m);
exit(-1);
}
@@ -651,7 +674,7 @@
for(m=0;m <= l;m++)
for(j=0;j < shps;j++){
if(fscanf(in,HC_TWO_FLT_FORMAT,value,(value+1))!=2){
- fprintf(stderr,"sh_read_coefficients_from_file: read error: set %i l %i m %i, last val: %g %g\n",
+ fprintf(stderr,"sh_read_coefficients_from_stream: read error: set %i l %i m %i, last val: %g %g\n",
j+1,l,m,(double)value[0],(double)value[1]);
exit(-1);
}
@@ -713,9 +736,9 @@
instead
*/
-void sh_read_spatial_data_from_file(struct sh_lms *exp, FILE *in,
+void sh_read_spatial_data_from_stream(struct sh_lms *exp, FILE *in,
my_boolean use_3d,
- int shps, HC_PREC *data,
+ int shps, HC_PREC *data,
HC_PREC *z)
{
struct ggrd_gt *gdummy=(struct ggrd_gt *)NULL;
@@ -1327,9 +1350,9 @@
shps is the number of scalars that are passed in the data[shps * npoints] array
*/
-void sh_print_spatial_data_to_file(struct sh_lms *exp, int shps,
- HC_PREC *data, hc_boolean use_3d,
- HC_PREC z, FILE *out)
+void sh_print_spatial_data_to_stream(struct sh_lms *exp, int shps,
+ HC_PREC *data, hc_boolean use_3d,
+ HC_PREC z, FILE *out)
{
int j,k;
HC_PREC lon,lat;
@@ -1404,11 +1427,11 @@
regular grid version
*/
-void sh_print_reg_spatial_data_to_file(struct sh_lms *exp, int shps,
- HC_PREC *data, hc_boolean use_3d,
- HC_PREC z, HC_PREC *theta,int ntheta,
- HC_PREC *phi,int nphi,
- FILE *out)
+void sh_print_reg_spatial_data_to_stream(struct sh_lms *exp, int shps,
+ HC_PREC *data, hc_boolean use_3d,
+ HC_PREC z, HC_PREC *theta,int ntheta,
+ HC_PREC *phi,int nphi,
+ FILE *out)
{
int i,j,k,l,npoints;
HC_PREC lon,lat;
@@ -1420,7 +1443,7 @@
#ifdef HC_USE_HEALPIX
case SH_HEALPIX:
- HC_ERROR("sh_print_reg_spatial_data_to_file","healpix not implemented");
+ HC_ERROR("sh_print_reg_spatial_data_to_stream","healpix not implemented");
break;
#endif
case SH_RICK:
@@ -1446,7 +1469,7 @@
#ifdef HC_USE_SPHEREPACK
case SH_SPHEREPACK_GAUSS:
case SH_SPHEREPACK_EVEN:
- HC_ERROR("sh_print_reg_spatial_data_to_file","spherepack not implemented");
+ HC_ERROR("sh_print_reg_spatial_data_to_stream","spherepack not implemented");
break;
#endif
default:
@@ -1459,10 +1482,10 @@
irregular, arbitrary version
*/
-void sh_print_irreg_spatial_data_to_file(struct sh_lms *exp, int shps,
- HC_PREC *data, hc_boolean use_3d,
- HC_PREC z, HC_PREC *theta,HC_PREC *phi,int npoints,
- FILE *out)
+void sh_print_irreg_spatial_data_to_stream(struct sh_lms *exp, int shps,
+ HC_PREC *data, hc_boolean use_3d,
+ HC_PREC z, HC_PREC *theta,HC_PREC *phi,int npoints,
+ FILE *out)
{
int i,k;
HC_PREC lon,lat;
@@ -1472,7 +1495,7 @@
switch(exp->type){
#ifdef HC_USE_HEALPIX
case SH_HEALPIX:
- HC_ERROR("sh_print_irreg_spatial_data_to_file","healpix not implemented");
+ HC_ERROR("sh_print_irreg_spatial_data_to_stream","healpix not implemented");
break;
#endif
case SH_RICK:
@@ -1495,7 +1518,7 @@
#ifdef HC_USE_SPHEREPACK
case SH_SPHEREPACK_GAUSS:
case SH_SPHEREPACK_EVEN:
- HC_ERROR("sh_print_irreg_spatial_data_to_file","spherepack not implemented");
+ HC_ERROR("sh_print_irreg_spatial_data_to_stream","spherepack not implemented");
break;
#endif
default:
@@ -1782,7 +1805,7 @@
break;
#endif
default:
- sh_exp_type_error("sh_read_coefficients_from_file",exp);
+ sh_exp_type_error("sh_read_coefficients_from_stream",exp);
break;
}
}
Modified: mc/1D/hc/trunk/sh_extract_layer.c
===================================================================
--- mc/1D/hc/trunk/sh_extract_layer.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/sh_extract_layer.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -56,8 +56,8 @@
in = ggrd_open(argv[1],"r","sh_extract_layer");
/* start loop */
i = 0;
- sh_read_parameters_from_file(&type,&lmax, &shps,&i,&nset,(zdepth+i),
- &ivec,in,short_format_in,binary,verbose);
+ sh_read_parameters_from_stream(&type,&lmax, &shps,&i,&nset,(zdepth+i),
+ &ivec,in,short_format_in,binary,verbose);
if(verbose)
fprintf(stderr,"sh_extract_layer: detected %i layers, vec: %i, type %i SH file, shps: %i, layer %i at depth %g\n",
nset,ivec,type,shps,i+1,(double)zdepth[i]);
@@ -76,10 +76,10 @@
}
for(;i<nset;i++){
if(i != 0)
- sh_read_parameters_from_file(&type,&lmax,&shps,&i,&nset,(zdepth+i),
- &ivec,in,short_format_in,binary,verbose);
- sh_read_coefficients_from_file((exp+i),shps,-1,in,binary,unitya,
- verbose);
+ sh_read_parameters_from_stream(&type,&lmax,&shps,&i,&nset,(zdepth+i),
+ &ivec,in,short_format_in,binary,verbose);
+ sh_read_coefficients_from_stream((exp+i),shps,-1,in,binary,unitya,
+ verbose);
if(i == ilayer){ /* output */
if(verbose)
fprintf(stderr,"%s: printing SH from %s at layer %i out of %i to stdout (depth: %g)\n",
@@ -87,9 +87,9 @@
/*
output will remove the layer number information
*/
- sh_print_parameters_to_file((exp+i),shps,ilayer,1,zdepth[i],
- stdout,short_format_output,FALSE,verbose);
- sh_print_coefficients_to_file((exp+i),shps,stdout,unitya,FALSE,verbose);
+ sh_print_parameters_to_stream((exp+i),shps,ilayer,1,zdepth[i],
+ stdout,short_format_output,FALSE,verbose);
+ sh_print_coefficients_to_stream((exp+i),shps,stdout,unitya,FALSE,verbose);
}
}
fclose(in);
Modified: mc/1D/hc/trunk/sh_model.c
===================================================================
--- mc/1D/hc/trunk/sh_model.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/sh_model.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -167,12 +167,12 @@
HC_ERROR("sh_print_model_coefficients","coefficients not initialized expansion");
for(i = 0; i < model->nset;i++){ /* loop through sets */
/* output of parameters */
- sh_print_parameters_to_file((model->exp+i*model->shps),
+ sh_print_parameters_to_stream((model->exp+i*model->shps),
model->shps,i,model->nset,
model->z[i],out,FALSE,binary,verbose);
/* output of coefficient(s) (1 (ivec=0) or 3 (ivec=1)) */
- sh_print_coefficients_to_file((model->exp+i*model->shps),
- model->shps,out,fac,binary,verbose);
+ sh_print_coefficients_to_stream((model->exp+i*model->shps),
+ model->shps,out,fac,binary,verbose);
} /* end set loop */
free(fac);
}
@@ -261,11 +261,11 @@
}
for(i=0;i < model->nset;i++)
/* read in the data for this layer */
- sh_read_spatial_data_from_file((model->exp+i*model->shps),
- in,(model->nset!=1)?(TRUE):(FALSE),
- model->shps,
- (*data+model->shps*model->exp[i*model->shps].npoints),
- (model->z+i));
+ sh_read_spatial_data_from_stream((model->exp+i*model->shps),
+ in,(model->nset!=1)?(TRUE):(FALSE),
+ model->shps,
+ (*data+model->shps*model->exp[i*model->shps].npoints),
+ (model->z+i));
model->spatial_init = TRUE;
}
@@ -372,10 +372,10 @@
os = 0;
for(i=0;i < model->nset;i++){
/* print out data for the set */
- sh_print_spatial_data_to_file((model->exp+i*model->shps),
- (model->ivec)?(3):(1),(data+os),
- (model->nset == 1)?(FALSE):(TRUE),
- model->z[i],out);
+ sh_print_spatial_data_to_stream((model->exp+i*model->shps),
+ (model->ivec)?(3):(1),(data+os),
+ (model->nset == 1)?(FALSE):(TRUE),
+ model->z[i],out);
if(model->ivec)
for(j=0;j<3;j++)
os += model->exp[i*model->shps+j].npoints;
Modified: mc/1D/hc/trunk/sh_power.c
===================================================================
--- mc/1D/hc/trunk/sh_power.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/sh_power.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -65,16 +65,16 @@
}
fprintf(stderr,"%s: awaiting spherical harmonics expansion (%s) from stdin (use %s -h for help)\n",
argv[0],short_format ? "short format" : "long format",argv[0]);
- while(sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer,&nset,
- &zlabel,&ivec,stdin,short_format,
- binary,verbose)){
+ while(sh_read_parameters_from_stream(&type,&lmax,&shps,&ilayer,&nset,
+ &zlabel,&ivec,stdin,short_format,
+ binary,verbose)){
fprintf(stderr,"%s: computing power per degree and unit area for lmax %i ivec: %i at z: %g\n",
argv[0],lmax,ivec,(double)zlabel);
/*
input and init
*/
sh_allocate_and_init(&exp,shps,lmax, type,ivec,verbose,FALSE);
- sh_read_coefficients_from_file(exp,shps,-1,stdin,binary,fac,verbose);
+ sh_read_coefficients_from_stream(exp,shps,-1,stdin,binary,fac,verbose);
/* get space */
hc_vecrealloc(&power,exp->lmaxp1 * shps,"sh_power");
/* compute the powers */
Modified: mc/1D/hc/trunk/sh_syn.c
===================================================================
--- mc/1D/hc/trunk/sh_syn.c 2013-06-21 18:00:48 UTC (rev 22399)
+++ mc/1D/hc/trunk/sh_syn.c 2013-06-21 23:55:36 UTC (rev 22400)
@@ -85,9 +85,9 @@
if(verbose)
fprintf(stderr,"%s: waiting to read spherical harmonic coefficients from stdin (use %s -h for help)\n",
argv[0],argv[0]);
- while(sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer,&nset,
- &zlabel,&ivec,stdin,short_format,
- binary,verbose)){
+ while(sh_read_parameters_from_stream(&type,&lmax,&shps,&ilayer,&nset,
+ &zlabel,&ivec,stdin,short_format,
+ binary,verbose)){
if(short_format_ivec){
ivec = 1;
shps = 2;
@@ -98,7 +98,7 @@
/* input and init */
sh_allocate_and_init(&exp,shps,lmax,type,ivec,verbose,((regular_basis != 0)?(1):(0)));
- sh_read_coefficients_from_file(exp,shps,-1,stdin,binary,fac,verbose);
+ sh_read_coefficients_from_stream(exp,shps,-1,stdin,binary,fac,verbose);
if(regular_basis == 1){
/*
regular basis output on regular grid
@@ -137,10 +137,10 @@
theta,ntheta,phi,nphi,data,
verbose,FALSE);
/* output */
- sh_print_reg_spatial_data_to_file(exp,shps,data,
- (nset>1)?(TRUE):(FALSE),
- zlabel, theta,ntheta,
- phi,nphi,stdout);
+ sh_print_reg_spatial_data_to_stream(exp,shps,data,
+ (nset>1)?(TRUE):(FALSE),
+ zlabel, theta,ntheta,
+ phi,nphi,stdout);
}else if(regular_basis == -1){
/* output on locations input lon lat file */
if(verbose)
@@ -165,15 +165,15 @@
fclose(in);
hc_vecalloc(&data,npoints * shps,"sh_shsyn data");
sh_compute_spatial_irreg(exp,ivec,theta,phi,npoints,data,verbose);
- sh_print_irreg_spatial_data_to_file(exp,shps,data,(nset>1)?(TRUE):(FALSE),
- zlabel,theta,phi,npoints,stdout);
+ sh_print_irreg_spatial_data_to_stream(exp,shps,data,(nset>1)?(TRUE):(FALSE),
+ zlabel,theta,phi,npoints,stdout);
}else{ /* use the built in spatial basis (Gaussian) */
/* expansion */
hc_vecalloc(&data,exp[0].npoints * shps,"sh_syn");
sh_compute_spatial(exp,ivec,FALSE,&dummy,data,verbose);
/* output */
- sh_print_spatial_data_to_file(exp,shps,data,(nset>1)?(TRUE):(FALSE),
- zlabel,stdout);
+ sh_print_spatial_data_to_stream(exp,shps,data,(nset>1)?(TRUE):(FALSE),
+ zlabel,stdout);
}
free(exp);free(data);
}
Added: mc/1D/hc/trunk/visc.dat
===================================================================
--- mc/1D/hc/trunk/visc.dat (rev 0)
+++ mc/1D/hc/trunk/visc.dat 2013-06-21 23:55:36 UTC (rev 22400)
@@ -0,0 +1,4 @@
+0.546225 5e22
+0.8964 1e21
+0.9356 1e20
+0.9843 5e22
More information about the CIG-COMMITS
mailing list