Difficulty in reaction rate calculation

[ktravis213]

Your name

Katherine Travis

Your affiliation

NASA LaRC

What happened? What did you expect to happen?

I wrote a reaction rate in fullchem_RateLaws.F90 as follows:

       k4 = 7.45e+7_dp *exp(-4430_dp*(1/TEMP-1/298_dp))

I also tried to separate it into parts.

       k4b = (TEMP-1/298_dp)
       k4c = exp(-4430_dp*k4b)
       k4 = 7.45e+7_dp *k4c 

In either case, if T is 298.55, the answer should be k4 = 7.65E7. However the code outputs k4 = 16.89. Printing out k4b and k4c above, the issue is in k4b which is calculated by fortran as 3.45E-3 when it should be 9.79E-5. Is there something wrong with my organization in fortran that causes this?

289.55118 --

What are the steps to reproduce the bug?

Print out the following somewhere in fullchem_RateLaws.F90

       k4b = (TEMP-1/298_dp)
       k4c = exp(-4430_dp*k4b)
       k4 = 7.45e+7_dp *k4c 

Calculate in excel or somewhere else the same code. Hopefully that reproduces my issue.

Please attach any relevant configuration and log files.

No response

What GEOS-Chem version were you using?

14.4.1

What environment were you running GEOS-Chem on?

Local cluster

What compiler and version were you using?

ifort: intel/19.1.3.304

Will you be addressing this bug yourself?

Yes, but I will need some help

In what configuration were you running GEOS-Chem?

GCClassic

What simulation were you running?

Full chemistry

As what resolution were you running GEOS-Chem?

0.25 Asia

What meterology fields did you use?

GEOS-FP

Additional information

No response

[ktravis213]

Noticed that I had a typo when I made this issue - k4b = (1/TEMP-1/298_dp).

[yantosca]

Thanks for writing @ktravis213. I've made a test program with your equations as you have them:

PROGRAM Katie

  IMPLICIT NONE

  INTEGER, PARAMETER :: dp = KIND( 0.0d0 )   ! Double precision

  ! Test with T = 298.55
  CALL PrintValues( 298.55_dp )

  CONTAINS

    SUBROUTINE PrintValues( temp )

      ! Argument
      REAL(dp), INTENT(IN) :: temp

      ! Local variables
      REAL(dp) :: k4, k4b, k4c

      ! Temparature
      PRINT*, REPEAT( '#', 60 )
      PRINT*, '### temp               : ', temp
      PRINT*, '###'

      ! 1-line calculation
      k4 = 7.45e+7_dp * exp(-4430_dp*(1/TEMP-1/298_dp))
      PRINT*, '### 1/TEMP                 : ', 1/TEMP
      print*, '### 1/298_dp               : ', 1/298_dp
      PRINT*, '### 1/TEMP-1/298_dp        : ', 1/TEMP - 1/298_dp
      PRINT*, '### exponential term       : ', exp(-4430_dp*(1/TEMP - 1/298_dp))
      PRINT*, '### k4, in one line        : ', k4
      PRINT*, '###'

      ! 1-line calculation, with e.g. 298.0_dp
      k4 = 7.45e+7_dp * exp(-4430_dp*(1.0_dp/TEMP - 1.0_dp/298.0_dp))
      PRINT*, '### 1.0_dp/TEMP             : ', 1.0_dp/TEMP            
      PRINT*, '### 1.0_dp/298_dp           : ', 1.0_dp/298.0_dp
      PRINT*, '### 1.0_dp/TEMP-1.0_/298_dp : ', 1.0_dp/TEMP - 1.0_dp/298.0_dp
      PRINT*, '### exponential term        : ', exp(-4430_dp*(1.0_dp/TEMP - 1.0_dp/298.0_dp))
      PRINT*, '### k4, in one line         : ', k4

    END SUBROUTINE PrintValues

END PROGRAM Katie

Which gives as output:

 ############################################################
 ### temp               :    298.55000000000001     
 ###
 ### 1/TEMP                 :    3.3495226930162453E-003
 ### 1/298_dp               :                     0
 ### 1/TEMP-1/298_dp        :    3.3495226930162453E-003
 ### exponential term       :    3.5955972859675217E-007
 ### k4, in one line        :    26.787199780458035     
 ###
 ### 1.0_dp/TEMP             :    3.3495226930162453E-003
 ### 1.0_dp/298_dp           :    3.3557046979865771E-003
 ### 1.0_dp/TEMP-1.0_/298_dp :   -6.1820049703318475E-006
 ### exponential term        :    1.0277647331307509     
 ### k4, in one line         :    76568472.618240938     

So what is happening is that this is integer division. I think you need to use e.g. 298.0_dp instead of 298_dp (since the compiler will think that is an integer because it doesn't have a decimal place). Otherwise this will do an integer division, where 1/298 = 0.

[ktravis213]

Thank you so much! I can't believe I couldn't figure this out!!!

[yantosca]

No worries @ktravis213. I guessed that it had something to do with operator precedence but I wouldn't have found it unless I made the toy example program. I'll close this out.