Changeset 2433


Ignore:
Timestamp:
Nov 14, 2020, 4:51:23 PM (4 years ago)
Author:
flefevre
Message:

1) ajout de la photodissociation de HDO (commentee pour l'instant)

  • necessite la presence des sections efficaces de HDO dans le datadir:


datadir/hdo_composite_295K.txt

2) CO + OH -> CO2 + H : cinetique JPL 2015 (attention ce changement affecte les resultats)
3) chimie heterogene : inactive par defaut.

Location:
trunk/LMDZ.MARS/libf/aeronomars
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90

    r2321 r2433  
    718718
    719719            call photochemistry(nlayer,nq,nbq,ionchem,nb_reaction_3_max,  &
    720                            nb_reaction_4_max, nb_phot_max, nphotion,      &
    721                            jonline, ig,lswitch,zycol,szacol,ptimestep,    &
     720                           nb_reaction_4_max,nphot,nb_phot_max,nphotion,  &
     721                           jonline,ig,lswitch,zycol,szacol,ptimestep,     &
    722722                           zpress,zlocal,ztemp,ztelec,zdens,zmmean,       &
    723723                           dist_sol,zday,surfdust1d,surfice1d,            &

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

    r2321 r2433  
    33!     Photochemical routine
    44!
    5 !     Author: Franck Lefevre
    6 !     ------
     5!     Authors: Franck Lefevre, Francisco Gonzalez-Galindo
     6!     -------
    77!
    8 !     Version: 27/04/2017
     8!     Version: 14/11/2020
    99!
    1010!     ASIS scheme : for details on the method see
    … …  
    1414
    1515subroutine photochemistry(nlayer, nq, nesp, ionchem, nb_reaction_3_max,        &
    16                           nb_reaction_4_max, nb_phot_max, nphotion,            &
     16                          nb_reaction_4_max, nphot, nb_phot_max, nphotion,     &
    1717                          jonline, ig, lswitch, zycol, sza, ptimestep, press,  &
    1818                          alt, temp, temp_elect, dens, zmmean,                 &
    … …  
    3838integer, intent(in) :: nb_reaction_4_max
    3939                              ! number of bimolecular reactions
     40integer, intent(in) :: nphot
     41                              ! number of photodissociations
    4042integer, intent(in) :: nb_phot_max
    4143                              ! number of reactions treated numerically as photodissociations
    … …  
    121123integer,parameter :: i_elec    = 32
    122124
     125!Tracer indexes for photionization coeffs
     126
     127integer,parameter :: induv_co2 = 1
     128integer,parameter :: induv_o2  = 2
     129integer,parameter :: induv_o   = 3
     130integer,parameter :: induv_h2o = 4
     131integer,parameter :: induv_h2  = 5
     132integer,parameter :: induv_h2o2= 6
     133integer,parameter :: induv_o3  = 7
     134integer,parameter :: induv_n2  = 8
     135integer,parameter :: induv_n   = 9
     136integer,parameter :: induv_no  = 10
     137integer,parameter :: induv_co  = 11
     138integer,parameter :: induv_h   = 12
     139integer,parameter :: induv_no2 = 13
     140
    123141integer :: ilay
    124142
    … …  
    193211                             tau, sza, dist_sol, v_phot)
    194212
     213      !Calculation of photoionization rates, if needed
    195214      if (jionos .and. ionchem) then
    196215         call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday)
    197216         do ilay=1,lswitch-1
    198217            call phdisrate(ig,nlayer,2,sza,ilay)
    199          enddo
    200          v_phot(:,14)=jion(1,:,1)
    201          v_phot(:,15)=jion(1,:,2)
    202          v_phot(:,16)=jion(1,:,2)
    203          v_phot(:,17)=jion(1,:,3)
    204          v_phot(:,18)=jion(1,:,3)
    205          v_phot(:,19)=jion(1,:,4)
    206          v_phot(:,20)=jion(1,:,4)
    207          v_phot(:,21)=jion(2,:,1)
    208          v_phot(:,22)=jion(3,:,1)
    209          v_phot(:,23)=jion(10,:,1)
    210          v_phot(:,24)=jion(11,:,1)
    211          v_phot(:,25)=jion(11,:,2)
    212          v_phot(:,26)=jion(11,:,2)
    213          v_phot(:,27)=jion(8,:,1)
    214          v_phot(:,28)=jion(8,:,2)
    215          v_phot(:,29)=jion(8,:,2)
    216          v_phot(:,30)=jion(9,:,1)
    217          v_phot(:,31)=jion(12,:,1)
    218       endif
     218         end do
     219         !CO2 photoionization
     220         v_phot(:,nphot+ 1) = jion(induv_co2,:,1)
     221         v_phot(:,nphot+ 2) = jion(induv_co2,:,2)
     222         v_phot(:,nphot+ 3) = jion(induv_co2,:,2)
     223         v_phot(:,nphot+ 4) = jion(induv_co2,:,3)
     224         v_phot(:,nphot+ 5) = jion(induv_co2,:,3)
     225         v_phot(:,nphot+ 6) = jion(induv_co2,:,4)
     226         v_phot(:,nphot+ 7) = jion(induv_co2,:,4)
     227         !O2 photoionization
     228         v_phot(:,nphot+ 8) = jion(induv_o2,:,1)
     229         !O photoionization
     230         v_phot(:,nphot+ 9) = jion(induv_o,:,1)
     231         !NO photoionization
     232         v_phot(:,nphot+10) = jion(induv_no,:,1)
     233         !CO photoionization
     234         v_phot(:,nphot+11) = jion(induv_co,:,1)
     235         v_phot(:,nphot+12) = jion(induv_co,:,2)
     236         v_phot(:,nphot+13) = jion(induv_co,:,2)
     237         !N2 photoionization
     238         v_phot(:,nphot+14) = jion(induv_n2,:,1)
     239         v_phot(:,nphot+15) = jion(induv_n2,:,2)
     240         v_phot(:,nphot+16) = jion(induv_n2,:,2)
     241         !N photoionization
     242         v_phot(:,nphot+17) = jion(induv_n,:,1)
     243         !H photoionization
     244         v_phot(:,nphot+18) = jion(induv_h,:,1)
     245      end if
    219246     
    220247   else ! night
    … …  
    277304
    278305hetero_dust = .false.
    279 hetero_ice  = .true.
     306hetero_ice  = .false.
    280307
    281308call reactionrates(nlayer, ionchem, nb_reaction_3_max, nb_reaction_4_max, &
    … …  
    758785      c002(:) = 1.8e-11*exp(180./t(:))
    759786
     787!     c002(:) = c002(:)*1.12 ! FL li et al. 2007
     788
    760789!     robertson and smith, j. chem. phys. a 110, 6673, 2006
    761790
    … …  
    809838
    810839      c007(:) = 4.8e-11*exp(250./t(:))
     840
     841!     c007(:) = c007(:)*0.9 ! FL li et al. 2007
    811842
    812843      nb_reaction_4 = nb_reaction_4 + 1
    … …  
    851882
    852883      do ilev = 1,lswitch-1
     884!        ak0 = 3.1*2.4*4.4e-32*(t(ilev)/300.)**(-1.3) ! FL li et al 2017
    853885         ak0 = 2.4*4.4e-32*(t(ilev)/300.)**(-1.3)
    854886         ak1 = 7.5e-11*(t(ilev)/300.)**(0.2)
    … …  
    10581090!     jpl 2015
    10591091
     1092      do ilev = 1,lswitch-1
     1093
     1094!        branch 1 : oh + co -> h + co2
     1095
     1096         rate1 = 1.5e-13*(t(ilev)/300.)**(0.0)
     1097
     1098!        branch 2 : oh + co + m -> hoco + m
     1099
     1100         ak0 = 5.9e-33*(t(ilev)/300.)**(-1.0)
     1101         ak1 = 1.1e-12*(t(ilev)/300.)**(1.3)
     1102         rate2 = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
     1103         xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
     1104
     1105         e001(ilev) = rate1 + rate2*0.6**xpo
     1106      end do
     1107
     1108!     joshi et al., 2006
     1109
    10601110!     do ilev = 1,lswitch-1
    1061 
    1062 !        branch 1 : oh + co -> h + co2
    1063 
    1064 !        rate1 = 1.5e-13*(t(ilev)/300.)**(0.0)
    1065 
    1066 !        branch 2 : oh + co + m -> hoco + m
    1067 
    1068 !        ak0 = 5.9e-33*(t(ilev)/300.)**(-1.0)
    1069 !        ak1 = 1.1e-12*(t(ilev)/300.)**(1.3)
    1070 !        rate2 = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1)
    1071 !        xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2)
    1072 
    1073 !        e001(ilev) = rate1 + rate2*0.6**xpo
     1111!        k1a0 = 1.34*2.5*dens(ilev)                                  &
     1112!              *1/(1/(3.62e-26*t(ilev)**(-2.739)*exp(-20./t(ilev)))  &
     1113!              + 1/(6.48e-33*t(ilev)**(0.14)*exp(-57./t(ilev))))     ! typo in paper corrected
     1114!        k1b0 = 1.17e-19*t(ilev)**(2.053)*exp(139./t(ilev))          &
     1115!             + 9.56e-12*t(ilev)**(-0.664)*exp(-167./t(ilev))
     1116!        k1ainf = 1.52e-17*t(ilev)**(1.858)*exp(28.8/t(ilev))        &
     1117!               + 4.78e-8*t(ilev)**(-1.851)*exp(-318./t(ilev))
     1118!        x = k1a0/(k1ainf - k1b0)
     1119!        y = k1b0/(k1ainf - k1b0)
     1120!        fc = 0.628*exp(-1223./t(ilev)) + (1. - 0.628)*exp(-39./t(ilev))  &
     1121!           + exp(-t(ilev)/255.)
     1122!        fx = fc**(1./(1. + (alog(x))**2))                           ! typo in paper corrected
     1123!        k1a = k1a0*((1. + y)/(1. + x))*fx
     1124!        k1b = k1b0*(1./(1.+x))*fx
     1125!        e001(ilev) = k1a + k1b
    10741126!     end do
    1075 
    1076 !     joshi et al., 2006
    1077 
    1078       do ilev = 1,lswitch-1
    1079          k1a0 = 1.34*2.5*dens(ilev)                                  &
    1080                *1/(1/(3.62e-26*t(ilev)**(-2.739)*exp(-20./t(ilev)))  &
    1081                + 1/(6.48e-33*t(ilev)**(0.14)*exp(-57./t(ilev))))     ! typo in paper corrected
    1082          k1b0 = 1.17e-19*t(ilev)**(2.053)*exp(139./t(ilev))          &
    1083               + 9.56e-12*t(ilev)**(-0.664)*exp(-167./t(ilev))
    1084          k1ainf = 1.52e-17*t(ilev)**(1.858)*exp(28.8/t(ilev))        &
    1085                 + 4.78e-8*t(ilev)**(-1.851)*exp(-318./t(ilev))
    1086          x = k1a0/(k1ainf - k1b0)
    1087          y = k1b0/(k1ainf - k1b0)
    1088          fc = 0.628*exp(-1223./t(ilev)) + (1. - 0.628)*exp(-39./t(ilev))  &
    1089             + exp(-t(ilev)/255.)
    1090          fx = fc**(1./(1. + (alog(x))**2))                           ! typo in paper corrected
    1091          k1a = k1a0*((1. + y)/(1. + x))*fx
    1092          k1b = k1b0*(1./(1.+x))*fx
    1093          e001(ilev) = k1a + k1b
    1094       end do
    10951127
    10961128      nb_reaction_4 = nb_reaction_4 + 1
    … …  
    20162048
    20172049indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2d, 1.0, i_n)
     2050
     2051!===========================================================
     2052!      HDO + hv -> products
     2053!===========================================================
     2054
     2055!nb_phot = nb_phot + 1
     2056
     2057!indice_phot(nb_phot) = z3spec(1.0, i_hdo, 0.0, i_dummy, 0.0, i_dummy)
    20182058
    20192059!Only if ion chemistry included

  • trunk/LMDZ.MARS/libf/aeronomars/photolysis_mod.F90

    r2170 r2433  
    66
    77  integer, parameter :: nphot = 13        ! number of photolysis
     8! integer, parameter :: nphot = 14        ! number of photolysis (with hdo)
    89  integer, parameter :: nabs  = 10        ! number of absorbing gases
    910
    … …  
    3738  real, dimension(nw), save :: xsn2                                   ! n2 absorption cross-section (cm2)
    3839  real, dimension(nw), save :: yieldn2                                ! n2 photodissociation yield
     40  real, dimension(nw), save :: xshdo                                  ! hdo absorption cross-section (cm2)
    3941  real, dimension(nw), save :: albedo                                 ! surface albedo
    4042
    … …  
    9799
    98100  call rdxsn2(nw,wl,xsn2,yieldn2)
     101
     102! read and grid hdo cross-sections
     103
     104  call rdxshdo(nw,wl,xshdo)
    99105
    100106! set surface albedo
    … …  
    12471253!==============================================================================
    12481254
     1255      subroutine rdxshdo(nw, wl, yg)
     1256
     1257!-----------------------------------------------------------------------------*
     1258!=  PURPOSE:                                                                 =*
     1259!=  Read HDO molecular absorption cross section.  Re-grid data to match      =*
     1260!=  specified wavelength working grid.                                       =*
     1261!-----------------------------------------------------------------------------*
     1262!=  PARAMETERS:                                                              =*
     1263!=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
     1264!=           wavelength grid                                                 =*
     1265!=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
     1266!=           working wavelength grid                                         =*
     1267!=  YG     - REAL, molecular absoprtion cross section (cm^2) of HDO at    (O)=*
     1268!=           each specified wavelength                                       =*
     1269!-----------------------------------------------------------------------------*
     1270
     1271      use datafile_mod, only: datadir
     1272
     1273      IMPLICIT NONE
     1274
     1275!     input
     1276
     1277      integer :: nw               ! number of wavelength grid points
     1278      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
     1279
     1280!     output
     1281
     1282      real, dimension(nw) :: yg   ! hdo cross-sections (cm2)
     1283
     1284!     local
     1285
     1286      integer, parameter :: kdata = 900
     1287      real, parameter :: deltax = 1.e-4
     1288      REAL x1(kdata)
     1289      REAL y1(kdata)
     1290      INTEGER ierr
     1291      INTEGER i, n
     1292      CHARACTER*100 fil
     1293      integer :: kin, kout ! input/output logical units
     1294
     1295      kin = 10
     1296
     1297      fil = trim(datadir)//'/cross_sections/hdo_composite_295K.txt'
     1298      print*, 'section efficace HDO: ', fil
     1299
     1300      OPEN(UNIT=kin,FILE=fil,STATUS='old')
     1301
     1302      DO i = 1,17
     1303         read(kin,*)
     1304      END DO
     1305
     1306      n = 806
     1307      DO i = 1, n
     1308         READ(kin,*) x1(i), y1(i)
     1309      END DO
     1310      CLOSE (kin)
     1311
     1312      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
     1313      CALL addpnt(x1,y1,kdata,n,          0.,0.)
     1314      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
     1315      CALL addpnt(x1,y1,kdata,n,      1.e+38,0.)
     1316      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
     1317      IF (ierr .NE. 0) THEN
     1318         WRITE(*,*) ierr, fil
     1319         STOP
     1320      ENDIF
     1321     
     1322      end subroutine rdxshdo
     1323
     1324!==============================================================================
     1325
    12491326      subroutine rdxsh2o2(nw, wl, xsh2o2)
    12501327

  • trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F

    r2170 r2433  
    6060
    6161      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o,
    62      $           j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2
     62     $           j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2, j_hdo
    6363
    6464      integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no,
    … …  
    9696      j_no2     = 12     ! no2 + hv    -> no + o
    9797      j_n2      = 13     ! n2 + hv     -> n + n
     98      j_hdo     = 14     ! hdo + hv    -> products
    9899
    99100!==== air column increments and rayleigh optical depth
    … …  
    167168            sj(ilay,iw,j_no)  = xsno(iw)*yieldno(iw)     ! no
    168169            sj(ilay,iw,j_n2)  = xsn2(iw)*yieldn2(iw)     ! n2
     170!           sj(ilay,iw,j_hdo) = xshdo(iw)                ! hdo
    169171         end do
    170172      end do
Note: See TracChangeset for help on using the changeset viewer.