anchor.diff:    use initial_model_module, only: initial_model
   
     ! Probin file
  +
     character (len=:), allocatable, save :: probin
   
     ! For determining if we are the I/O processor.
  +  
     integer, save :: ioproc
     
  -  ! Initial binary star characteristics
  +  ! Initial stellar properties
     ! Note that the envelope mass is included within the total mass of the star
  -  double precision, save :: mass_P, mass_S
  -  double precision, save :: central_density_P, central_density_S
  -  double precision, save :: stellar_temp
  +  
  +  double precision, save :: mass_P = ONE
  +  double precision, save :: mass_S = ONE
  +  double precision, save :: central_density_P = -ONE
  +  double precision, save :: central_density_S = -ONE
  +  double precision, save :: stellar_temp = 1.0d7
     double precision, save :: primary_envelope_mass, secondary_envelope_mass
     double precision, save :: primary_envelope_comp(nspec), secondary_envelope_comp(nspec)
   
  +
  +  
     ! Ambient medium
  -  double precision, save :: ambient_density, ambient_temp, ambient_comp(nspec)
  -
  +  
  +  double precision, save :: ambient_density = 1.0d-4
  +  double precision, save :: ambient_temp = 1.0d7
  +  double precision, save :: ambient_comp(nspec)
  +
  +
  +  
     ! Smallest allowed mass fraction
  -  double precision, save :: smallx
  +  
  +  double precision, save :: smallx = 1.0d-12
   
     ! Smallest allowed velocity on the grid
  -  double precision, save :: smallu
  -  
  -  ! Controls interpolation from 1D model to 3D model
  -  integer, save :: nsub
  -  logical, save :: interp_temp
  -
  +  
  +  double precision, save :: smallu = 1.0d-12
  +
  +
  +  
  +  ! Parameters for nterpolation from 1D model to 3D model:
  +
  +  ! Number of sub-grid-scale zones to use
  +  
  +  integer, save :: nsub = 1
  +
  +  ! Default to interpolation that preserves temperature; otherwise, use pressure
  +  
  +  logical, save :: interp_temp = .true.
  +
  +
  +  
     ! Whether or not to give stars an initial orbital velocity
     ! consistent with their Keplerian orbit speed.
  -  logical, save :: no_orbital_kick
  -
  +  
  +  logical, save :: no_orbital_kick = .false.
  +
  +
  +  
  +  ! Collision parameters
  +  
     ! Whether or not to give the stars an initial velocity
     ! consistent with the free-fall speed.
  -  logical, save :: collision
  +  
  +  logical, save :: collision = .false.
   
     ! For a collision, number of (secondary) WD radii to 
     ! separate the WDs by.
  -  double precision, save :: collision_separation
  +  
  +  double precision, save :: collision_separation = 4.0d0
   
     ! For a collision, the impact parameter measured in
     ! units of the primary's initial radius.
  -  double precision, save :: collision_impact_parameter
  -  
  -  ! Binary properties
  +  
  +  double precision, save :: collision_impact_parameter = 0.0d0
  +
  +
  +  
  +  ! Binary orbit properties
  +  
     double precision, save :: r_P_initial, r_S_initial, a_P_initial, a_S_initial, a  
     double precision, save :: v_P_r, v_S_r, v_P_phi, v_S_phi
     double precision, save :: center_P_initial(3), center_S_initial(3)
  -  double precision, save :: orbital_eccentricity, orbital_angle
  -  
  -  integer, save :: axis_1, axis_2, axis_3
  +  double precision, save :: orbital_eccentricity = 0.0d0
  +  double precision, save :: orbital_angle = 0.0d0
  +    
  +
  +  
  +  ! Axis is in orbital plane; we measure angle with respect to this axis. Normally the x axis.
  +  
  +  integer, save :: axis_1 = 1
  +
  +  ! Perpendicular axis in the orbital plane. Normally the y axis.
  +  
  +  integer, save :: axis_2 = 2
  +
  +  ! Perpendicular to both other axies. Normally the z axis and also the rotation axis.
  +  
  +  integer, save :: axis_3 = 3
     
     ! Location of the physical center of the problem, as a fraction of domain size
  -  double precision, save :: center_fracx, center_fracy, center_fracz
  +  
  +  double precision, save :: center_fracx = HALF
  +  double precision, save :: center_fracy = HALF
  +  double precision, save :: center_fracz = HALF
   
     ! Bulk system motion
  -  double precision, save :: bulk_velx, bulk_vely, bulk_velz
  +  
  +  double precision, save :: bulk_velx = ZERO
  +  double precision, save :: bulk_vely = ZERO
  +  double precision, save :: bulk_velz = ZERO
   
     ! Whether we're doing an initialization or a restart
  +  
     integer, save :: init
   
     ! Are we doing a single star simulation?
  +  
     logical, save :: single_star
   
     ! Should we override the domain boundary conditions with
     ! ambient material?
  -  logical, save :: fill_ambient_bc
  -  
  +  
  +  logical, save :: fill_ambient_bc = .false.
  +
  +
  +
     ! 1D initial models
   
     type (initial_model) :: model_P, model_S
   
  -  double precision, save :: initial_model_dx
  -  integer, save          :: initial_model_npts
  -  double precision, save :: initial_model_mass_tol
  -  double precision, save :: initial_model_hse_tol
  -
  -  ! Composition properties of initial models
  -  double precision, save :: max_he_wd_mass
  -  double precision, save :: max_hybrid_co_wd_mass, hybrid_co_wd_he_shell_mass
  -  double precision, save :: max_co_wd_mass, co_wd_he_shell_mass
  -
  +  ! For the grid spacing for our model, we'll use 
  +  ! 6.25 km. No simulation we do is likely to have a resolution
  +  ! higher than that inside the stars (it represents
  +  ! three jumps by a factor of four compared to our 
  +  ! normal coarse grid resolution). By using 4096
  +  ! grid points, the size of the 1D domain will be 2.56e9 cm,
  +  ! which is larger than any reasonable mass white dwarf.
  +
  +  double precision, save :: initial_model_dx = 6.25d5
  +  integer, save          :: initial_model_npts = 4096
  +
  +  ! initial_model_mass_tol is tolerance used for getting the total WD mass 
  +  ! equal to the desired mass. It can be reasonably small, since there
  +  ! will always be a central density value that can give the desired
  +  ! WD mass on the grid we use.
  +  
  +  double precision, save :: initial_model_mass_tol = 1.d-6
  +
  +  ! hse_tol is the tolerance used when iterating over a zone to force
  +  ! it into HSE by adjusting the current density (and possibly
  +  ! temperature).  hse_tol should be very small (~ 1.e-10).
  +  
  +  double precision, save :: initial_model_hse_tol = 1.d-10
  +
  +
  +  
  +  ! Composition properties of initial models.
  +  ! We follow the prescription of Dan et al. 2012 for determining
  +  ! the composition of the WDs. In this approach, below a certain 
  +  ! mass there are pure He WDs; just above that are hybrid WDs
  +  ! with pure CO cores and a He mantle; above that are pure CO WDs
  +  ! with slightly more oxygen than carbon; and above that are 
  +  ! ONeMg WDs. All masses are in solar masses.
  +
  +  double precision, save :: max_he_wd_mass = 0.45d0
  +  double precision, save :: max_hybrid_co_wd_mass = 0.6d0
  +  double precision, save :: hybrid_co_wd_he_shell_mass = 0.1d0
  +  double precision, save :: max_co_wd_mass = 1.05d0
  +  double precision, save :: co_wd_he_shell_mass = 0.0d0
  +
  +
  +  
     ! Tagging criteria
  -  double precision, save :: max_tagging_radius
  -  double precision, save :: stellar_density_threshold
  -
  +  
  +  double precision, save :: max_tagging_radius = 0.75d0
  +  double precision, save :: stellar_density_threshold = 1.0d0
  +
  +
  +  
     ! Stores the center of mass location of the stars throughout the run
  +  
     double precision, save :: com_P(3), com_S(3)
     double precision, save :: vel_P(3), vel_S(3)
   
     ! Stores the effective Roche radii
  +  
     double precision, save :: roche_rad_P, roche_rad_S
   
     ! Lagrange point zone indices
  +  
     integer, save          :: L1_idx(3)
   
     ! Relaxation parameters
  -  logical, save          :: do_initial_relaxation
  -  double precision, save :: relaxation_timescale
  -  double precision, save :: relaxation_density_cutoff
  +  
  +  logical, save          :: do_initial_relaxation = .false.
  +  double precision, save :: relaxation_timescale = 0.001
  +  double precision, save :: relaxation_density_cutoff = 1.0d0
   
     ! Input parameters for SCF relaxation
  -  logical, save          :: do_scf_initial_models
  -  double precision, save :: scf_d_A, scf_d_B, scf_d_C
  -  double precision, save :: scf_relax_tol
  +  
  +  logical, save          :: do_scf_initial_models = .false.
  +  double precision, save :: scf_d_A = 1.0d9
  +  double precision, save :: scf_d_B = 1.0d9
  +  double precision, save :: scf_d_C = 1.8d9
  +  double precision, save :: scf_relax_tol = 1.d-3
     
     ! Internal data for SCF relaxation
  +  
     double precision, save :: scf_h_max_P, scf_h_max_S
     double precision, save :: scf_enthalpy_min
     double precision, save :: scf_rpos(3,3)         ! Relative position of points A, B, and C
  
positive.code:          gw_dist, &
           fill_ambient_bc
  
      nsub = 1
  
      central_density_P = -ONE
      central_density_S = -ONE
  
      mass_P = 1.0d0
      mass_S = 1.0d0
  
      stellar_temp = 1.0d7
      ambient_temp = 1.0d7
  
      ! We follow the prescription of Dan et al. 2012 for determining
      ! the composition of the WDs. In this approach, below a certain 
      ! mass there are pure He WDs; just above that are hybrid WDs
      ! with pure CO cores and a He mantle; above that are pure CO WDs
      ! with slightly more oxygen than carbon; and above that are 
      ! ONeMg WDs. All masses are in solar masses.
  
      max_he_wd_mass = 0.45d0
      max_hybrid_co_wd_mass = 0.6d0
      hybrid_co_wd_he_shell_mass = 0.1d0
      max_co_wd_mass = 1.05d0
      co_wd_he_shell_mass = 0.0d0
  
      smallx = 1.d-12
      smallu = 1.d-12
  
      ambient_density = 1.d-4
  
      no_orbital_kick = .false.
  
      collision = .false.
      collision_separation = 4.0
      collision_impact_parameter = 0.0
  
      interp_temp = .false.
  
      orbital_eccentricity = 0.0d0
      orbital_angle = 0.0d0
      
      axis_1 = 1 ! Axis is in orbital plane; we measure angle with respect to this axis. Normally the x axis.
      axis_2 = 2 ! Perpendicular axis in the orbital plane. Normally the y axis.
      axis_3 = 3 ! Perpendicular to both other axies. Normally the z axis and also the rotation axis.
      
      max_tagging_radius = 0.75d0
      stellar_density_threshold = 1.0d0
  
      do_initial_relaxation = .false.
      relaxation_timescale = 0.001
      relaxation_density_cutoff = 1.0d0
  
      do_scf_initial_models = .false.
      scf_d_A = 1.0d9
      scf_d_B = 1.0d9
      scf_d_C = 1.8d9
      scf_relax_tol = 1.d-3
  
      bulk_velx = ZERO
      bulk_vely = ZERO
      bulk_velz = ZERO
  
      center_fracx = HALF
      center_fracy = HALF
      center_fracz = HALF
  
      ! For the grid spacing for our model, we'll use 
      ! 6.25 km. No simulation we do is likely to have a resolution
      ! higher than that inside the stars (it represents
      ! three jumps by a factor of four compared to our 
      ! normal coarse grid resolution). By using 2048
      ! grid points, the size of the 1D domain will be 2.56e9 cm,
      ! which is larger than any reasonable mass white dwarf.
  
      initial_model_dx = 6.25d5
      initial_model_npts = 4096
  
      ! initial_model_mass_tol is tolerance used for getting the total WD mass 
      ! equal to the desired mass. It can be reasonably small, since there
      ! will always be a central density value that can give the desired
      ! WD mass on the grid we use.
  
      initial_model_mass_tol = 1.d-6
  
      ! hse_tol is the tolerance used when iterating over a zone to force
      ! it into HSE by adjusting the current density (and possibly
      ! temperature).  hse_tol should be very small (~ 1.e-10).
  
      initial_model_hse_tol = 1.d-10
  
      gw_dist = 10.0 ! kpc
  
      fill_ambient_bc = .false.
      
      ! Read namelist to override the defaults.
  
      untin = 9 
      open(untin,file=probin,form='formatted',status='old')
  
negative.code: 
    ! Distance (in kpc) used for calculation of the gravitational wave amplitude
    ! (this wil be calculated along all three coordinate axes).
    
    double precision, save :: gw_dist = 10.0d0
  
  contains
  
  