/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
*
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
* Clint Conrad, Michael Gurnis, and Eun-seo Choi.
* Copyright (C) 1994-2005, California Institute of Technology.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include
#include
#include "global_defs.h"
#include "lith_age.h"
#include "parsing.h"
void parallel_process_termination();
void temperatures_conform_bcs();
#include "initial_temperature.h"
void debug_tic(struct All_variables *);
void read_tic_from_file(struct All_variables *);
#ifdef USE_GZDIR
void restart_tic_from_gzdir_file(struct All_variables *);
#endif
void tic_input(struct All_variables *E)
{
int m = E->parallel.me;
int noz = E->lmesh.noz;
int n;
input_int("tic_method", &(E->convection.tic_method), "0,0,2", m);
/* When tic_method is 0 (default), the temperature is a linear profile +
perturbation at some layers.
When tic_method is -1, the temperature is read in from the
[datafile_old].velo.[rank].[solution_cycles_init] files.
When tic_method is 1, the temperature is isothermal (== bottom b.c.) +
uniformly cold plate (thickness specified by 'half_space_age').
When tic_method is 2, (tic_method==1) + a hot blob. A user can specify
the location and radius of the blob, and also the amplitude of temperature
change in the blob relative to the ambient mantle temperautre
(E->control.lith_age_mantle_temp).
- blob_center: A comma-separated list of three float numbers.
- blob_radius: A dmensionless length, typically a fraction
of the Earth's radius.
- blob_dT : Dimensionless temperature.
When tic_method is 3, the temperature is a linear profile + perturbation
for whole mantle.
tic_method is 4: read in initial temperature distribution from a set of netcdf grd
files. this required the GGRD extension to be compiled in
*/
switch(E->convection.tic_method){
case -1: /* read from file, no other options needed */
break;
case 0:
case 3:
/* This part put a temperature anomaly at depth where the global
node number is equal to load_depth. The horizontal pattern of
the anomaly is given by spherical harmonic ll & mm. */
input_int("num_perturbations", &n, "0,0,PERTURB_MAX_LAYERS", m);
if (n > 0) {
E->convection.number_of_perturbations = n;
if (! input_float_vector("perturbmag", n, E->convection.perturb_mag, m) ) {
fprintf(stderr,"Missing input parameter: 'perturbmag'\n");
parallel_process_termination();
}
if (! input_int_vector("perturbm", n, E->convection.perturb_mm, m) ) {
fprintf(stderr,"Missing input parameter: 'perturbm'\n");
parallel_process_termination();
}
if (! input_int_vector("perturbl", n, E->convection.perturb_ll, m) ) {
fprintf(stderr,"Missing input parameter: 'perturbl'\n");
parallel_process_termination();
}
if (! input_int_vector("perturblayer", n, E->convection.load_depth, m) ) {
fprintf(stderr,"Missing input parameter: 'perturblayer'\n");
parallel_process_termination();
}
}
else {
E->convection.number_of_perturbations = 1;
E->convection.perturb_mag[0] = 1;
E->convection.perturb_mm[0] = 2;
E->convection.perturb_ll[0] = 2;
E->convection.load_depth[0] = (noz+1)/2;
}
break;
case 1: /* case 1 */
input_float("half_space_age", &(E->convection.half_space_age), "40.0,1e-3,nomax", m);
break;
case 2: /* case 2 */
input_float("half_space_age", &(E->convection.half_space_age), "40.0,1e-3,nomax", m);
if( ! input_float_vector("blob_center", 3, E->convection.blob_center, m)) {
assert( E->sphere.caps == 12 || E->sphere.caps == 1 );
if(E->sphere.caps == 12) { /* Full version: just quit here */
fprintf(stderr,"Missing input parameter: 'blob_center'.\n");
parallel_process_termination();
}
else if(E->sphere.caps == 1) { /* Regional version: put the blob at the center */
fprintf(stderr,"Missing input parameter: 'blob_center'. The blob will be placed at the center of the domain.\n");
E->convection.blob_center[0] = 0.5*(E->control.theta_min+E->control.theta_max);
E->convection.blob_center[1] = 0.5*(E->control.fi_min+E->control.fi_max);
E->convection.blob_center[2] = 0.5*(E->sphere.ri+E->sphere.ro);
}
}
input_float("blob_radius", &(E->convection.blob_radius), "0.063,0.0,1.0", m);
input_float("blob_dT", &(E->convection.blob_dT), "0.18,nomin,nomax", m);
break;
case 4:
/*
case 4: initial temp from grd files
*/
#ifdef USE_GGRD
/* read in some more parameters */
/* scale the anomalies with PREM densities */
input_boolean("ggrd_tinit_scale_with_prem",&(E->convection.ggrd_tinit_scale_with_prem),"off",E->parallel.me);
/* limit T to 0...1 */
input_boolean("ggrd_tinit_limit_trange",&(E->convection.ggrd_tinit_limit_trange),"on",E->parallel.me);
/* scaling factor for the grids */
input_double("ggrd_tinit_scale",&(E->convection.ggrd_tinit_scale),"1.0",E->parallel.me); /* scale */
/* temperature offset factor */
input_double("ggrd_tinit_offset",&(E->convection.ggrd_tinit_offset),"0.0",E->parallel.me); /* offset */
/* grid name, without the .i.grd suffix */
input_string("ggrd_tinit_gfile",E->convection.ggrd_tinit_gfile,"",E->parallel.me); /* grids */
input_string("ggrd_tinit_dfile",E->convection.ggrd_tinit_dfile,"",E->parallel.me); /* depth.dat layers of grids*/
/* override temperature boundary condition? */
input_boolean("ggrd_tinit_override_tbc",&(E->convection.ggrd_tinit_override_tbc),"off",E->parallel.me);
input_string("ggrd_tinit_prem_file",E->convection.prem.model_filename,"", E->parallel.me); /* PREM model filename */
#else
fprintf(stderr,"tic_method 4 only works for USE_GGRD compiled code\n");
parallel_process_termination();
#endif
break;
default: /* unknown option */
fprintf(stderr,"Invalid value of 'tic_method'\n");
parallel_process_termination();
break;
}
return;
}
void convection_initial_temperature(struct All_variables *E)
{
void report();
report(E,"Initialize temperature field");
if (E->convection.tic_method == -1) {
/* read temperature from file */
#ifdef USE_GZDIR
if(strcmp(E->output.format, "ascii-gz") == 0)
restart_tic_from_gzdir_file(E);
else
#endif
read_tic_from_file(E);
}
else if (E->control.lith_age)
lith_age_construct_tic(E);
else
(E->solver.construct_tic_from_input)(E);
/* Note: it is the callee's responsibility to conform tbc. */
/* like a call to temperatures_conform_bcs(E); */
if (E->control.verbose)
debug_tic(E);
return;
}
void debug_tic(struct All_variables *E)
{
int m, j;
fprintf(E->fp_out,"output_temperature\n");
for(m=1;m<=E->sphere.caps_per_proc;m++) {
fprintf(E->fp_out,"for cap %d\n",E->sphere.capid[m]);
for (j=1;j<=E->lmesh.nno;j++)
fprintf(E->fp_out,"X = %.6e Z = %.6e Y = %.6e T[%06d] = %.6e \n",E->sx[m][1][j],E->sx[m][2][j],E->sx[m][3][j],j,E->T[m][j]);
}
fflush(E->fp_out);
return;
}
void read_tic_from_file(struct All_variables *E)
{
void temperatures_conform_bcs();
int ii, ll, mm;
float tt;
int i, m;
char output_file[255], input_s[1000];
FILE *fp;
float v1, v2, v3, g;
ii = E->monitor.solution_cycles_init;
sprintf(output_file,"%s.velo.%d.%d",E->control.old_P_file,E->parallel.me,ii);
fp=fopen(output_file,"r");
if (fp == NULL) {
fprintf(E->fp,"(Initial_temperature.c #1) Cannot open %s\n",output_file);
parallel_process_termination();
}
if (E->parallel.me==0)
fprintf(E->fp,"Reading %s for initial temperature\n",output_file);
fgets(input_s,1000,fp);
sscanf(input_s,"%d %d %f",&ll,&mm,&tt);
for(m=1;m<=E->sphere.caps_per_proc;m++) {
fgets(input_s,1000,fp);
sscanf(input_s,"%d %d",&ll,&mm);
for(i=1;i<=E->lmesh.nno;i++) {
fgets(input_s,1000,fp);
sscanf(input_s,"%g %g %g %f",&(v1),&(v2),&(v3),&(g));
/* Truncate the temperature to be within (0,1). */
/* This might not be desirable in some situations. */
E->T[m][i] = max(0.0,min(g,1.0));
}
}
fclose (fp);
temperatures_conform_bcs(E);
return;
}