forked from geodynamics/hc
-
Notifications
You must be signed in to change notification settings - Fork 0
/
hc_visc_scan.c
321 lines (276 loc) · 9.9 KB
/
hc_visc_scan.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
#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@ig.utexas.edu) for
additional comments, see hc.c
this binary scans through viscosities and computes correlation with
the geoid
output viscosities are log10(eta/1e21), top to bottom
and three correlations: r_20, r_4-9, r_2-4
*/
/* indices for arrays */
#define IVMIN 0
#define IVMAX 1
#define IDV 2
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;
struct hc_parameters p[1]; /* parameters */
hc_boolean solved = FALSE; /* init with FALSE! */
hc_boolean vary_umlm = FALSE;
HC_PREC dv_use,vl[HC_VSCAN_NLAYER_MAX][3],v[HC_VSCAN_NLAYER_MAX]; /* for viscosity scans */
strncpy(p->main_program_name,argv[0],HC_CHAR_LENGTH);
/*
(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; /* compute geoid at surface */
/*
handle command line arguments
*/
hc_handle_command_line(argc,argv,1,p);
fprintf(stderr,"%s: starting scan using reference %s, dv_ref: %g, nlayer: %i/%i, z_ulm: %g z_asth: %g\n",
argv[0],p->ref_geoid_file,(double)p->vscan_dv,p->vscan_n,
HC_VSCAN_NLAYER_MAX,
(double)HC_Z_DEPTH(p->rlayer[0]),
(double)HC_Z_DEPTH(p->rlayer[1]));
if(p->vscan_n < 0){
p->vscan_n = - p->vscan_n;
fprintf(stderr,"%s: i.e. %i layers, but additionally varying upper/lower mantle boundary\n",
argv[0],p->vscan_n);
vary_umlm = TRUE;
}
/*
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);
/* select plate velocity if not free slip */
if(!p->free_slip)
hc_select_pvel(p->pvel_time,&model->pvel,pvel,p->verbose);
/* parameter space log bounds */
dv_use = p->vscan_dv; /* don't modify */
switch(p->vscan_n){
case 4:
/*
for layer case
*/
/* uniform "priors" */
vl[0][IVMIN]= -p->vscan_em;vl[0][IVMAX]=p->vscan_em+1e-5;vl[0][IDV]=dv_use; /* 0..100
layer
log
bounds
and
spacing */
vl[1][IVMIN]= -p->vscan_em;vl[1][IVMAX]=p->vscan_em+1e-5;vl[1][IDV]=dv_use; /* 100..410 */
if(p->free_slip){
vl[2][IVMIN]= 0;vl[2][IVMAX]=0+1e-5;vl[2][IDV]=dv_use; /* for free
slip, only relative
viscosisites matter
for correlation */
fprintf(stderr,"%s: for free slip, we set upper mantle (layer 2) to unity (only relative viscosities matter)\n",argv[0]);
}else{
vl[2][IVMIN]= -p->vscan_em;vl[2][IVMAX]=p->vscan_em+1e-5;vl[2][IDV]=dv_use; /* need to actually
loop 410 .660 */
}
vl[3][IVMIN]= -p->vscan_em;vl[3][IVMAX]=p->vscan_em+1e-5;vl[3][IDV]=dv_use; /* 660 ... 2871 */
/* loop */
for(v[0]=vl[0][IVMIN];v[0] <= vl[0][IVMAX];v[0] += vl[0][IDV])
for(v[1]=vl[1][IVMIN];v[1] <= vl[1][IVMAX];v[1] += vl[1][IDV])
for(v[2]=vl[2][IVMIN];v[2] <= vl[2][IVMAX];v[2] += vl[2][IDV])
for(v[3]=vl[3][IVMIN];v[3] <= vl[3][IVMAX];v[3] += vl[3][IDV])
visc_scan_out(v,geoid,sol_spectral,pvel,p,model,&solved,vary_umlm);
break;
case 3:
dv_use /= 3; /* refine */
/*
three layer case
*/
vl[0][IVMIN]= -p->vscan_em;vl[0][IVMAX]=p->vscan_em+1e-5;vl[0][IDV]=dv_use; /* 0..100
layer
log
bounds
and
spacing */
/* vl[1] will be same as vl[2] */
if(p->free_slip){
vl[2][IVMIN]= 0;vl[2][IVMAX]=0+1e-5;vl[2][IDV]=dv_use; /* for free
slip, only relative
viscosisites matter
for correlation */
fprintf(stderr,"%s: for free slip, we set upper mantle (layer 2) to unity (only relative viscosities matter)\n",argv[0]);
}else{
vl[2][IVMIN]= -p->vscan_em;vl[2][IVMAX]=p->vscan_em+1e-5;vl[2][IDV]=dv_use; /* need to actually
loop 410 .660 */
}
vl[3][IVMIN]= -p->vscan_em;vl[3][IVMAX]=p->vscan_em+1e-5;vl[3][IDV]=dv_use; /* 660 ... 2871 */
/* loop */
for(v[0]=vl[0][IVMIN];v[0] <= vl[0][IVMAX];v[0] += vl[0][IDV])
for(v[2]=vl[2][IVMIN];v[2] <= vl[2][IVMAX];v[2] += vl[2][IDV]){
v[1] = v[2];
for(v[3]=vl[3][IVMIN];v[3] <= vl[3][IVMAX];v[3] += vl[3][IDV])
visc_scan_out(v,geoid,sol_spectral,pvel,p,model,&solved,vary_umlm);
}
break;
case 2:
/*
two layer case
vl[0] and vl[1] will be same as vl[2]
*/
dv_use /= 10; /* finer */
if(p->free_slip){
vl[2][IVMIN]= 0;vl[2][IVMAX]=0+1e-5;vl[2][IDV]=dv_use; /* for
free
slip,
only
relative
viscosisites
matter
for
correlation */
fprintf(stderr,"%s: for free slip, we set upper mantle (layer 2) to unity (only relative viscosities matter)\n",argv[0]);
}else{
vl[2][IVMIN]= -p->vscan_em;vl[2][IVMAX]=p->vscan_em+1e-5;vl[2][IDV]=dv_use; /* need to actually
loop 410 .660 */
}
vl[3][IVMIN]= -p->vscan_em;vl[3][IVMAX]=p->vscan_em+1e-5;vl[3][IDV]=dv_use; /* 660 ... 2871 */
/* loop */
for(v[2]=vl[2][IVMIN];v[2] <= vl[2][IVMAX];v[2] += vl[2][IDV]){
v[0] = v[1] = v[2];
for(v[3]=vl[3][IVMIN];v[3] <= vl[3][IVMAX];v[3] += vl[3][IDV])
visc_scan_out(v,geoid,sol_spectral,pvel,p,model,&solved,vary_umlm);
}
break;
default:
fprintf(stderr,"%s: not set up for %i layers\n",argv[0],p->vscan_n);
exit(-1);
break;
}
/*
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;
}
/*
print out a four layer viscosity structure geoid correlation suite,
or additionally scan through the upper/lower mantle depths
*/
void visc_scan_out (v, geoid, sol_spectral, pvel, p, model, solved, vary_umlm)
HC_PREC *v;
struct sh_lms *geoid;
struct sh_lms *sol_spectral;
struct sh_lms *pvel;
struct hc_parameters *p;
struct hcs *model;
hc_boolean *solved;
hc_boolean vary_umlm;
{
HC_PREC corr[3],r660=660,rms;
const HC_PREC rtop = 300.1, rbot = 1800+1e-5, dr = 25;
if(p->vscan_rlv){
if((v[0] < v[1])||(v[0] < v[2])) /* lithosphere should be > asth or upper mantle */
return;
if(v[1] > v[2]) /* asthenosphere should be < upper mantle */
return;
if(v[2] > v[3])
return;
}
if(vary_umlm){
for(r660=rtop;r660 <= rbot;r660 += dr){
/* overwrite 660 as first non CMB boundary from the bottom */
p->rlayer[0] = HC_ND_RADIUS(r660);
/* print viscosities of 0...100, 100...410, 410 ... 660 and
660...2871 layer in log space */
fprintf(stdout,"%14.7e %14.7e %14.7e %14.7e\t",
(double)v[0],(double)v[1],(double)v[2],(double)v[3]);
hc_calc_geoid_corr_four_layer(v,geoid,sol_spectral,pvel,p,model,solved,corr,&rms);
fprintf(stdout,"%10.7f %10.7f %10.7f\t%8.3f\t%.4e\n",
(double)corr[0],(double)corr[1],(double)corr[2],(double)r660,(double)rms);
}
}else{
/* no radius scan */
fprintf(stdout,"%14.7e %14.7e %14.7e %14.7e\t",
(double)v[0],(double)v[1],(double)v[2],(double)v[3]);
hc_calc_geoid_corr_four_layer(v,geoid,sol_spectral,pvel,p,model,solved,corr,&rms);
fprintf(stdout,"%10.7f %10.7f %10.7f\tNaN\t%.4e\n",
(double)corr[0],(double)corr[1],(double)corr[2],(double)rms);
}
}