Changeset 2795 for trunk/LMDZ.VENUS


Ignore:
Timestamp:
Sep 14, 2022, 1:40:00 PM (2 years ago)
Author:
flefevre
Message:

VENUS:

1) Ajout de la chimie de l'azote. 4 nouvelles especes: n, n2d, no, no2

Photolyse online (jonline = .true.) obligatoire pour le moment.

2) Mise à jour de la vitesse de CO + OH -> CO2 + H

Location:
trunk/LMDZ.VENUS/libf/phyvenus
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.VENUS/libf/phyvenus/chemparam_mod.F90

    r2464 r2795  
    1313                 i_cocl2, i_s, i_so, i_so2, i_so3,       &
    1414                 i_s2o2, i_ocs, i_hso3, i_h2so4, i_s2,   &
    15                  i_clso2, i_oscl, i_n2, i_he
     15                 i_clso2, i_oscl, i_n2, i_he, i_no2,     &
     16                 i_no, i_n, i_n2d
    1617                 
    1718INTEGER, SAVE :: i_h2oliq, i_h2so4liq
    … …  
    2728REAL, DIMENSION(:), SAVE, ALLOCATABLE :: M_tr
    2829 
    29 
     30REAL, DIMENSION(:,:),SAVE, ALLOCATABLE :: no_emission
     31REAL, DIMENSION(:,:),SAVE, ALLOCATABLE :: o2_emission
    3032!----------------------------------------------------------------------------
    3133! DEF FOR CL_SCHEME = 1 (AURELIEN)
    … …  
    837839                i_he=i
    838840                M_tr(i_he)=4.0026
     841                CASE('no2')
     842                i_no2=i
     843                M_tr(i_no2)=46.0055
     844                CASE('no')
     845                i_no = i
     846                M_tr(i_no)=30.0061
     847                CASE('n')
     848                i_n = i
     849                M_tr(i_n)=14.0067
     850                CASE('n2d')
     851                i_n2d = i
     852                M_tr(i_n2d)=14.0067
    839853! MICROPHYSICAL TRACERS FOR CL_SCHEME=1
    840854                CASE('h2oliq')

  • trunk/LMDZ.VENUS/libf/phyvenus/moldiff_red.F90

    r2791 r2795  
    6060      real*8,dimension(:),allocatable,save ::  wi,Wad,Uthermal,Lambdaexo,Hspecie
    6161      real*8,dimension(:),allocatable,save :: Mtot1,Mtot2,Mraf1,Mraf2
    62       integer,parameter :: nb_esp_diff=15
     62      integer,parameter :: nb_esp_diff=16
    6363      character(len=20),dimension(nb_esp_diff) :: ListeDiff
    6464!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    … …  
    125125        ListeDiff(14)='n'
    126126        ListeDiff(15)='he'
     127        ListeDiff(16)='n2d'
    127128        i_esc_h=1000
    128129        i_esc_h2=1000

  • trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90

    r2780 r2795  
    1 subroutine photochemistry_venus(nz, n_lon, zlocal, ptimestep, p, t, tr, mumean, sza_input, lon, lat, nesp, iter, prod_tr, loss_tr)
     1subroutine photochemistry_venus(nz, n_lon, zlocal, ptimestep, p, t, tr, mumean, sza_input, lon, lat, nesp, iter, prod_tr, loss_tr, em_no, em_o2)
    22
    33use chemparam_mod
    … …  
    3030real, dimension(nz,nesp) :: prod_tr ! production (cm-3.s-1)
    3131real, dimension(nz,nesp) :: loss_tr ! loss       (cm-3.s-1)
     32real, dimension (nz) :: em_no       ! volume emission rate of no
     33real, dimension (nz) :: em_o2       ! volume emission rate of o2(deltag)
    3234
    3335!===================================================================
    … …  
    3840
    3941!===================================================================
    40 !     local:
     42!     local: 
    4143!===================================================================
    4244
    … …  
    4446! true : on-line calculation of photodissociation rates ! false : lookup table
    4547
    46 logical, save :: jonline = .false.
     48logical, save :: jonline = .true.
    4749
    4850logical, save :: firstcall = .true.
    … …  
    5153real, dimension(nz)  :: surfice1d, surfdust1d
    5254
     55integer :: ind_norec
     56integer :: ind_orec
     57
    5358! photolysis lookup table (case jonline = .false.)
    5459! if prior to jvenus.20211025, set nztable = 201 below
    5560
    56 integer, parameter :: nj = 19, nztable = 281, nsza = 27, nso2 = 13
     61integer, parameter :: nj = 22, nztable = 281, nsza = 27, nso2 = 13
    5762real, dimension(nso2,nsza,nztable,nj), save :: jphot
    5863real, dimension(nztable), save :: table_colair
    … …  
    7782! reaction rates
    7883
    79 integer, parameter :: nb_phot_max       = 30
     84integer, parameter :: nb_phot_max       = 35
     85integer, parameter :: nb_phot           = 22
    8086integer, parameter :: nb_reaction_3_max = 12
    81 integer, parameter :: nb_reaction_4_max = 87
     87integer, parameter :: nb_reaction_4_max = 98
    8288
    8389real, dimension(nz,nb_phot_max)       :: v_phot
    … …  
    104110if (firstcall) then
    105111!===================================================================
    106 !     initialisation of photolysis
     112!     initialisation of  photolysis
    107113!===================================================================
    108114
    109115   if (jonline) then
    110       print*, 'Photochemistry: Read UV absorption cross-sections:'
     116      print*, 'photochemistry: Read UV absorption cross-sections:'
    111117      call init_photolysis
    112118   else
    113       print*, 'Photochemistry: Read photolysis lookup table:'
     119      print*, 'photochemistry: Read photolysis lookup table:'
    114120      call init_chimie(nj, nztable, nsza, nso2, jphot, table_colair, table_colso2, table_sza)
    115121   end if
    … …  
    135141do iz = 1,nz
    136142   conc(iz) = p(iz)/(1.38E-19*t(iz))
    137    c(iz,:) = tr(iz,:)*conc(iz) 
     143   c(iz,:) = tr(iz,:)*conc(iz)
    138144end do
    139145     
    … …  
    147153
    148154if (jonline) then
    149    if (sza_input <= 95.) then !day at 30 km
     155   if (sza_input <= 95.) then ! day at 30 km
    150156      call photolysis_online(nz, nb_phot_max, zlocal, p,                 &
    151157                             t, mumean, i_co2,i_co, i_o, i_o1d,          &
    … …  
    153159                             i_h, i_hcl, i_cl2, i_hocl, i_so2, i_so,     &
    154160                             i_so3, i_clo, i_ocs, i_cocl2, i_h2so4, i_cl,&
     161                             i_no2, i_no, i_n2, i_n2d,                   &
    155162                             nesp, tr, sza_input, dist_sol, v_phot)
    156163   else ! night
    … …  
    158165   end if
    159166else
    160    call phot(nj, nztable, nsza, nso2, sza_input, dist_sol, mumean, tr(:,i_co2), tr(:,i_so2),      &
     167   call phot(nj, nztable, nsza, nso2, sza_input, dist_sol, mumean, tr(:,i_co2), tr(:,i_so2),         &
    161168             jphot, table_colair, table_colso2, table_sza, nz, nb_phot_max, t, p, v_phot)
    162169end if
    … …  
    166173!===================================================================
    167174                   
    168 call krates(hetero_ice, hetero_dust, nz, nesp, nj, c, conc, t, p, nb_phot_max, nb_reaction_3_max, &
    169             nb_reaction_4_max, v_3, v_4, v_phot, sza_input)
     175call krates(hetero_ice, hetero_dust, nz, nesp, nb_phot, c, conc, t, p, nb_phot_max, nb_reaction_3_max, &
     176            nb_reaction_4_max, v_3, v_4, v_phot, sza_input, ind_norec, ind_orec)
    170177
    171178!===================================================================
    … …  
    252259
    253260   tr(iz,:)  = max(c(iz,:)/conc(iz), 1.e-30)
    254    
    255 !  save prod and loss espece
     261
     262!  save production and loss terms (for diagnostic only)
    256263   
    257264   prod_tr(iz,:) = prod(:)
    … …  
    259266                       
    260267end do  ! end of loop over vertical levels
    261    
    262 !==================
    263 !!!!! MODEL 1D !!!! ==> n_lon = 1 !!!!
    264 !==================
    265 
    266 !IF(n_lon .EQ. 1) THEN
    267 !PRINT*,'On est en 1D'
    268 !PRINT*,"DEBUT rate_save"
    269 !CALL rate_save(nz,p(:),t(:),tr(:,:),nesp,v_phot(:,:),v_3(:,:),v_4(:,:))
    270 !PRINT*,"FIN rate_save"
    271 !END IF
    272      
     268 
     269! no and o2(delta) emissions
     270
     271em_no(:) = c(:,i_o)*c(:,i_n)*v_4(:,ind_norec)   !2.8e-17*(300./temp(:)))**0.5
     272em_o2(:) = c(:,i_o2dg)/4470.    ! 4470 s : lafferty et al., 1998
     273
    273274end subroutine photochemistry_venus
    274275
    … …  
    513514
    514515!===========================================================
     516!      NO2 + hv -> NO + O
     517!===========================================================
     518
     519nb_phot = nb_phot + 1
     520
     521indice_phot(nb_phot) = z3spec(1.0, i_no2, 1.0, i_no, 1.0, i_o)
     522
     523!===========================================================
     524!      NO + hv -> N + O
     525!===========================================================
     526
     527nb_phot = nb_phot + 1
     528
     529indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_n, 1.0, i_o)
     530
     531!===========================================================
     532!      N2 + hv -> N(2D) + N
     533!===========================================================
     534
     535nb_phot = nb_phot + 1
     536
     537indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n, 1.0, i_n2d)
     538
     539!===========================================================
    515540!      a001 : O + O2 + CO2 -> O3 + CO2
    516541!===========================================================
    … …  
    727752
    728753indice_3(nb_reaction_3) = z3spec(2.0, i_h, 1.0, i_h2, 0.0, i_dummy)
     754
     755!===========================================================
     756!      d001 : NO2 + O -> NO + O2
     757!===========================================================
     758
     759nb_reaction_4 = nb_reaction_4 + 1
     760
     761indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_o, 1.0, i_no, 1.0, i_o2)
     762
     763!===========================================================
     764!      d002 : NO + O3 -> NO2 + O2
     765!===========================================================
     766
     767nb_reaction_4 = nb_reaction_4 + 1
     768
     769indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_o3, 1.0, i_no2, 1.0, i_o2)
     770
     771!===========================================================
     772!      d003 : NO + HO2 -> NO2 + OH
     773!===========================================================
     774
     775nb_reaction_4 = nb_reaction_4 + 1
     776
     777indice_4(nb_reaction_4) = z4spec(1.0, i_no, 1.0, i_ho2, 1.0, i_no2, 1.0, i_oh)
     778
     779!===========================================================
     780!      d004 : N + NO -> N2 + O
     781!===========================================================
     782
     783nb_reaction_4 = nb_reaction_4 + 1
     784
     785indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_no, 1.0, i_n2, 1.0, i_o)
     786
     787!===========================================================
     788!      d005 : N + O2 -> NO + O
     789!===========================================================
     790
     791nb_reaction_4 = nb_reaction_4 + 1
     792
     793indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o2, 1.0, i_no, 1.0, i_o)
     794
     795!===========================================================
     796!      d006 : NO2 + H -> NO + OH
     797!===========================================================
     798
     799nb_reaction_4 = nb_reaction_4 + 1
     800
     801indice_4(nb_reaction_4) = z4spec(1.0, i_no2, 1.0, i_h, 1.0, i_no, 1.0, i_oh)
     802
     803!===========================================================
     804!      d007 : N + O -> NO
     805!===========================================================
     806
     807nb_reaction_4 = nb_reaction_4 + 1
     808
     809indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy)
     810
     811!===========================================================
     812!      d008 : N + HO2 -> NO + OH
     813!===========================================================
     814
     815nb_reaction_4 = nb_reaction_4 + 1
     816
     817indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_ho2, 1.0, i_no, 1.0, i_oh)
     818
     819!===========================================================
     820!      d009 : N + OH -> NO + H
     821!===========================================================
     822
     823nb_reaction_4 = nb_reaction_4 + 1
     824
     825indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_oh, 1.0, i_no, 1.0, i_h)
     826
     827!===========================================================
     828!      d010 : N(2D) + O -> N + O
     829!===========================================================
     830
     831nb_phot = nb_phot + 1
     832
     833indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
     834
     835!===========================================================
     836!      d011 : N(2D) + N2 -> N + N2
     837!===========================================================
     838
     839nb_phot = nb_phot + 1
     840
     841indice_phot(nb_phot) = z3spec(1.0, i_n2d, 1.0, i_n, 0.0, i_dummy)
     842
     843!===========================================================
     844!      d012 : N(2D) + CO2 -> NO + CO
     845!===========================================================
     846
     847nb_reaction_4 = nb_reaction_4 + 1
     848
     849indice_4(nb_reaction_4) = z4spec(1.0, i_n2d, 1.0, i_co2, 1.0, i_no, 1.0, i_co)
     850
     851!===========================================================
     852!      d013 : N + O + CO2 -> NO + CO2
     853!===========================================================
     854
     855nb_reaction_4 = nb_reaction_4 + 1
     856
     857indice_4(nb_reaction_4) = z4spec(1.0, i_n, 1.0, i_o, 1.0, i_no, 0.0, i_dummy)
    729858
    730859!===========================================================
    … …  
    13741503
    13751504!===========================================================
    1376 !      i001: O2(Dg) + CO2 -> O2 + CO2 + hv
     1505!      i001: O2(Dg) + CO2 -> O2 + CO2
    13771506!===========================================================
    13781507
    … …  
    14011530    (nb_reaction_3 /= nb_reaction_3_max) .or.  &
    14021531    (nb_reaction_4 /= nb_reaction_4_max)) then
    1403   print*, 'wrong dimensions in indice'
    1404   stop
     1532   print*, 'wrong dimensions in indice'
     1533   stop
    14051534end if 
    14061535
    … …  
    14381567integer :: indcol, indsza, indso2
    14391568integer :: isza, iz, i, iso2, ij
    1440 
    14411569integer :: nb_phot_max
    14421570real, dimension(nz,nb_phot_max), INTENT(INOUT) :: v_phot
    … …  
    15761704!    18    cocl2 + hv  -> cl + cl + co
    15771705!    19    h2so4 + hv  -> so3 + h2o
     1706!    20    no2 + hv    -> no + o
     1707!    21    no + hv     -> n + o 
     1708!    22    n2 + hv     -> n + n
    15781709
    15791710! fill v_phot array
    … …  
    16221753!======================================================================
    16231754
    1624  subroutine krates(hetero_ice,hetero_dust, nz, nesp, nj, c, conc, t, p, nb_phot_max, nb_reaction_3_max, nb_reaction_4_max, v_3, v_4, v_phot,sza_input)
     1755 subroutine krates(hetero_ice,hetero_dust, nz, nesp, nj, c, conc, t, p, nb_phot_max, nb_reaction_3_max, nb_reaction_4_max, v_3, v_4, v_phot,sza_input, ind_norec, ind_orec)
    16251756 
    16261757!================================================================
    … …  
    16441775integer :: iz
    16451776logical, INTENT(IN) :: hetero_ice, hetero_dust
    1646 real    :: ak0, ak1, xpo, rate, pi, gam, bid
     1777integer :: ind_norec, ind_orec
     1778real    :: ak0, ak1, xpo, rate, rate1, rate2, pi, gam
    16471779real    :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y
    16481780real, dimension(nz)  :: surfice1d, surfdust1d
    … …  
    16541786                       c008, c009, c010, c011, c012, c013, c014,   &
    16551787                       c015, c016, c017, c018,                     &
    1656                        d001, d002, d003,                           &
     1788                       d001, d002, d003, d004, d005, d006, d007,   &
     1789                       d008, d009, d010, d011, d012, d013,         &
    16571790                       e001, e002, e003, e004, e005, e006, e007,   &
    16581791                       e008, e009, e010, e011, e012, e013, e014,   &
    … …  
    17161849      v_3(:,nb_reaction_3) = a002(:)
    17171850
     1851      ind_orec = nb_reaction_3
     1852
    17181853!---  a003: o + o3 -> o2 + o2
    17191854
    … …  
    19832118!     jpl 2006
    19842119
    1985       d001(:) = 5.1E-12*exp(210./t(:))
     2120      d001(:) = 5.1e-12*exp(210./t(:))
     2121
     2122      nb_reaction_4 = nb_reaction_4 + 1
     2123      v_4(:,nb_reaction_4) = d001(:)
    19862124
    19872125!---  d002: no + o3 -> no2 + o2
    … …  
    19892127!     jpl 2006
    19902128
    1991       d002(:) = 3.0E-12*exp(-1500./t(:))
     2129      d002(:) = 3.0e-12*exp(-1500./t(:))
     2130
     2131      nb_reaction_4 = nb_reaction_4 + 1
     2132      v_4(:,nb_reaction_4) = d002(:)
    19922133
    19932134!---  d003: no + ho2 -> no2 + oh
    19942135
    1995 !     jpl 2006
    1996 
    1997       d003(:) = 3.5E-12*exp(250./t(:))
     2136!     jpl 2011
     2137
     2138      d003(:) = 3.3e-12*exp(270./t(:))
     2139
     2140      nb_reaction_4 = nb_reaction_4 + 1
     2141      v_4(:,nb_reaction_4) = d003(:)
     2142
     2143
     2144!---  d004: n + no -> n2 + o
     2145
     2146!     jpl 2011
     2147
     2148      d004(:) = 2.1e-11*exp(100./t(:))
     2149
     2150      nb_reaction_4 = nb_reaction_4 + 1
     2151      v_4(:,nb_reaction_4) = d004(:)
     2152
     2153!---  d005: n + o2 -> no + o
     2154
     2155!     jpl 2011
     2156
     2157      d005(:) = 1.5e-11*exp(-3600./t(:))
     2158
     2159      nb_reaction_4 = nb_reaction_4 + 1
     2160      v_4(:,nb_reaction_4) = d005(:)
     2161
     2162!---  d006: no2 + h -> no + oh
     2163
     2164!     jpl 2011
     2165
     2166      d006(:) = 4.0e-10*exp(-340./t(:))
     2167
     2168      nb_reaction_4 = nb_reaction_4 + 1
     2169      v_4(:,nb_reaction_4) = d006(:)
     2170
     2171!---  d007: n + o -> no
     2172
     2173      d007(:) = 2.8e-17*(300./t(:))**0.5
     2174
     2175      nb_reaction_4 = nb_reaction_4 + 1
     2176      v_4(:,nb_reaction_4) = d007(:)
     2177
     2178      ind_norec = nb_reaction_4
     2179
     2180!---  d008: n + ho2 -> no + oh
     2181
     2182!     brune et al., j. chem. phys., 87, 1983
     2183
     2184      d008(:) = 2.19e-11
     2185
     2186      nb_reaction_4 = nb_reaction_4 + 1
     2187      v_4(:,nb_reaction_4) = d008(:)
     2188
     2189!---  d009: n + oh -> no + h
     2190
     2191!     atkinson et al., j. phys. chem. ref. data, 18, 881, 1989
     2192
     2193      d009(:) = 3.8e-11*exp(85./t(:))
     2194
     2195      nb_reaction_4 = nb_reaction_4 + 1
     2196      v_4(:,nb_reaction_4) = d009(:)
     2197
     2198!---  d010: n2d + o -> n + o
     2199
     2200!     herron, j. phys. chem. ref. data, 1999
     2201
     2202      d010(:) = 3.3e-12*exp(-260./t(:))
     2203
     2204      nb_phot = nb_phot + 1
     2205      v_phot(:,nb_phot) = d010(:)*c(:,i_o)
     2206
     2207!---  d011: n2d + n2 -> n + n2
     2208
     2209!     herron, j. phys. chem. ref. data, 1999
     2210
     2211      d011(:) = 1.7e-14
     2212
     2213      nb_phot = nb_phot + 1
     2214      v_phot(:,nb_phot) = d011(:)*c(:,i_n2)
     2215
     2216!---  d012: n2d + co2 -> no + co
     2217
     2218!     herron, j. phys. chem. ref. data, 1999
     2219
     2220      d012(:) = 3.6e-13
     2221
     2222      nb_reaction_4 = nb_reaction_4 + 1
     2223      v_4(:,nb_reaction_4) = d012(:)
     2224
     2225!---  d013: n + o + co2 -> no + co2
     2226     
     2227      d013(:) = 1.83e-32 * (298./t(:))**0.5
     2228     
     2229      nb_reaction_4 = nb_reaction_4 + 1
     2230      v_4(:,nb_reaction_4) = d013(:)
    19982231
    19992232!----------------------------------------------------------------------
    … …  
    20032236!---  e001: oh + co -> co2 + h
    20042237
    2005 !     jpl 2003
    2006 
    2007 !     e001(:) = 1.5E-13*(1 + 0.6*p(:)/1013.)
    2008 
    2009 !     mccabe et al., grl, 28, 3135, 2001
    2010 
    2011 !     e001(:) = 1.57E-13 + 3.54E-33*conc(:)
    2012 
    2013 !     jpl 2006
    2014 
    2015 !     ak0 = 1.5E-13*(t(:)/300.)**(0.6)
    2016 !     ak1 = 2.1E-9*(t(:)/300.)**(6.1)
    2017 !     rate1 = ak0/(1. + ak0/(ak1/conc(:)))
    2018 !     xpo1 = 1./(1. + alog10(ak0/(ak1/conc(:)))**2)
    2019 
    2020 !     ak0 = 5.9E-33*(t(:)/300.)**(-1.4)
    2021 !     ak1 = 1.1E-12*(t(:)/300.)**(1.3)
    2022 !     rate2 = (ak0*conc(:))/(1. + ak0*conc(:)/ak1)
    2023 !     xpo2 = 1./(1. + alog10((ak0*conc(:))/ak1)**2)
    2024 
    2025 !     e001(:) = rate1*0.6**xpo1 + rate2*0.6**xpo2
     2238!     jpl 2015
     2239
     2240      do iz = 1,nz
     2241
     2242!        branch 1 : oh + co -> h + co2
     2243
     2244         rate1 = 1.5e-13*(t(iz)/300.)**(0.0)
     2245
     2246!        branch 2 : oh + co + m -> hoco + m
     2247
     2248         ak0 = 5.9e-33*(t(iz)/300.)**(-1.0)
     2249         ak1 = 1.1e-12*(t(iz)/300.)**(1.3)
     2250         rate2 = (ak0*conc(iz))/(1. + ak0*conc(iz)/ak1)
     2251         xpo = 1./(1. + alog10((ak0*conc(iz))/ak1)**2)
     2252
     2253         e001(iz) = rate1 + rate2*0.6**xpo
     2254      end do
    20262255
    20272256!     joshi et al., 2006
    20282257
    2029       do iz = 1,nz
    2030          k1a0 = 1.34*2.5*conc(iz)                                &
    2031                *1/(1/(3.62E-26*t(iz)**(-2.739)*exp(-20./t(iz)))  &
    2032                + 1/(6.48E-33*t(iz)**(0.14)*exp(-57./t(iz))))     ! corrige de l'erreur publi
    2033          k1b0 = 1.17E-19*t(iz)**(2.053)*exp(139./t(iz))          &
    2034               + 9.56E-12*t(iz)**(-0.664)*exp(-167./t(iz))
    2035          k1ainf = 1.52E-17*t(iz)**(1.858)*exp(28.8/t(iz))        &
    2036                 + 4.78E-8*t(iz)**(-1.851)*exp(-318./t(iz))
    2037          x = k1a0/(k1ainf - k1b0)
    2038          y = k1b0/(k1ainf - k1b0)
    2039          fc = 0.628*exp(-1223./t(iz)) + (1. - 0.628)*exp(-39./t(iz))  &
    2040             + exp(-t(iz)/255.)
    2041          fx = fc**(1./(1. + (alog(x))**2))                   ! corrige de l'erreur publi
    2042          k1a = k1a0*((1. + y)/(1. + x))*fx
    2043          k1b = k1b0*(1./(1.+x))*fx
    2044 
    2045          e001(iz) = k1a + k1b
    2046       end do
     2258!     do iz = 1,nz
     2259!        k1a0 = 1.34*2.5*conc(iz)                                &
     2260!              *1/(1/(3.62E-26*t(iz)**(-2.739)*exp(-20./t(iz)))  &
     2261!              + 1/(6.48E-33*t(iz)**(0.14)*exp(-57./t(iz))))     ! typo in paper corrected
     2262!        k1b0 = 1.17E-19*t(iz)**(2.053)*exp(139./t(iz))          &
     2263!             + 9.56E-12*t(iz)**(-0.664)*exp(-167./t(iz))
     2264!        k1ainf = 1.52E-17*t(iz)**(1.858)*exp(28.8/t(iz))        &
     2265!               + 4.78E-8*t(iz)**(-1.851)*exp(-318./t(iz))
     2266!        x = k1a0/(k1ainf - k1b0)
     2267!        y = k1b0/(k1ainf - k1b0)
     2268!        fc = 0.628*exp(-1223./t(iz)) + (1. - 0.628)*exp(-39./t(iz))  &
     2269!           + exp(-t(iz)/255.)
     2270!        fx = fc**(1./(1. + (alog(x))**2))                       ! typo in paper corrected
     2271!        k1a = k1a0*((1. + y)/(1. + x))*fx
     2272!        k1b = k1b0*(1./(1.+x))*fx
     2273!        e001(iz) = k1a + k1b
     2274!     end do
    20472275
    20482276      nb_reaction_4 = nb_reaction_4 + 1
    … …  
    29783206! or reactions a + c -> b + c
    29793207! or reactions a + ice -> b + c
    2980 
    29813208do iphot = 1,nb_phot_max
    29823209
    … …  
    44704697                i_v=i_v+1 
    44714698!===============================
    4472 !    30    o2(Dg) + CO2 -> O2 + CO2 + hv
     4699!    30    o2(Dg) + CO2 -> O2 + CO2
    44734700!===============================
    44744701            rate_local = vphot(i_lev,i_v)*traceur(i_lev,i_o2dg)*concentration(i_lev)

  • trunk/LMDZ.VENUS/libf/phyvenus/photolysis_mod.F90

    r2780 r2795  
    55! photolysis
    66
    7   integer, save :: nphot = 19              ! number of photolysis
     7  integer, save :: nphot = 22             ! number of photolysis
    88
    99!$OMP THREADPRIVATE(nphot)
    1010
    11   integer, parameter :: nabs  = 16        ! number of absorbing gases
     11  integer, parameter :: nabs  = 19        ! number of absorbing gases
    1212
    1313! spectral grid
    … …  
    4646  real, dimension(nw), save :: xscocl2                                ! cocl2 absorption cross-section (cm2)
    4747  real, dimension(nw), save :: xsh2so4                                ! h2so4 absorption cross-section (cm2)
    48 !  real, dimension(nw), save :: xsh2                                   ! h2 absorption cross-section (cm2)
    49 !  real, dimension(nw), save :: yieldh2                                ! h2 photodissociation yield
    50 !  real, dimension(nw), save :: xsno2, xsno2_220, xsno2_294            ! no2 absorption cross-section at 220-294 k (cm2)
    51 !  real, dimension(nw), save :: yldno2_248, yldno2_298                 ! no2 quantum yield at 248-298 k
    52 !  real, dimension(nw), save :: xsno                                   ! no absorption cross-section (cm2)
    53 !  real, dimension(nw), save :: yieldno                                ! no photodissociation yield
    54 !  real, dimension(nw), save :: xsn2                                   ! n2 absorption cross-section (cm2)
    55 !  real, dimension(nw), save :: yieldn2                                ! n2 photodissociation yield
    56 !  real, dimension(nw), save :: xshdo                                  ! hdo absorption cross-section (cm2)
     48  real, dimension(nw), save :: xsh2                                   ! h2 absorption cross-section (cm2)
     49  real, dimension(nw), save :: yieldh2                                ! h2 photodissociation yield
     50  real, dimension(nw), save :: xsno2, xsno2_220, xsno2_294            ! no2 absorption cross-section at 220-294 k (cm2)
     51  real, dimension(nw), save :: yldno2_248, yldno2_298                 ! no2 quantum yield at 248-298 k
     52  real, dimension(nw), save :: xsno                                   ! no absorption cross-section (cm2)
     53  real, dimension(nw), save :: yieldno                                ! no photodissociation yield
     54  real, dimension(nw), save :: xsn2                                   ! n2 absorption cross-section (cm2)
     55  real, dimension(nw), save :: yieldn2                                ! n2 photodissociation yield
     56  real, dimension(nw), save :: xshdo                                  ! hdo absorption cross-section (cm2)
    5757  real, dimension(nw), save :: albedo                                 ! surface albedo
    5858
    … …  
    146146! read and grid h2 cross-sections
    147147
    148 !  call rdxsh2(nw,wl,wc,xsh2,yieldh2)
     148  call rdxsh2(nw,wl,wc,xsh2,yieldh2)
    149149
    150150! read and grid no cross-sections
    151151
    152 !  call rdxsno(nw,wl,xsno,yieldno)
     152  call rdxsno(nw,wl,xsno,yieldno)
    153153
    154154! read and grid no2 cross-sections
    155155
    156 !  call rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248,yldno2_298)
     156  call rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248,yldno2_298)
    157157
    158158! read and grid n2 cross-sections
    159159
    160 !  call rdxsn2(nw,wl,xsn2,yieldn2)
     160  call rdxsn2(nw,wl,xsn2,yieldn2)
    161161
    162162! read and grid hdo cross-sections
    163163
    164 !  call rdxshdo(nw,wl,xshdo)
     164  call rdxshdo(nw,wl,xshdo)
    165165
    166166! set surface albedo
    … …  
    15431543!==============================================================================
    15441544
    1545 !      subroutine rdxshdo(nw, wl, yg)
     1545       subroutine rdxshdo(nw, wl, yg)
    15461546
    15471547!-----------------------------------------------------------------------------*
    … …  
    15591559!-----------------------------------------------------------------------------*
    15601560
    1561 !      USE mod_phys_lmdz_para, ONLY: is_master
    1562 !      USE mod_phys_lmdz_transfert_para, ONLY: bcast
    1563 
    1564 !      IMPLICIT NONE
     1561      USE mod_phys_lmdz_para, ONLY: is_master
     1562      USE mod_phys_lmdz_transfert_para, ONLY: bcast
     1563
     1564      IMPLICIT NONE
    15651565
    15661566!     input
    15671567
    1568 !      integer :: nw               ! number of wavelength grid points
    1569 !      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
     1568      integer :: nw               ! number of wavelength grid points
     1569      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
    15701570
    15711571!     output
    15721572
    1573 !      real, dimension(nw) :: yg   ! hdo cross-sections (cm2)
     1573      real, dimension(nw) :: yg   ! hdo cross-sections (cm2)
    15741574
    15751575!     local
    15761576
    1577 !      integer, parameter :: kdata = 900
    1578 !      real, parameter :: deltax = 1.e-4
    1579 !      REAL x1(kdata)
    1580 !      REAL y1(kdata)
    1581 !      INTEGER ierr
    1582 !      INTEGER i, n
    1583 !      CHARACTER*100 fil
    1584 !      integer :: kin, kout ! input/output logical units
    1585 
    1586 !      kin = 10
    1587 
    1588 !      fil = '/cross_sections/hdo_composite_295K.txt'
    1589 !      print*, 'section efficace HDO: ', fil
    1590 
    1591 !      if(is_master) then
    1592 
    1593 !      OPEN(UNIT=kin,FILE=fil,STATUS='old')
    1594 
    1595 !      DO i = 1,17
    1596 !         read(kin,*)
    1597 !      END DO
    1598 
    1599 !      n = 806
    1600 !      DO i = 1, n
    1601 !         READ(kin,*) x1(i), y1(i)
    1602 !      END DO
    1603 !      CLOSE (kin)
    1604 
    1605 !      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
    1606 !      CALL addpnt(x1,y1,kdata,n,          0.,0.)
    1607 !      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
    1608 !      CALL addpnt(x1,y1,kdata,n,      1.e+38,0.)
    1609 !      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
    1610 !      IF (ierr .NE. 0) THEN
    1611 !         WRITE(*,*) ierr, fil
    1612 !         STOP
    1613 !      ENDIF
    1614 
    1615 !      endif !is_master
    1616 
    1617 !      call bcast(yg)
     1577      integer, parameter :: kdata = 900
     1578      real, parameter :: deltax = 1.e-4
     1579      REAL x1(kdata)
     1580      REAL y1(kdata)
     1581      INTEGER ierr
     1582      INTEGER i, n
     1583      CHARACTER*100 fil
     1584      integer :: kin, kout ! input/output logical units
     1585
     1586      kin = 10
     1587
     1588      fil = '/cross_sections/hdo_composite_295K.txt'
     1589      print*, 'section efficace HDO: ', fil
     1590
     1591      if(is_master) then
     1592
     1593      OPEN(UNIT=kin,FILE=fil,STATUS='old')
     1594
     1595      DO i = 1,17
     1596         read(kin,*)
     1597      END DO
     1598
     1599      n = 806
     1600      DO i = 1, n
     1601         READ(kin,*) x1(i), y1(i)
     1602      END DO
     1603      CLOSE (kin)
     1604
     1605      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
     1606      CALL addpnt(x1,y1,kdata,n,          0.,0.)
     1607      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
     1608      CALL addpnt(x1,y1,kdata,n,      1.e+38,0.)
     1609      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
     1610      IF (ierr .NE. 0) THEN
     1611         WRITE(*,*) ierr, fil
     1612         STOP
     1613      ENDIF
     1614
     1615      endif !is_master
     1616
     1617      call bcast(yg)
    16181618     
    1619 !      end subroutine rdxshdo
     1619      end subroutine rdxshdo
    16201620
    16211621!==============================================================================
    … …  
    26242624!==============================================================================
    26252625
    2626 !      subroutine rdxsh2(nw, wl, wc, yg, yieldh2)
     2626       subroutine rdxsh2(nw, wl, wc, yg, yieldh2)
    26272627
    26282628!-----------------------------------------------------------------------------*
    … …  
    26352635!-----------------------------------------------------------------------------*
    26362636
    2637 !      USE mod_phys_lmdz_para, ONLY: is_master
    2638 !      USE mod_phys_lmdz_transfert_para, ONLY: bcast
    2639 
    2640 !      implicit none
     2637      USE mod_phys_lmdz_para, ONLY: is_master
     2638      USE mod_phys_lmdz_transfert_para, ONLY: bcast
     2639
     2640      implicit none
    26412641
    26422642!     input
    26432643
    2644 !      integer :: nw                   ! number of wavelength grid points
    2645 !      real, dimension(nw) :: wl, wc   ! lower and central wavelength for each interval
     2644      integer :: nw                   ! number of wavelength grid points
     2645      real, dimension(nw) :: wl, wc   ! lower and central wavelength for each interval
    26462646
    26472647!     output
    26482648
    2649 !      real, dimension(nw) :: yg        ! h2 cross-sections (cm2)
    2650 !      real, dimension(nw) :: yieldh2   ! photodissociation yield
     2649      real, dimension(nw) :: yg        ! h2 cross-sections (cm2)
     2650      real, dimension(nw) :: yieldh2   ! photodissociation yield
    26512651
    26522652!     local
    26532653
    2654 !      integer, parameter :: kdata = 1000
    2655 !      real, parameter :: deltax = 1.e-4
    2656 !      real, dimension(kdata) :: x1, y1, x2, y2
    2657 !      real :: xl, xu
    2658 !      integer :: i, iw, n, ierr
    2659 !      integer :: kin, kout ! input/output logical units
    2660 !      character*100 fil
    2661 
    2662 !      kin = 10
     2654      integer, parameter :: kdata = 1000
     2655      real, parameter :: deltax = 1.e-4
     2656      real, dimension(kdata) :: x1, y1, x2, y2
     2657      real :: xl, xu
     2658      integer :: i, iw, n, ierr
     2659      integer :: kin, kout ! input/output logical units
     2660      character*100 fil
     2661
     2662      kin = 10
    26632663
    26642664!     h2 cross sections
    26652665
    2666 !      fil = '/cross_sections/h2secef.txt'
    2667 !      print*, 'section efficace H2: ', fil
    2668 
    2669 !      if(is_master) then
    2670 
    2671 !      OPEN(kin,FILE=fil,STATUS='OLD')
    2672 
    2673 !      n = 792
    2674 !      read(kin,*) ! avoid first line with wavelength = 0.
    2675 !      DO i = 1, n
    2676 !        READ(kin,*) x1(i), y1(i)
    2677 !      ENDDO
    2678 !      CLOSE(kin)
    2679 
    2680 !      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
    2681 !      CALL addpnt(x1,y1,kdata,n,          0.,0.)
    2682 !      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
    2683 !      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
    2684 !      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
     2666      fil = '/cross_sections/h2secef.txt'
     2667      print*, 'section efficace H2: ', fil
     2668
     2669      if(is_master) then
     2670
     2671      OPEN(kin,FILE=fil,STATUS='OLD')
     2672
     2673      n = 792
     2674      read(kin,*) ! avoid first line with wavelength = 0.
     2675      DO i = 1, n
     2676        READ(kin,*) x1(i), y1(i)
     2677      ENDDO
     2678      CLOSE(kin)
     2679
     2680      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
     2681      CALL addpnt(x1,y1,kdata,n,          0.,0.)
     2682      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
     2683      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
     2684      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
    26852685 
    2686 !      IF (ierr .NE. 0) THEN
    2687 !        WRITE(*,*) ierr, fil
    2688 !        STOP
    2689 !      ENDIF
    2690 
    2691 !      endif !is_master
    2692 
    2693 !      call bcast(yg)
     2686      IF (ierr .NE. 0) THEN
     2687        WRITE(*,*) ierr, fil
     2688        STOP
     2689      ENDIF
     2690
     2691      endif !is_master
     2692
     2693      call bcast(yg)
    26942694 
    26952695!     photodissociation yield
    26962696
    2697 !      fil = '/cross_sections/h2_ionef_schunknagy2000.txt'
    2698 
    2699 !      if(is_master) then
    2700 
    2701 !      OPEN(UNIT=kin,FILE=fil,STATUS='old')
    2702 
    2703 !      n = 19
    2704 !      read(kin,*)
    2705 !      DO i = 1, n
    2706 !         READ(kin,*) xl, xu, y2(i)
    2707 !         x2(i) = (xl + xu)/2.
    2708 !         y2(i) = max(1. - y2(i),0.)
    2709 !      END DO
    2710 !      CLOSE (kin)
    2711 
    2712 !      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
    2713 !      CALL addpnt(x2,y2,kdata,n,          0.,0.)
    2714 !      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
    2715 !      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
    2716 !      CALL inter2(nw,wl,yieldh2,n,x2,y2,ierr)
    2717 !      IF (ierr .NE. 0) THEN
    2718 !        WRITE(*,*) ierr, fil
    2719 !        STOP
    2720 !      ENDIF
    2721 
    2722 !      endif !is_master
    2723 
    2724 !      call bcast(yieldh2)
    2725 
    2726 !      end subroutine rdxsh2
     2697      fil = '/cross_sections/h2_ionef_schunknagy2000.txt'
     2698
     2699      if(is_master) then
     2700
     2701      OPEN(UNIT=kin,FILE=fil,STATUS='old')
     2702
     2703      n = 19
     2704      read(kin,*)
     2705      DO i = 1, n
     2706         READ(kin,*) xl, xu, y2(i)
     2707         x2(i) = (xl + xu)/2.
     2708         y2(i) = max(1. - y2(i),0.)
     2709      END DO
     2710      CLOSE (kin)
     2711
     2712      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
     2713      CALL addpnt(x2,y2,kdata,n,          0.,0.)
     2714      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
     2715      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
     2716      CALL inter2(nw,wl,yieldh2,n,x2,y2,ierr)
     2717      IF (ierr .NE. 0) THEN
     2718        WRITE(*,*) ierr, fil
     2719        STOP
     2720      ENDIF
     2721
     2722      endif !is_master
     2723
     2724      call bcast(yieldh2)
     2725
     2726      end subroutine rdxsh2
    27272727
    27282728!==============================================================================
    2729 !
    2730 !      subroutine rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248, yldno2_298)
    2731 !
     2729
     2730      subroutine rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248, yldno2_298)
     2731
    27322732!-----------------------------------------------------------------------------*
    27332733!=  PURPOSE:                                                                 =*
    … …  
    27472747!=           at each defined altitude level                                  =*
    27482748!-----------------------------------------------------------------------------*
    2749 !
    2750 !      use datafile_mod, only: datadir
    2751 !      USE mod_phys_lmdz_para, ONLY: is_master
    2752 !      USE mod_phys_lmdz_transfert_para, ONLY: bcast
    2753 !
    2754 !      implicit none
    2755 !
     2749
     2750      USE mod_phys_lmdz_para, ONLY: is_master
     2751      USE mod_phys_lmdz_transfert_para, ONLY: bcast
     2752
     2753      implicit none
     2754
    27562755!     input
    2757 !
    2758 !      integer :: nw               ! number of wavelength grid points
    2759 !      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
    2760 !     
     2756
     2757      integer :: nw               ! number of wavelength grid points
     2758      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
     2759     
    27612760!     output
    2762 !
    2763 !      real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 cross-sections (cm2)
    2764 !      real, dimension(nw) :: yldno2_248, yldno2_298      ! quantum yields at 248-298 k
    2765 !
     2761
     2762      real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 cross-sections (cm2)
     2763      real, dimension(nw) :: yldno2_248, yldno2_298      ! quantum yields at 248-298 k
     2764
    27662765!     local
    2767 !
    2768 !      integer, parameter :: kdata = 28000
    2769 !      real, parameter :: deltax = 1.e-4
    2770 !      real, dimension(kdata) :: x1, x2, x3, x4, x5, y1, y2, y3, y4, y5
    2771 !      real, dimension(nw) :: yg1, yg2, yg3, yg4, yg5
    2772 !      real :: dum, qy
    2773 !      integer :: i, iw, n, n1, n2, n3, n4, n5, ierr
    2774 !      character*100 fil
    2775 !      integer :: kin, kout ! input/output logical units
    2776 !
    2777 !      kin = 10
    2778 !
     2766
     2767      integer, parameter :: kdata = 28000
     2768      real, parameter :: deltax = 1.e-4
     2769      real, dimension(kdata) :: x1, x2, x3, x4, x5, y1, y2, y3, y4, y5
     2770      real, dimension(nw) :: yg1, yg2, yg3, yg4, yg5
     2771      real :: dum, qy
     2772      integer :: i, iw, n, n1, n2, n3, n4, n5, ierr
     2773      character*100 fil
     2774      integer :: kin, kout ! input/output logical units
     2775
     2776      kin = 10
     2777
    27792778!*************** NO2 photodissociation
    2780 !
     2779
    27812780!  Jenouvrier 1996 + Vandaele 1998 (JPL 2006) 
    2782 !
    2783 !      fil = trim(datadir)//'/cross_sections/no2_xs_jenouvrier.txt'
    2784 !      print*, 'section efficace NO2: ', fil
    2785 !
    2786 !      if(is_master) then
    2787 !
    2788 !      OPEN(UNIT=kin,FILE=fil,status='old')
    2789 !      DO i = 1, 3
    2790 !         READ(kin,*)
    2791 !      ENDDO
    2792 !      n1 = 10001
    2793 !      DO i = 1, n1
    2794 !         READ(kin,*) x1(i), y1(i)
    2795 !      end do
    2796 !
    2797 !      CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax), 0.)
    2798 !      CALL addpnt(x1,y1,kdata,n1,               0., 0.)
    2799 !      CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
    2800 !      CALL addpnt(x1,y1,kdata,n1,           1.e+38, 0.)
    2801 !      CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
    2802 !
    2803 !      endif !is_master
    2804 !
    2805 !      call bcast(yg1)
    2806 !
    2807 !      fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_294K.txt'
    2808 !      print*, 'section efficace NO2: ', fil
    2809 !
    2810 !      if(is_master) then
    2811 !
    2812 !      OPEN(UNIT=kin,FILE=fil,status='old')
    2813 !      DO i = 1, 3
    2814 !         READ(kin,*)
    2815 !      ENDDO
    2816 !      n2 = 27993
    2817 !      DO i = 1, n2
    2818 !         READ(kin,*) x2(i), y2(i)
    2819 !      end do
    2820 !
    2821 !      CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax), 0.)
    2822 !      CALL addpnt(x2,y2,kdata,n2,               0., 0.)
    2823 !      CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
    2824 !      CALL addpnt(x2,y2,kdata,n2,           1.e+38, 0.)
    2825 !      CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
    2826 !
    2827 !      endif !is_master
    2828 !
    2829 !      call bcast(yg2)
    2830 !
    2831 !      fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_220K.txt'
    2832 !      print*, 'section efficace NO2: ', fil
    2833 !
    2834 !      if(is_master) then
    2835 !
    2836 !      OPEN(UNIT=kin,FILE=fil,status='old')
    2837 !      DO i = 1, 3
    2838 !         READ(kin,*)
    2839 !      ENDDO
    2840 !      n3 = 27993
    2841 !      DO i = 1, n3
    2842  !        READ(kin,*) x3(i), y3(i)
    2843   !    end do
    2844 !
    2845  !     CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax), 0.)
    2846   !    CALL addpnt(x3,y3,kdata,n3,               0., 0.)
    2847    !   CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
    2848     !  CALL addpnt(x3,y3,kdata,n3,           1.e+38, 0.)
    2849      ! CALL inter2(nw,wl,yg3,n3,x3,y3,ierr)
    2850 !
    2851 !      do iw = 1, nw - 1
    2852 !         xsno2(iw)     = yg1(iw)
    2853 !         xsno2_294(iw) = yg2(iw)
    2854 !         xsno2_220(iw) = yg3(iw)
    2855 !      end do
    2856 !
    2857 !      endif !is_master
    2858 !
    2859 !      call bcast(yg3)
    2860 !      call bcast(xsno2)
    2861 !      call bcast(xsno2_294)
    2862 !      call bcast(xsno2_220)
    2863 !
     2781
     2782      fil = 'cross_sections/no2_xs_jenouvrier.txt'
     2783      print*, 'section efficace NO2: ', fil
     2784
     2785      if (is_master) then
     2786
     2787      OPEN(UNIT=kin,FILE=fil,status='old')
     2788      DO i = 1, 3
     2789         READ(kin,*)
     2790      END DO
     2791      n1 = 10001
     2792      DO i = 1, n1
     2793         READ(kin,*) x1(i), y1(i)
     2794      end do
     2795
     2796      CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax), 0.)
     2797      CALL addpnt(x1,y1,kdata,n1,               0., 0.)
     2798      CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
     2799      CALL addpnt(x1,y1,kdata,n1,           1.e+38, 0.)
     2800      CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
     2801
     2802      end if !is_master
     2803
     2804      call bcast(yg1)
     2805
     2806      fil = 'cross_sections/no2_xs_vandaele_294K.txt'
     2807      print*, 'section efficace NO2: ', fil
     2808
     2809      if (is_master) then
     2810
     2811      OPEN(UNIT=kin,FILE=fil,status='old')
     2812      DO i = 1, 3
     2813        READ(kin,*)
     2814      END DO
     2815      n2 = 27993
     2816      DO i = 1, n2
     2817         READ(kin,*) x2(i), y2(i)
     2818      end do
     2819
     2820      CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax), 0.)
     2821      CALL addpnt(x2,y2,kdata,n2,               0., 0.)
     2822      CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
     2823      CALL addpnt(x2,y2,kdata,n2,           1.e+38, 0.)
     2824      CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
     2825
     2826      end if !is_master
     2827
     2828      call bcast(yg2)
     2829
     2830      fil = 'cross_sections/no2_xs_vandaele_220K.txt'
     2831      print*, 'section efficace NO2: ', fil
     2832
     2833      if (is_master) then
     2834
     2835      OPEN(UNIT=kin,FILE=fil,status='old')
     2836      DO i = 1, 3
     2837         READ(kin,*)
     2838      END DO
     2839      n3 = 27993
     2840      do i = 1, n3
     2841         READ(kin,*) x3(i), y3(i)
     2842      end do
     2843
     2844      CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax), 0.)
     2845      CALL addpnt(x3,y3,kdata,n3,               0., 0.)
     2846      CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
     2847      CALL addpnt(x3,y3,kdata,n3,           1.e+38, 0.)
     2848      CALL inter2(nw,wl,yg3,n3,x3,y3,ierr)
     2849
     2850      do iw = 1, nw - 1
     2851         xsno2(iw)     = yg1(iw)
     2852         xsno2_294(iw) = yg2(iw)
     2853         xsno2_220(iw) = yg3(iw)
     2854      end do
     2855
     2856      end if !is_master
     2857
     2858      call bcast(yg3)
     2859      call bcast(xsno2)
     2860      call bcast(xsno2_294)
     2861      call bcast(xsno2_220)
     2862
    28642863!     photodissociation efficiency from jpl 2006
    2865 !
    2866 !      fil = trim(datadir)//'/cross_sections/no2_yield_jpl2006.txt'
    2867 !      print*, 'quantum yield NO2: ', fil
    2868 !
    2869 !      if(is_master) then
    2870 !
    2871 !      OPEN(UNIT=kin,FILE=fil,STATUS='old')
    2872 !      DO i = 1, 5
    2873 !         READ(kin,*)
    2874 !      ENDDO
    2875 !      n = 25
    2876 !      n4 = n
    2877 !      n5 = n
    2878 !      DO i = 1, n
    2879 !         READ(kin,*) x4(i), y4(i), y5(i)
    2880 !         x5(i) = x4(i)
    2881 !      ENDDO
    2882 !      CLOSE(kin)
    2883 !
    2884 !      CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),y4(1))
    2885 !      CALL addpnt(x4,y4,kdata,n4,               0.,y4(1))
    2886 !      CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),  0.)
    2887 !      CALL addpnt(x4,y4,kdata,n4,           1.e+38,   0.)
    2888 !      CALL inter2(nw,wl,yg4,n4,x4,y4,ierr)
    2889 !      IF (ierr .NE. 0) THEN
    2890 !         WRITE(*,*) ierr, fil
    2891 !         STOP
    2892 !      ENDIF
    2893 !
    2894 !      endif !is_master
    2895 !
    2896 !      call bcast(yg4)
    2897 !
    2898 !      if(is_master) then
    2899 !
    2900 !      CALL addpnt(x5,y5,kdata,n5,x5(1)*(1.-deltax),y5(1))
    2901 !      CALL addpnt(x5,y5,kdata,n5,               0.,y5(1))
    2902 !      CALL addpnt(x5,y5,kdata,n5,x5(n5)*(1.+deltax),  0.)
    2903 !      CALL addpnt(x5,y5,kdata,n5,           1.e+38,   0.)
    2904 !      CALL inter2(nw,wl,yg5,n5,x5,y5,ierr)
    2905 !      IF (ierr .NE. 0) THEN
    2906 !         WRITE(*,*) ierr, fil
    2907 !         STOP
    2908 !      ENDIF
    2909 !
    2910 !      do iw = 1, nw - 1
    2911 !         yldno2_298(iw) = yg4(iw)
    2912 !         yldno2_248(iw) = yg5(iw)
    2913 !      end do
    2914 !
    2915 !      endif !is_master
    2916 !
    2917 !      call bcast(yg5)
    2918 !      call bcast(yldno2_298)
    2919 !      call bcast(yldno2_248)
    2920 !     
    2921 !      end subroutine rdxsno2
    2922 !
     2864
     2865      fil = 'cross_sections/no2_yield_jpl2006.txt'
     2866      print*, 'quantum yield NO2: ', fil
     2867
     2868      if (is_master) then
     2869
     2870      OPEN(UNIT=kin,FILE=fil,STATUS='old')
     2871      DO i = 1, 5
     2872         READ(kin,*)
     2873      END DO
     2874      n = 25
     2875      n4 = n
     2876      n5 = n
     2877      DO i = 1, n
     2878         READ(kin,*) x4(i), y4(i), y5(i)
     2879         x5(i) = x4(i)
     2880      END DO
     2881      CLOSE(kin)
     2882
     2883      CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),y4(1))
     2884      CALL addpnt(x4,y4,kdata,n4,               0.,y4(1))
     2885      CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),  0.)
     2886      CALL addpnt(x4,y4,kdata,n4,           1.e+38,   0.)
     2887      CALL inter2(nw,wl,yg4,n4,x4,y4,ierr)
     2888      IF (ierr .NE. 0) THEN
     2889         WRITE(*,*) ierr, fil
     2890         STOP
     2891      END IF
     2892
     2893      end if !is_master
     2894
     2895      call bcast(yg4)
     2896
     2897      if (is_master) then
     2898
     2899      CALL addpnt(x5,y5,kdata,n5,x5(1)*(1.-deltax),y5(1))
     2900      CALL addpnt(x5,y5,kdata,n5,               0.,y5(1))
     2901      CALL addpnt(x5,y5,kdata,n5,x5(n5)*(1.+deltax),  0.)
     2902      CALL addpnt(x5,y5,kdata,n5,           1.e+38,   0.)
     2903      CALL inter2(nw,wl,yg5,n5,x5,y5,ierr)
     2904      IF (ierr .NE. 0) THEN
     2905         WRITE(*,*) ierr, fil
     2906         STOP
     2907      END IF
     2908
     2909      do iw = 1, nw - 1
     2910         yldno2_298(iw) = yg4(iw)
     2911         yldno2_248(iw) = yg5(iw)
     2912      end do
     2913
     2914      end if !is_master
     2915
     2916      call bcast(yg5)
     2917      call bcast(yldno2_298)
     2918      call bcast(yldno2_248)
     2919     
     2920      end subroutine rdxsno2
     2921
    29232922!==============================================================================
    29242923
    2925  !     subroutine rdxsno(nw, wl, yg, yieldno)
     2924      subroutine rdxsno(nw, wl, yg, yieldno)
    29262925
    29272926!-----------------------------------------------------------------------------*
    … …  
    29352934!-----------------------------------------------------------------------------*
    29362935
    2937 !      use datafile_mod, only: datadir
    2938 !      USE mod_phys_lmdz_para, ONLY: is_master
    2939 !      USE mod_phys_lmdz_transfert_para, ONLY: bcast
    2940 
    2941 !      implicit none
     2936      USE mod_phys_lmdz_para, ONLY: is_master
     2937      USE mod_phys_lmdz_transfert_para, ONLY: bcast
     2938
     2939      implicit none
    29422940
    29432941!     input
    29442942
    2945 !      integer :: nw               ! number of wavelength grid points
    2946 !      real, dimension(nw) :: wl   ! lower wavelength for each interval
     2943      integer :: nw               ! number of wavelength grid points
     2944      real, dimension(nw) :: wl   ! lower wavelength for each interval
    29472945
    29482946!     output
    29492947
    2950 !      real, dimension(nw) :: yg        ! no cross-sections (cm2)
    2951 !      real, dimension(nw) :: yieldno   ! no photodissociation efficiency
     2948      real, dimension(nw) :: yg        ! no cross-sections (cm2)
     2949      real, dimension(nw) :: yieldno   ! no photodissociation efficiency
    29522950
    29532951!     local
    29542952
    2955 !      integer, parameter :: kdata = 110
    2956 !      real, parameter :: deltax = 1.e-4
    2957 !      real, dimension(kdata) :: x1, y1, x2, y2
    2958 !      integer :: i, iw, n, ierr
    2959 !      character*100 fil
    2960 !      integer :: kin, kout ! input/output logical units
    2961 
    2962 !      kin = 10
     2953      integer, parameter :: kdata = 110
     2954      real, parameter :: deltax = 1.e-4
     2955      real, dimension(kdata) :: x1, y1, x2, y2
     2956      integer :: i, iw, n, ierr
     2957      character*100 fil
     2958      integer :: kin, kout ! input/output logical units
     2959
     2960      kin = 10
    29632961
    29642962!     no cross-sections
    29652963
    2966 !      fil = trim(datadir)//'/cross_sections/no_xs_francisco.txt'
    2967 !      print*, 'section efficace NO: ', fil
    2968 
    2969 !      if(is_master) then
    2970 
    2971 !      OPEN(kin,FILE=fil,STATUS='OLD')
    2972 
    2973 !      n = 99
    2974 !      DO i = 1, n
    2975 !        READ(kin,*) x1(i), y1(i)
    2976 !      ENDDO
    2977 !      CLOSE(kin)
    2978 
    2979 !      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
    2980 !      CALL addpnt(x1,y1,kdata,n,          0.,0.)
    2981 !      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
    2982 !      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
    2983 
    2984 !      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
    2985 !      IF (ierr .NE. 0) THEN
    2986 !         WRITE(*,*) ierr, fil
    2987 !         STOP
    2988 !      ENDIF
    2989 
    2990 !      endif !is_master
    2991 
    2992 !      call bcast(yg)
     2964      fil = 'cross_sections/no_xs_francisco.txt'
     2965      print*, 'section efficace NO: ', fil
     2966
     2967      if (is_master) then
     2968
     2969      OPEN(kin,FILE=fil,STATUS='OLD')
     2970
     2971      n = 99
     2972      DO i = 1, n
     2973        READ(kin,*) x1(i), y1(i)
     2974      END DO
     2975      CLOSE(kin)
     2976
     2977      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
     2978      CALL addpnt(x1,y1,kdata,n,          0.,0.)
     2979      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
     2980      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
     2981
     2982      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
     2983      IF (ierr .NE. 0) THEN
     2984         WRITE(*,*) ierr, fil
     2985         STOP
     2986      END IF
     2987
     2988      end if !is_master
     2989
     2990      call bcast(yg)
    29932991 
    29942992!     photodissociation yield
    29952993
    2996 !      fil = trim(datadir)//'/cross_sections/noefdis.txt'
    2997 
    2998 !      if(is_master) then
    2999 
    3000 !      OPEN(UNIT=kin,FILE=fil,STATUS='old')
    3001 
    3002 !      n = 33
    3003 !      DO i = 1, n
    3004 !         READ(kin,*) x2(n-i+1), y2(n-i+1)
    3005 !      END DO
    3006 !      CLOSE (kin)
    3007 
    3008 !      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
    3009 !      CALL addpnt(x2,y2,kdata,n,          0.,0.)
    3010 !      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
    3011 !      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
    3012 !      CALL inter2(nw,wl,yieldno,n,x2,y2,ierr)
    3013 !      IF (ierr .NE. 0) THEN
    3014 !        WRITE(*,*) ierr, fil
    3015 !        STOP
    3016 !      ENDIF
    3017 
    3018 !      endif !is_master
    3019 
    3020 !      call bcast(yieldno)
    3021 
    3022 !      end subroutine rdxsno
     2994      fil = 'cross_sections/noefdis.txt'
     2995
     2996      if (is_master) then
     2997
     2998      OPEN(UNIT=kin,FILE=fil,STATUS='old')
     2999
     3000      n = 33
     3001      DO i = 1, n
     3002         READ(kin,*) x2(n-i+1), y2(n-i+1)
     3003      END DO
     3004      CLOSE (kin)
     3005
     3006      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
     3007      CALL addpnt(x2,y2,kdata,n,          0.,0.)
     3008      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
     3009      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
     3010      CALL inter2(nw,wl,yieldno,n,x2,y2,ierr)
     3011      IF (ierr .NE. 0) THEN
     3012        WRITE(*,*) ierr, fil
     3013        STOP
     3014      END IF
     3015
     3016      end if !is_master
     3017
     3018      call bcast(yieldno)
     3019
     3020      end subroutine rdxsno
    30233021
    30243022!==============================================================================
    30253023
    3026 !      subroutine rdxsn2(nw, wl, yg, yieldn2)
     3024      subroutine rdxsn2(nw, wl, yg, yieldn2)
    30273025
    30283026!-----------------------------------------------------------------------------*
    … …  
    30353033!-----------------------------------------------------------------------------*
    30363034
    3037 !      use datafile_mod, only: datadir
    3038 !      USE mod_phys_lmdz_para, ONLY: is_master
    3039 !      USE mod_phys_lmdz_transfert_para, ONLY: bcast
    3040 
    3041 !      implicit none
     3035      USE mod_phys_lmdz_para, ONLY: is_master
     3036      USE mod_phys_lmdz_transfert_para, ONLY: bcast
     3037
     3038      implicit none
    30423039
    30433040!     input
    30443041
    3045 !      integer :: nw               ! number of wavelength grid points
    3046 !      real, dimension(nw) :: wl   ! lower wavelength for each interval
     3042      integer :: nw               ! number of wavelength grid points
     3043      real, dimension(nw) :: wl   ! lower wavelength for each interval
    30473044
    30483045!     output
    30493046
    3050 !      real, dimension(nw) :: yg        ! n2 cross-sections (cm2)
    3051 !      real, dimension(nw) :: yieldn2   ! n2 photodissociation yield
     3047      real, dimension(nw) :: yg        ! n2 cross-sections (cm2)
     3048      real, dimension(nw) :: yieldn2   ! n2 photodissociation yield
    30523049
    30533050!     local
    30543051
    3055 !      integer, parameter :: kdata = 1100
    3056 !      real, parameter :: deltax = 1.e-4
    3057 !      real, dimension(kdata) :: x1, y1, x2, y2
    3058 !      real :: xl, xu
    3059 !      integer :: i, iw, n, ierr
    3060 !      integer :: kin, kout ! input/output logical units
    3061 !      character*100 fil
    3062 
    3063 !      kin = 10
     3052      integer, parameter :: kdata = 1100
     3053      real, parameter :: deltax = 1.e-4
     3054      real, dimension(kdata) :: x1, y1, x2, y2
     3055      real :: xl, xu
     3056      integer :: i, iw, n, ierr
     3057      integer :: kin, kout ! input/output logical units
     3058      character*100 fil
     3059
     3060      kin = 10
    30643061
    30653062!     n2 cross sections
    30663063
    3067 !      fil = trim(datadir)//'/cross_sections/n2secef_01nm.txt'
    3068 !      print*, 'section efficace N2: ', fil
    3069 
    3070 !      if(is_master) then
    3071 
    3072 !      OPEN(kin,FILE=fil,STATUS='OLD')
    3073 
    3074 !      n = 1020
    3075 !      DO i = 1, n
    3076 !        READ(kin,*) x1(i), y1(i)
    3077 !      ENDDO
    3078 !      CLOSE(kin)
    3079 
    3080 !      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
    3081 !      CALL addpnt(x1,y1,kdata,n,          0.,0.)
    3082 !      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
    3083 !      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
    3084 
    3085 !      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
     3064      fil = 'cross_sections/n2secef_01nm.txt'
     3065      print*, 'section efficace N2: ', fil
     3066
     3067      if (is_master) then
     3068
     3069      OPEN(kin,FILE=fil,STATUS='OLD')
     3070
     3071      n = 1020
     3072      DO i = 1, n
     3073        READ(kin,*) x1(i), y1(i)
     3074      END DO
     3075      CLOSE(kin)
     3076
     3077      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
     3078      CALL addpnt(x1,y1,kdata,n,          0.,0.)
     3079      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
     3080      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
     3081
     3082      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
    30863083 
    3087 !      IF (ierr .NE. 0) THEN
    3088 !        WRITE(*,*) ierr, fil
    3089 !        STOP
    3090 !      ENDIF
    3091 
    3092 !      endif !is_master
    3093 
    3094 !      call bcast(yg)
     3084      IF (ierr .NE. 0) THEN
     3085        WRITE(*,*) ierr, fil
     3086        STOP
     3087      END IF
     3088
     3089      end if !is_master
     3090
     3091      call bcast(yg)
    30953092 
    30963093!     photodissociation yield
    30973094
    3098 !      fil = trim(datadir)//'/cross_sections/n2_ionef_schunknagy2000.txt'
    3099 
    3100 !      if(is_master) then
    3101 
    3102 !      OPEN(UNIT=kin,FILE=fil,STATUS='old')
    3103 
    3104 !      n = 19
    3105 !      read(kin,*)
    3106 !      DO i = 1, n
    3107 !         READ(kin,*) xl, xu, y2(i)
    3108 !         x2(i) = (xl + xu)/2.
    3109 !         y2(i) = 1. - y2(i)
    3110 !      END DO
    3111 !      CLOSE (kin)
    3112 
    3113 !      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
    3114 !      CALL addpnt(x2,y2,kdata,n,          0.,0.)
    3115 !      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
    3116 !      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
    3117 !      CALL inter2(nw,wl,yieldn2,n,x2,y2,ierr)
    3118 !      IF (ierr .NE. 0) THEN
    3119 !        WRITE(*,*) ierr, fil
    3120 !        STOP
    3121 !      ENDIF
    3122 
    3123 !      endif !is_master
    3124 
    3125 !      call bcast(yieldn2)
    3126 
    3127 !      end subroutine rdxsn2
     3095      fil = 'cross_sections/n2_ionef_schunknagy2000.txt'
     3096
     3097      if (is_master) then
     3098
     3099      OPEN(UNIT=kin,FILE=fil,STATUS='old')
     3100
     3101      n = 19
     3102      read(kin,*)
     3103      DO i = 1, n
     3104         READ(kin,*) xl, xu, y2(i)
     3105         x2(i) = (xl + xu)/2.
     3106         y2(i) = 1. - y2(i)
     3107      END DO
     3108      CLOSE (kin)
     3109
     3110      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
     3111      CALL addpnt(x2,y2,kdata,n,          0.,0.)
     3112      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
     3113      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
     3114      CALL inter2(nw,wl,yieldn2,n,x2,y2,ierr)
     3115      IF (ierr .NE. 0) THEN
     3116        WRITE(*,*) ierr, fil
     3117        STOP
     3118      END IF
     3119
     3120      end if !is_master
     3121
     3122      call bcast(yieldn2)
     3123
     3124      end subroutine rdxsn2
    31283125
    31293126!==============================================================================

  • trunk/LMDZ.VENUS/libf/phyvenus/photolysis_online.F

    r2780 r2795  
    77     $           i_cl2, i_hocl, i_so2, i_so, i_so3,
    88     $           i_clo, i_ocs, i_cocl2, i_h2so4, i_cl,                 
    9      $           nesp, rm,
     9     $           i_no2, i_no, i_n2, i_n2d, nesp, rm,
    1010     $           sza, dist_sol, v_phot)
    1111
    … …  
    1515
    1616!     input
    17 
    18 !     logical, intent(in) :: deutchem
    1917
    2018      integer, intent(in) :: nesp                                    ! total number of chemical species
    … …  
    2422     $                       i_oh, i_ho2, i_h2o2, i_h2o, i_h, i_hcl,
    2523     $                       i_cl2, i_hocl, i_so2, i_so, i_so3, i_clo,
    26      $                       i_ocs, i_cocl2, i_h2so4 , i_cl
     24     $                       i_ocs, i_cl, i_cocl2, i_h2so4, i_no2,
     25     $                       i_no, i_n2, i_n2d
    2726
    2827      real, dimension(nlayer), intent(in) :: press, temp,  zmmean    ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
    … …  
    6766      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o,
    6867     $           j_h2o, j_h2o2, j_ho2, j_h, j_hcl, j_cl2, j_hocl, j_so2,
    69      $           j_so, j_so3, j_clo, j_ocs, j_cocl2, j_h2so4
     68     $           j_so, j_so3, j_clo, j_ocs, j_cocl2, j_h2so4, j_no2,
     69     $           j_no, j_n2
    7070
    7171      integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_hcl, a_cl2,
    7272     $           a_hocl, a_so2, a_so, a_so3, a_clo, a_ocs, a_cocl2,
    73      $           a_h2so4
     73     $           a_h2so4, a_no2, a_no, a_n2
    7474      integer :: i, ilay, iw, ialt
    7575      real    :: deltaj
     76      logical :: deutchem
    7677
    7778!     absorbing gas numbering
    … …  
    9394      a_cocl2   = 15     ! cocl2
    9495      a_h2so4   = 16     ! h2so4
    95 !     a_h2      = 7      ! h2
    96 !     a_no      = 8      ! no
    97 !     a_no2     = 9      ! no2
    98 !     a_n2      = 10     ! n2
     96      a_no2     = 17     ! no2
     97      a_no      = 18     ! no
     98      a_n2      = 19     ! n2
    9999
    100100!     photodissociation rates numbering.
    … …  
    120120      j_cocl2   = 18     ! cocl2 + hv  -> 2cl + co
    121121      j_h2so4   = 19     ! h2so4 + hv  -> so3 + h2o
    122 !     j_h2      = 10     ! h2 + hv     -> h + h
    123 !     j_no      = 11     ! no + hv     -> n + o
    124 !     j_no2     = 12     ! no2 + hv    -> no + o
    125 !     j_n2      = 13     ! n2 + hv     -> n + n
    126 !     j_hdo_od  = 14     ! hdo + hv    -> od + h
    127 !     j_hdo_d   = 15     ! hdo + hv    -> d + oh
     122      j_no2     = 20     ! no2 + hv    -> no + o
     123      j_no      = 21     ! no + hv     -> n + o
     124      j_n2      = 22     ! n2 + hv     -> n(2d) + n
     125
     126!     j_hdo_od  =        ! hdo + hv    -> od + h
     127!     j_hdo_d   =        ! hdo + hv    -> d + oh
    128128
    129129!==== air column increments and rayleigh optical depth
    … …  
    164164!     no2:
    165165
    166 !     call setno2(nphot, nlayer, nw, wc, temp, xsno2, xsno2_220,
    167 !    $            xsno2_294, yldno2_248, yldno2_298, j_no2,
    168 !    $            colinc, rm(:,i_no2), dtgas(:,:,a_no2), sj)
     166      call setno2(nphot, nlayer, nw, wc, temp, xsno2, xsno2_220,
     167     $            xsno2_294, yldno2_248, yldno2_298, j_no2,
     168     $            colinc, rm(:,i_no2), dtgas(:,:,a_no2), sj)
    169169
    170170!==== temperature independent optical depths and cross-sections
    … …  
    188188            dtgas(ilay,iw,a_h2so4) = colinc(ilay)*rm(ilay,i_h2so4)
    189189     $      *xsh2so4(iw)
    190 !           dtgas(ilay,iw,a_h2)  = colinc(ilay)*rm(ilay,i_h2)*xsh2(iw)
    191 !           dtgas(ilay,iw,a_no)  = colinc(ilay)*rm(ilay,i_no)*xsno(iw)
    192 !           dtgas(ilay,iw,a_n2)  = colinc(ilay)*rm(ilay,i_n2)*xsn2(iw)
     190            dtgas(ilay,iw,a_no) = colinc(ilay)*rm(ilay,i_no)*xsno(iw)
     191            dtgas(ilay,iw,a_n2) = colinc(ilay)*rm(ilay,i_n2)*xsn2(iw)
    193192         end do
    194193      end do
    … …  
    221220            sj(ilay,iw,j_cocl2) = xscocl2(iw)            ! cocl2
    222221            sj(ilay,iw,j_h2so4) = xsh2so4(iw)            ! h2so4
    223 !           sj(ilay,iw,j_h2)  = xsh2(iw)*yieldh2(iw)     ! h2
    224 !           sj(ilay,iw,j_no)  = xsno(iw)*yieldno(iw)     ! no
    225 !           sj(ilay,iw,j_n2)  = xsn2(iw)*yieldn2(iw)     ! n2
     222            sj(ilay,iw,j_no)  = xsno(iw)*yieldno(iw)     ! no
     223            sj(ilay,iw,j_n2)  = xsn2(iw)*yieldn2(iw)     ! n2
    226224         end do
    227225      end do
    228226
    229 !      !HDO cross section
     227!      hdo cross section
     228
    230229!      if (deutchem) then
    231230!         do ilay=1,nlayer
    … …  
    238237!      endif
    239238
    240 !==== set solid aerosol properties and optical depth
    241 
    242       tau = 0.  ! no solid aerosols
     239!==== set aerosol properties and optical depth
     240
     241      tau = 0. ! no solid aerosols for the present time
    243242
    244243      call setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
    … …  
    293292         end do
    294293
    295 !        write(45,*) wc(iw), v_phot(nlayer,3), v_phot(nlayer,4)
    296 
    297294!     eliminate small values
    298295
    … …  
    300297            v_phot(:,1:nphot) = 0.
    301298         end where
    302 
     299         
    303300      end do ! iw
    304301
    … …  
    438435
    439436            yieldco2(:) = 1
     437
    440438            if (wc(l) >= 167.) then
    441439               sj(i,l,j_co2_o) = sco2*yieldco2(l)
    … …  
    666664     $                              /(360. - 298.)*(temp - 298.)
    667665            end if
    668 
    669666!     219 nm photolysis treshold
    670 
    671667            if (wc(iw) <= 219.) then
    672668               sj(ilev,iw,j_so2) = xsso2
    … …  
    763759
    764760!==============================================================================
    765 !
    766 !      subroutine setno2(nd, nlayer, nw, wc, tlay, xsno2, xsno2_220,
    767 !     $                  xsno2_294, yldno2_248, yldno2_298, j_no2,
    768 !     $                  colinc, rm, dt, sj)
    769 !
     761
     762      subroutine setno2(nd, nlayer, nw, wc, tlay, xsno2, xsno2_220,
     763     $                  xsno2_294, yldno2_248, yldno2_298, j_no2,
     764     $                  colinc, rm, dt, sj)
     765
    770766!-----------------------------------------------------------------------------*
    771767!=  PURPOSE:                                                                 =*
    772768!=  Set up the no2 temperature-dependent cross-sections and optical depth    =*
    773769!-----------------------------------------------------------------------------*
    774 !
    775 !      implicit none
    776 !
     770
     771      implicit none
     772
    777773!     input:
    778 !
    779 !      integer :: nd                                            ! number of photolysis rates
    780 !      integer :: nlayer                                        ! number of vertical layers
    781 !      integer :: nw                                            ! number of wavelength grid points
    782 !      integer :: j_no2                                         ! photolysis index
    783 !      real, dimension(nw)     :: wc                            ! central wavelength for each interval
    784 !      real, dimension(nw) :: xsno2, xsno2_220, xsno2_294       ! no2 absorption cross-section at 220-294 k (cm2)
    785 !      real, dimension(nw) :: yldno2_248, yldno2_298            ! no2 quantum yield at 248-298 k
    786 !      real, dimension(nlayer) :: tlay                          ! temperature (k)
    787 !      real, dimension(nlayer) :: rm                            ! no2 mixing ratio
    788 !      real, dimension(nlayer) :: colinc                        ! air column increment (molecule.cm-2)
     774
     775      integer :: nd                                            ! number of photolysis rates
     776      integer :: nlayer                                        ! number of vertical layers
     777      integer :: nw                                            ! number of wavelength grid points
     778      integer :: j_no2                                         ! photolysis index
     779      real, dimension(nw)     :: wc                            ! central wavelength for each interval
     780      real, dimension(nw) :: xsno2, xsno2_220, xsno2_294       ! no2 absorption cross-section at 220-294 k (cm2)
     781      real, dimension(nw) :: yldno2_248, yldno2_298            ! no2 quantum yield at 248-298 k
     782      real, dimension(nlayer) :: tlay                          ! temperature (k)
     783      real, dimension(nlayer) :: rm                            ! no2 mixing ratio
     784      real, dimension(nlayer) :: colinc                        ! air column increment (molecule.cm-2)
    789785
    790786!     output:
    791787
    792 !      real, dimension(nlayer,nw)    :: dt                      !  optical depth
    793 !      real, dimension(nlayer,nw,nd) :: sj                      !  cross-section array (cm2)
     788      real, dimension(nlayer,nw)    :: dt                      !  optical depth
     789      real, dimension(nlayer,nw,nd) :: sj                      !  cross-section array (cm2)
    794790
    795791!     local:
    796792
    797 !      integer :: ilev, iw
    798 !      real    :: temp, qy
     793      integer :: ilev, iw
     794      real    :: temp, qy
    799795
    800796!     temperature dependance: jpl 2006
    801797
    802 !      do ilev = 1,nlayer
    803 !         temp = max(220.,min(tlay(ilev),294.))
    804 !         do iw = 1, nw - 1
    805 !            if (wc(iw) < 238.) then
    806 !               sj(ilev,iw,j_no2) = xsno2(iw)
    807 !            else
    808 !               sj(ilev,iw,j_no2) = xsno2_220(iw)
    809 !     $                           + (xsno2_294(iw) - xsno2_220(iw))
    810 !     $                            /(294. - 220.)*(temp - 220.)
    811 !            end if
     798      do ilev = 1,nlayer
     799         temp = max(220.,min(tlay(ilev),294.))
     800         do iw = 1, nw - 1
     801            if (wc(iw) < 238.) then
     802               sj(ilev,iw,j_no2) = xsno2(iw)
     803            else
     804               sj(ilev,iw,j_no2) = xsno2_220(iw)
     805     $                           + (xsno2_294(iw) - xsno2_220(iw))
     806     $                            /(294. - 220.)*(temp - 220.)
     807            end if
    812808
    813809!     optical depth
    814810
    815 !            dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_no2)
    816 !         end do
    817 !      end do
     811            dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_no2)
     812         end do
     813      end do
    818814
    819815!     quantum yield: jpl 2006
    820816
    821 !      do ilev = 1,nlayer
    822 !         temp = max(248.,min(tlay(ilev),298.))
    823 !         do iw = 1, nw - 1
    824 !            qy = yldno2_248(iw) + (yldno2_298(iw) - yldno2_248(iw))
    825 !     $                           /(298. - 248.)*(temp - 248.)
    826 !            sj(ilev,iw,j_no2) = sj(ilev,iw,j_no2)*qy
    827 !         end do
    828 !      end do
    829 !
    830 !     end subroutine setno2
     817      do ilev = 1,nlayer
     818         temp = max(248.,min(tlay(ilev),298.))
     819         do iw = 1, nw - 1
     820            qy = yldno2_248(iw) + (yldno2_298(iw) - yldno2_248(iw))
     821     $                           /(298. - 248.)*(temp - 248.)
     822            sj(ilev,iw,j_no2) = sj(ilev,iw,j_no2)*qy
     823         end do
     824      end do
     825
     826      end subroutine setno2
    831827
    832828!==============================================================================
    … …  
    867863      gamma  = 0.03  ! conrath parameter
    868864
    869       do ilay = 1, nlayer - 1
    870          do iw = 1, nw - 1
    871             dtaer(ilay,iw) = 0.       ! no aer for the moment
    872             omaer(ilay,iw) = 0.622
    873             gaer(ilay,iw)  = 0.88
    874          end do
    875       end do
     865      dtaer(:,:) = 0.
     866      omaer(:,:) = 0.
     867      gaer(:,:)  = 0.
     868     
     869      omaer(:,:) = omega
     870      gaer(:,:) =g
     871
     872!     optical depth profile:
     873
     874       tautot = 0.
     875       do ilay = 1, nlayer-1
     876          dz = alt(ilay+1) - alt(ilay)
     877          tautot = tautot + exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz
     878       end do
     879     
     880       q0 = tau/tautot
     881       do ilay = 1, nlayer-1
     882          dz = alt(ilay+1) - alt(ilay)
     883          dtaer(ilay,:) = q0*exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz
     884          omaer(ilay,:) = omega
     885          gaer(ilay,:)  = g
     886       end do
    876887
    877888      end subroutine setaer
    … …  
    960971            omcld(i,iw) = 0.
    961972            gcld(i,iw) = 0.
    962         end do
     973         end do
    963974      end do
    964 
     975c
    965976c     if you do not want any clouds, set clouds = .false.
    966 
     977c
    967978      clouds = .true.
    968 
     979c
    969980      if (.not. clouds) then
    970981         write(kout,*) 'no clouds'
    971982         return
    972983      end if
    973 
     984c
    974985* cloud properties are set for each layer (not each level)
    975986
    … …  
    10151026* for g and omega, use averages weigthed by optical depth
    10161027
    1017 !      DO 11, i = 1, n    !***** CHANGED!!See header!!*****
     1028!     DO 11, i = 1, n    !***** CHANGED!!See header!!*****
    10181029      DO 11, i = 1, n-1
    10191030         womd(i) = omd(i) * cd(i)
    … …  
    10241035      CALL inter3(nz,z,gz , n, zd,wgd,0)
    10251036
     1037
    10261038! Imprimer Cz et imprimer cd pour comparer
     1039
    10271040
    10281041      DO 15, i = 1, nz-1
    … …  
    10471060   17 CONTINUE
    10481061
     1062*_______________________________________________________________________
     1063
    10491064      RETURN
    1050       END
     1065      END                                                                                                                                                                                                       
    10511066
    10521067!==============================================================================
    … …  
    14121427! CUPTN and CDNTN = calc. when TAU is TAUN
    14131428! DIVISR = prevents division by zero
    1414 
    14151429      do j = 0, nlev
    14161430         tauc(j) = 0.
    … …  
    17161730      fdn(lev) = edn(lev)/mu1(lev)
    17171731      fup(lev) = eup(lev)/mu1(lev)
    1718 
    17191732      DO 60, lev = 2, nlayer + 1
    17201733         fdr(lev) = EXP(-tausla(lev-1))

  • trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F

    r2788 r2795  
    456456         allocate(prod_tr(klon,klev,nqmax))
    457457         allocate(loss_tr(klon,klev,nqmax))
     458         allocate(no_emission(klon,klev))
     459         allocate(o2_emission(klon,klev))
    458460
    459461#ifdef CPP_XIOS
    … …  
    10701072     $                       iter,
    10711073     $                       prod_tr,
    1072      $                       loss_tr)
     1074     $                       loss_tr,
     1075     $                       no_emission,
     1076     $                       o2_emission)
    10731077
    10741078         if (ok_sedim) then
    … …  
    21232127
    21242128! tracers in gas phase, column densities
     2129
    21252130         do iq = 1,nqmax - nmicro
    21262131            col_dens_tr(:,iq)=0.
    … …  
    21312136            call send_xios_field("col_"//tname(iq),col_dens_tr(:,iq))
    21322137         end do
    2133 
    21342138
    21352139! tracers in liquid phase, volume mixing ratio
    … …  
    21442148            end if
    21452149         end if
     2150
     2151! aeronomical emissions
     2152
     2153!        call send_xios_field("no_emission",no_emission)
     2154!        call send_xios_field("o2_emission",o2_emission)
    21462155
    21472156! chemical iterations

  • trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F

    r2780 r2795  
    1515     $                    iter,
    1616     $                    prod_tr,
    17      $                    loss_tr)
     17     $                    loss_tr,
     18     $                    no_emi,
     19     $                    o2_emi)
    1820
    1921      use chemparam_mod
    … …  
    4850!===================================================================
    4951
    50       real, dimension(nlon,nlev,nqmax) :: d_tr_chem  ! chemical tendency for each tracer
    51       integer, dimension(nlon,nlev) :: iter          ! chemical iterations
    52       real, dimension(nlon,nlev,nqmax) :: prod_tr, loss_tr  ! produc and loss tracer
     52      integer, dimension(nlon,nlev) :: iter                 ! chemical iterations
     53      real, dimension(nlon,nlev,nqmax) :: d_tr_chem         ! chemical tendency for each tracer
     54      real, dimension(nlon,nlev,nqmax) :: prod_tr, loss_tr  ! production and loss terms (for info)
     55      real, dimension(nlon,nlev)    :: no_emi               ! no emission
     56      real, dimension(nlon,nlev)    :: o2_emi               ! o2 emission
    5357
    5458!===================================================================
    … …  
    7074      real, dimension(nlon,nlev) :: trac_sum
    7175      real, dimension(nlon,nlev,nqmax) :: ztrac  ! local tracer mixing ratio
     76      real, dimension(nlev) :: no_em
     77      real, dimension(nlev) :: o2_em
    7278     
    7379!===================================================================
    … …  
    299305     $                                nqmax, iter(ilon,:),
    300306     $                                prod_tr(ilon,:,:),
    301      $                                loss_tr(ilon,:,:))
     307     $                                loss_tr(ilon,:,:),
     308     $                                no_em, o2_em)
     309
     310             no_emi(ilon,:) = no_em(:)
     311             o2_emi(ilon,:) = o2_em(:)
    302312
    303313         end do         
Note: See TracChangeset for help on using the changeset viewer.