anchor.diff:           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.
  +    ! Read namelist to override the module defaults.
   
       untin = 9 
       open(untin,file=probin,form='formatted',status='old')
  
positive.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
  
  contains
  
  
negative.code:   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 stellar properties
    ! Note that the envelope mass is included within the total mass of the star
    
    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 = 1.0d-4
    double precision, save :: ambient_temp = 1.0d7
    double precision, save :: ambient_comp(nspec)
  
  
    
    ! Smallest allowed mass fraction
    
    double precision, save :: smallx = 1.0d-12
  
    ! Smallest allowed velocity on the grid
    
    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 = .false.
  
  
    
    ! Collision parameters
    
    ! Whether or not to give the stars an initial velocity
    ! consistent with the free-fall speed.
    
    logical, save :: collision = .false.
  
    ! For a collision, number of (secondary) WD radii to 
    ! separate the WDs by.
    
    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 = 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 = 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 = HALF
    double precision, save :: center_fracy = HALF
    double precision, save :: center_fracz = HALF
  
    ! Bulk system motion
    
    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 = .false.
  
  
  
    ! 1D initial models
  
    type (initial_model) :: model_P, model_S
  
    ! 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 = 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 = .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 = .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
  