Merge pull request #98 from brittonsmith/newH2selfshielding

New h2selfshielding (reissue of PR #72)
grackle-project · Mar 9, 2022 · d0957c6 · d0957c6
2 parents a135b32 + e2fe9ca
commit d0957c6
Show file tree

Hide file tree

Showing 2 changed files with 30 additions and 4 deletions.
diff --git a/doc/source/Parameters.rst b/doc/source/Parameters.rst
@@ -369,7 +369,10 @@ For all on/off integer flags, 0 is off and 1 is on.
    density. Default: 0.
 
     - 1: Use a Sobolev-like, spherically averaged method from
-      `Wolcott-Green et. al. 2011 <http://adsabs.harvard.edu/abs/2011MNRAS.418..838W>`_.
+      `Wolcott-Green \& Haiman (2019)
+      <https://ui.adsabs.harvard.edu/abs/2019MNRAS.484.2467W/>`__. Prior to
+      Grackle version 3.2, this option used the method of `Wolcott-Green et. al.
+      (2011) <https://ui.adsabs.harvard.edu/abs/2011MNRAS.418..838W/>`__.
       This option is only valid for Cartesian grid codes in 3D.
     - 2: Supply an array of lengths using the :c:data:`H2_self_shielding_length`
       field.

diff --git a/src/clib/solve_rate_cool_g.F b/src/clib/solve_rate_cool_g.F
@@ -1176,6 +1176,10 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k,
      &       f_shield, b_doppler, l_H2shield
       real*8 k13_CID, k13_DT
 
+!     locals for H2 self-shielding as WG+19
+
+      real*8 tgas_touse, ngas_touse, aWG2019
+
       real*8 nSSh, nratio
 
 !     Set log values of start and end of lookup tables
@@ -1429,13 +1433,32 @@ subroutine lookup_cool_rates1d_g(temstart, temend, nratec, j, k,
 
                N_H2 = dom*H2I(i,j,k) * l_H2shield
 
+! update: self-shielding following Wolcott-Green & Haiman (2019)
+! range of validity: T=100-8000 K, n<=1e7 cm^-3
+
+               tgas_touse = max(tgas1d(i),1E2_DKIND)
+               tgas_touse = min(tgas_touse,8E3_DKIND)
+               ngas_touse = d(i,j,k) * dom / mmw(i)
+               ngas_touse = min(ngas_touse,1E7_DKIND)
+
+               aWG2019 = (0.8711_DKIND *
+     &              log10(tgas_touse) - 1.928_DKIND) *
+     &              exp(-0.2856_DKIND * log10(ngas_touse)) +
+     &              (-0.9639_DKIND * log10(tgas_touse) + 3.892_DKIND)
+
                x = 2.0E-15_DKIND * N_H2
                b_doppler = 1E-5_DKIND *
      &                 sqrt(2._DKIND * kboltz *
      &                      tgas1d(i) / mass_h)
-               f_shield = 0.965_DKIND / (1._DKIND + x/b_doppler)**1.1+
-     &           0.035_DKIND * exp(-8.5E-4_DKIND * sqrt(1._DKIND +x))/
-     &           sqrt(1._DKIND + x)
+               f_shield = 0.965_DKIND /
+     &              (1._DKIND + x/b_doppler)**aWG2019 +
+     &              0.035_DKIND * exp(-8.5e-4_DKIND *
+     &              sqrt(1._DKIND + x)) /
+     &              sqrt(1._DKIND + x)
+
+! avoid f>1
+               f_shield = min(f_shield, 1._DKIND)
+
                k31shield(i) = f_shield * k31shield(i)
             endif
          enddo