Changeset 2833 for trunk/LMDZ.MARS/libf

Timestamp:

Nov 28, 2022, 1:02:29 PM (2 years ago)

Author:

flefevre

Message:

Kinetics data: JPL 2019 update

File:

• trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90 (modified) (9 diffs)

trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90

@@ -380 +380 @@
 
 hetero_dust = .false.
-hetero_ice  = .false.
+hetero_ice  = .true.
 
 call reactionrates(nlayer, ionchem, deutchem,                             &
@@ -706 +706 @@
 real :: ak0, ak1, xpo, rate, rate1, rate2
 real :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y, gam
+real :: k0, kinf, kf, kint, kca
 real, dimension(nlayer) :: deq
 real, dimension(nlayer) :: a001, a002, a003,                           &
@@ -752 +753 @@
 !     from sehested et al., j. geophys. res., 100, 1995.
 
-      a001(:) = 2.075*6.0e-34*(t(:)/300.)**(-2.4)*dens(:)
+!     a001(:) = 2.075*6.0e-34*(t(:)/300.)**(-2.4)*dens(:)
+
+!     jpl 2019
+
+      a001(:) = 2.075*6.1e-34*(t(:)/298.)**(-2.4)*dens(:)
 
       nb_reaction_4 = nb_reaction_4 + 1
@@ -1003 +1008 @@
       do ilev = 1,lswitch-1
 !        ak0 = 3.1*2.4*4.4e-32*(t(ilev)/300.)**(-1.3) ! FL li et al 2017
-         ak0 = 2.4*4.4e-32*(t(ilev)/300.)**(-1.3)
-         ak1 = 7.5e-11*(t(ilev)/300.)**(0.2)
+!        ak0 = 2.4*4.4e-32*(t(ilev)/300.)**(-1.3)
+!        ak1 = 7.5e-11*(t(ilev)/300.)**(0.2)
+
+!     jpl 2019
+
+         ak0 = 2.4*5.3e-32*(t(ilev)/298.)**(-1.8)
+         ak1 = 9.5e-11*(t(ilev)/298.)**(0.4)
 
          rate = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
@@ -1093 +1103 @@
 !     jpl 2006
 
-      d001(:) = 5.1e-12*exp(210./t(:))
-
+!     d001(:) = 5.1e-12*exp(210./t(:))
+
+!     jpl 2019
+
+!     For the sake of simplicity, it is assumed that the association yield (kf) 
+!     gives the same product as the chemical activation yield (kca). 
+!     Thus the only products are no + o2. There is no production of no3.
+
+      do ilev = 1,lswitch-1
+
+!        association
+
+         k0 = 2.5*3.4e-31*(298./t(ilev))**(1.6)
+         kinf = 2.3e-11*(298./t(ilev))**(0.2)
+
+         kf = (kinf*k0*dens(ilev)/(kinf + k0*dens(ilev)))                           &
+                *0.6**(1. + (log10(k0*dens(ilev)/kinf))**2.)**(-1.0)
+
+!        chemical activation
+
+         kint = 5.3e-12*exp(200./t(ilev))
+
+         kca = kint*(1. - kf/kinf)
+
+!        total : association + chemical activation
+
+         d001(ilev) = kf + kca
+
+      end do
+    
       nb_reaction_4 = nb_reaction_4 + 1
       v_4(:,nb_reaction_4) = d001(:)
@@ -1111 +1149 @@
 !     jpl 2011
 
-      d003(:) = 3.3e-12*exp(270./t(:))
+!     d003(:) = 3.3e-12*exp(270./t(:))
+
+!     jpl 2019
+
+      d003(:) = 3.44e-12*exp(260./t(:))
 
       nb_reaction_4 = nb_reaction_4 + 1
@@ -1129 +1171 @@
 !     jpl 2011
 
-      d005(:) = 1.5e-11*exp(-3600./t(:))
+!     d005(:) = 1.5e-11*exp(-3600./t(:))
+
+!     jpl 2019
+
+      d005(:) = 3.3e-12*exp(-3150./t(:))
 
       nb_reaction_4 = nb_reaction_4 + 1
@@ -1138 +1184 @@
 !     jpl 2011
 
-      d006(:) = 4.0e-10*exp(-340./t(:))
+!     d006(:) = 4.0e-10*exp(-340./t(:))
+
+!     jpl 2019
+
+      d006(:) = 1.35e-10
 
       nb_reaction_4 = nb_reaction_4 + 1
@@ -1202 +1252 @@
 !---  e001: oh + co -> co2 + h
 
-!     jpl 2003
-
-!     e001(:) = 1.5e-13*(1 + 0.6*press(:)/1013.)
-
-!     mccabe et al., grl, 28, 3135, 2001
-
-!     e001(:) = 1.57e-13 + 3.54e-33*dens(:)
-
 !     jpl 2015
+!
+!     do ilev = 1,lswitch-1
+!
+!        branch 1 : oh + co -> h + co2
+!
+!        rate1 = 1.5e-13*(t(ilev)/300.)**(0.0)
+!
+!        branch 2 : oh + co + m -> hoco + m
+!
+!        ak0 = 5.9e-33*(t(ilev)/300.)**(-1.0)
+!        ak1 = 1.1e-12*(t(ilev)/300.)**(1.3)
+!        rate2 = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
+!        xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
+!
+!        e001(ilev) = rate1 + rate2*0.6**xpo
+!     end do
+
+!     jpl 2019
+
+!     For the sake of simplicity, it is assumed that the association yield (kf)
+!     gives the same product as the chemical activation yield (kca).
+!     Thus the only products are h + co2. There is no production of hoco.
 
       do ilev = 1,lswitch-1
 
-!        branch 1 : oh + co -> h + co2
-
-         rate1 = 1.5e-13*(t(ilev)/300.)**(0.0)
-
-!        branch 2 : oh + co + m -> hoco + m
-
-         ak0 = 5.9e-33*(t(ilev)/300.)**(-1.0)
-         ak1 = 1.1e-12*(t(ilev)/300.)**(1.3)
-         rate2 = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
-         xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
-
-         e001(ilev) = rate1 + rate2*0.6**xpo
+!        association
+
+         k0 = 2.5*6.9e-33*(298./t(ilev))**(2.1)
+         kinf = 1.1e-12*(298./t(ilev))**(-1.3)
+
+         kf = (kinf*k0*dens(ilev)/(kinf + k0*dens(ilev)))                           &
+                *0.6**(1. + (log10(k0*dens(ilev)/kinf))**2.)**(-1.0)
+
+!        chemical activation
+
+         kint = 1.85e-13*exp(-65./t(ilev))
+
+         kca = kint*(1. - kf/kinf)
+
+!        total : association + chemical activation
+
+         e001(ilev) = kf + kca
+
       end do
-
-!     joshi et al., 2006
-
-!     do ilev = 1,lswitch-1
-!        k1a0 = 1.34*2.5*dens(ilev)                                  &
-!              *1/(1/(3.62e-26*t(ilev)**(-2.739)*exp(-20./t(ilev)))  &
-!              + 1/(6.48e-33*t(ilev)**(0.14)*exp(-57./t(ilev))))     ! typo in paper corrected
-!        k1b0 = 1.17e-19*t(ilev)**(2.053)*exp(139./t(ilev))          &
-!             + 9.56e-12*t(ilev)**(-0.664)*exp(-167./t(ilev))
-!        k1ainf = 1.52e-17*t(ilev)**(1.858)*exp(28.8/t(ilev))        &
-!               + 4.78e-8*t(ilev)**(-1.851)*exp(-318./t(ilev))
-!        x = k1a0/(k1ainf - k1b0)
-!        y = k1b0/(k1ainf - k1b0)
-!        fc = 0.628*exp(-1223./t(ilev)) + (1. - 0.628)*exp(-39./t(ilev))  &
-!           + exp(-t(ilev)/255.)
-!        fx = fc**(1./(1. + (alog(x))**2))                           ! typo in paper corrected
-!        k1a = k1a0*((1. + y)/(1. + x))*fx
-!        k1b = k1b0*(1./(1.+x))*fx
-!        e001(ilev) = k1a + k1b
-!     end do
 
       nb_reaction_4 = nb_reaction_4 + 1