New model for xrb_layered with a lot more carbon (#474)

The new model is at Exec/science/xrb_layered/create_initial_model/toy_atm/toy_atm_hot_carbon_3cm.hse

Some other updates to documentation in the readme's.

Exec/science/xrb_layered/README.md

@@ -4,6 +4,13 @@ Models the convection and burning in the envelope of a neutron star, where hydro
 has previously been depleted by the hot CNO cycle, leaving a pure helium layer at 
 large depths, where the burst will ignite.
 
+Appears in:
+
+  * `Hydrodynamical simulations of proton ingestion flashes in Type I X-ray Bursts`
+     Guichandut, S.; Zingale, M.; Cumming, A, 2023. 
+     https://arxiv.org/abs/2405.08952
+
+
 This is the multidimensional version of a calculation previously done in MESA in:
 
   * `The imprint of convection on Type I X-ray bursts: Pauses in photospheric
@@ -13,7 +20,6 @@ This is the multidimensional version of a calculation previously done in MESA in
 
 ## initial model
 
-* `xrb_layered_3cm.hse`
+* `toy_atm_hot_3cm.hse`
 
-An initial MESA model is placed into eulerian HSE using the lagrangian_planar routines
-from AMREX/initial_models. More detail in create_initial_model/README.md.
+Build using a MESA model as a guide. More details in `create_initial_model/`.

Exec/science/xrb_layered/create_initial_model/README.md

@@ -1,13 +1,11 @@
 ## MESA
 
-We first use MESA to create a model of a neutron star with mass $M=1.4M_\odot$ , radius $R=12$ km, and surface gravity $g=GM/R^2=1.29\times10^{14}$ g/cm $^2$ (MESA is Newtonian, and gravity has not been redshift-corrected). 
-
-We then accrete a column $y\sim 2\times 10^4$ g/cm $^{2}$ of iron. This will be the substrate/buffer below the burning region.
+We first use MESA to create a model of a neutron star with mass $M=1.4M_\odot$ , radius $R=12$ km, and surface gravity $g=GM/R^2=1.29\times10^{14}$ g/cm $^2$ (MESA is Newtonian, and gravity has not been redshift-corrected). Initially, the model is just an atmosphere of pure iron with a column depth of 5e9 g/cm2.
 
 We then accrete a "solar" mixture (70% hydrogen, 28% helium, 2% carbon-12) for a few days, until ignition in a helium layer. 
 
-The last two profiles of the MESA run (`profile37.data` & `profile38.data`) are copied in `mesa/`. These steps are all done in [github.com/simonguichandut/mesa\_mixed\_burst](https://github.com/simonguichandut/mesa_mixed_burst). The relevant run is `run/make_ignition_model_for_hydro/`. 
+The last two models of the MESA run (`profile4.data` & `profile5.data`) are copied in `mesa/`. These steps are all done in [github.com/simonguichandut/mesa\_mixed\_burst](https://github.com/simonguichandut/mesa_mixed_burst). The relevant run is `run/make_ignition_model_for_hydro/`. 
 
 ## Toy model
 
-See the notebook `Initial_Model.ipynb` for a detailed explanation of how we use the MESA model (profile37) to construct a simplified toy model, with a temperature kick at the bottom of the fuel layer. This is done in `toy_atm/`, which is a modified version of https://github.com/AMReX-Astro/initial_models/tree/main/toy_atm .
+See the notebook `Initial_Model.ipynb` for a detailed explanation of how we use the pre-igntion MESA model (profile4) to construct a simplified toy model, with a temperature kick at the bottom of the fuel layer. This is done in `toy_atm/`, which is a modified version of https://github.com/AMReX-Astro/initial_models/tree/main/toy_atm .

Exec/science/xrb_layered/create_initial_model/toy_atm/_parameters

@@ -33,6 +33,9 @@ X0             real         0.7
 Zcno           real         0.02
 yd             real         1.e7
 ash_name       character    "iron-56"
+add_carbon     integer      0
+carbon_y       real         5.e8
+carbon_del     real         3.0
 
 # Other
 index_base_from_temp  integer   0

Exec/science/xrb_layered/create_initial_model/toy_atm/init_1d.H

@@ -182,7 +182,7 @@ AMREX_INLINE void init_1d()
     // We know the pressure (from y1) and the temperature
     Real pres_base = -problem_rp::g_const * problem_rp::y1;
 
-    // However there is know (p,T) input option, so we need to iterate
+    // However there is no (p,T) input option, so we need to iterate
     // to find the base density.
     // We know that it is roughly 1.3-1.5e6 (Initial_Model notebook)
 
@@ -621,4 +621,4 @@ AMREX_INLINE void init_1d()
 
     std::cout << "maximum HSE error = " << max_hse_error << std::endl;
 
-}
+}