Skip to content

Commit

Permalink
Add bug fix for NEW_HENRY_CONSTANTS
Browse files Browse the repository at this point in the history
Jenny Fisher wrote:

   I am a bit confused about the code used for NEW_HENRY_CONSTANTS, and I
   think there may be a bug here (not 100% sure).

   If you take acetone as an example, we have:
     #if defined( NEW_HENRY_CONSTANTS )
             Henry_K0      = 2.7e-1_f8 * To_M_atm
     #else
             Henry_K0      = 2.7e+1_f8

   Earlier in the code, we define the parameter:
     To_M_atm = 9.86923e-3_f8  ! mol/m3/Pa -> M/atm

   So the multiplication for NEW_HENRY_CONSTANTS of 2.7e-1_f8 * To_M_atm
   would make new Henry_K0 = 2.7e-3 not 2.7e+1.

   Should this code actually be:
     #if defined( NEW_HENRY_CONSTANTS )
             Henry_K0      = 2.7e-1_f8 / To_M_atm
     #else
             Henry_K0      = 2.7e+1_f8

   I think the latter also makes sense going back to the original unit
   analysis:

     X mol/m3/Pa * (1 m3/1e6 cm3) * (1 cm/1 mL) * (1e3 mL / 1 L) * (1.01325e5 Pa/1 atm)
      = X mol/m3/Pa * 1.01325e2
      = X mol/m3/Pa * 101.325
      = X mol/m3/Pa / 9.86923e-3
      = X mol/m3/Pa / To_M_atm

   I checked the current v12.4.0 code on GitHub and it doesn’t look like
   this has been changed... I know this flag currently isn’t used but it
   looks a little dangerous if it is wrong.

Jenny Fisher followed up:

   I double-checked with Katie Travis and she said she just provided the
   table of constants, not the actual code. She also confirmed that it is
   a bug and should be a division rather than multiplication.

   I understand that this code is going to be replaced, but it seems to me
   like in the interim a bug fix should go in for anyone who IS currently
   using the NEW_HENRY_CONSTANTS flag (since it is available).

   A very easy way to do this without having to wade through the code would
   be to replace the current parameter with its inverse, e.g. change:
     To_M_atm = 9.86923e-3_f8  ! mol/m3/Pa -> M/atm

   to
     To_M_atm = 101.325_f8  ! mol/m3/Pa -> M/atm

Signed-off-by: Melissa Sulprizio <mpayer@seas.harvard.edu>
  • Loading branch information
msulprizio committed Jul 19, 2019
1 parent e6911c7 commit a905a3d
Showing 1 changed file with 1 addition and 1 deletion.
2 changes: 1 addition & 1 deletion Headers/species_database_mod.F90
Expand Up @@ -561,7 +561,7 @@ SUBROUTINE Spc_Info( am_I_Root, iName, Input_Opt, &
! Local parameter
LOGICAL, PARAMETER :: T = .TRUE. ! Yes
LOGICAL, PARAMETER :: F = .FALSE. ! No
REAL(f8), PARAMETER :: To_M_atm = 9.86923e-3_f8 ! mol/m3/Pa -> M/atm
REAL(f8), PARAMETER :: To_M_atm = 101.325_f8 ! mol/m3/Pa -> M/atm

!=====================================================================
! Spc_Info begins here!
Expand Down

0 comments on commit a905a3d

Please sign in to comment.