Changeset 2836


Ignore:
Timestamp:
Nov 30, 2022, 4:46:57 PM (2 years ago)
Author:
abierjon
Message:

VENUS GCM:

INCLUDING THE IONOSPHERE CODE IN VENUS GCM

ATTENTION: INCREASED MEMORY DEMAND

NEEDS AT LEAST 135 GB of allocated memory

==============================================================
===> LIST OF APPENDED FILES AND INTERNAL ADDITIONS <==========
==============================================================

NEW MODEL SPECIES

  • deftank/traceur-chemistry-IONOSPHERE.def

Neutrals: 4 Species: N, NO, NO2, N(2D)

Ions: 15 species:

CO2+, CO+, O+, O2+, H+,

N2+, H2O+, OH+, C+, HCO+,

H3O+, HCO2+, N+, NO+, elec

FROM 36 to 55 chemical species

NEW KEYWORD OF PHYSIQ.DEF

  • deftank/physiq-96x96x90-chemistry-IONOSPHERE.def
  • ok_ionchem: keyword supposed to activate ion chemistry.

(be careful that n, no, no2, n2d and the ion species are in the deftank/traceur-chemistry-IONOSPHERE.def)

  • ok_jonline: keyword supposed to activate the online photochemistry

(be careful that n, no, no2 and n2d are in the traceur-chemistry-IONOSPHERE.def)

==============================================================
===> LIST OF MODIFIED PROGRAMS <=========================
==============================================================

nonoro_gwd_ran_mod.F90

  • Change EPFLUXMAX value from 5.E-3 to 1.E-3

photochemistry_venus.F90 (krates; photolysis_online; indices;)

  • Import of keywords: ok_ionchem, tuneupperatm & ok_jonline
  • addition of ion species in order
  • Forcing electroneutrality
  • Update of the reactions a001 and a002 with taking into account the other species
  • Change of formula for a002 with the formulation of Baulch et al., 1976 (confirmed by Smith and Robertson, 2008)

photolysis_mod.F90

  • Modification of the rdsolarflux subroutine to include interpolation with Atlas1 and Atlas3 in connection with E10.7

concentration2.F90

  • Added chemical species n2d, no, no2 and n in the conductivity calculation.

The ions have been excluded because their sum is 105 times less dense than the neutrals and
their thermal conductivity is unknown

iono_h.F90

  • Addition of the phdisrate routine
  • replace the electronic temperature of Mars by that of

origin = 1: Theis et al. 1980 (Venus) with bilinear interpolation altitude/cos(SZA)
origin = 2: Theis et al. 1984 (Venus) with the formula of the electronic temperature at high solar activity

  • addition of an ion temperature model based on VIRA

cleshphys.h

  • added ok_ionchem & ok_jonline in COMMON/clesphys_l/

conf_phys.f90

  • add ok_ionchem & ok_jonline as parameters to read from physiq.def file set to .false. by default

chemparam_mod.F90

  • Add species in M_tr and corresponding i_X. Set all i_X to zero before reading traceur-chemistry-IONOSPHERE.def
  • Added Type_tr table to differentiate species: 1 == neutral, 2 == ION, 3 == liquid, 10 == others


euvheat.F90; hrtherm.F; jthermcalc_e107.F; param_read_e107.F

  • Normalization with Mars

A.M

Location:
trunk/LMDZ.VENUS
Files:
3 added
15 edited

Legend:

Unmodified
Added
Removed

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

    r2795 r2836  
    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, i_no2,     &
    16                  i_no, i_n, i_n2d
     15                 i_clso2, i_oscl, i_n2, i_he, i_n, i_no, &
     16                 i_no2, i_n2d,                           &
     17                 i_co2plus, i_coplus, i_oplus, i_o2plus, &
     18                 i_n2plus, i_hplus, i_h2oplus, i_nplus,  &
     19                 i_ohplus, i_cplus, i_noplus, i_h3oplus, &
     20                 i_hcoplus, i_hco2plus, i_elec     
    1721                 
    1822INTEGER, SAVE :: i_h2oliq, i_h2so4liq
    … …  
    2630INTEGER, SAVE :: nmicro  ! number of microphysical tracers
    2731
    28 REAL, DIMENSION(:), SAVE, ALLOCATABLE :: M_tr
     32REAL, DIMENSION(:), SAVE, ALLOCATABLE :: M_tr     ! Molecular Mass of tracers
     33REAL, DIMENSION(:), SAVE, ALLOCATABLE :: type_tr  ! type of tracer
    2934 
    3035REAL, DIMENSION(:,:),SAVE, ALLOCATABLE :: no_emission
    … …  
    697702  INTEGER :: i
    698703 
    699   ALLOCATE(M_tr(nqtot))
    700  
     704  ALLOCATE(M_tr(nqtot))    ! Molecular Mass of Chemistry Tracers
     705  ALLOCATE(type_tr(nqtot)) ! Type of chemistry tracers 1: neutral, 2: ion, 3: liquide
     706
     707! Initialisation index chemistry tracers
     708
     709  ! Neutral gas Tracer
     710  i_co2      = 0
     711  i_co       = 0
     712  i_h2       = 0
     713  i_h2o      = 0
     714  i_o1d      = 0
     715  i_o        = 0
     716  i_o2       = 0
     717  i_o2dg     = 0
     718  i_o3       = 0
     719  i_h        = 0
     720  i_oh       = 0
     721  i_ho2      = 0
     722  i_h2o2     = 0
     723  i_cl       = 0
     724  i_clo      = 0
     725  i_cl2      = 0
     726  i_hcl      = 0
     727  i_hocl     = 0
     728  i_clco     = 0
     729  i_clco3    = 0
     730  i_cocl2    = 0
     731  i_s        = 0
     732  i_so       = 0
     733  i_so2      = 0
     734  i_so3      = 0
     735  i_s2o2     = 0
     736  i_ocs      = 0
     737  i_hso3     = 0
     738  i_h2so4    = 0
     739  i_s2       = 0
     740  i_clso2    = 0
     741  i_oscl     = 0
     742  i_n2       = 0
     743  i_he       = 0
     744  i_n2d      = 0
     745  i_n        = 0
     746  i_no       = 0
     747  i_no2      = 0
     748  ! ION TRACERS FOR IONCHEM = TRUE
     749  i_co2plus  = 0
     750  i_coplus   = 0
     751  i_oplus    = 0
     752  i_o2plus   = 0
     753  i_n2plus   = 0
     754  i_hplus    = 0
     755  i_h2oplus  = 0
     756  i_nplus    = 0
     757  i_ohplus   = 0
     758  i_cplus    = 0
     759  i_noplus   = 0
     760  i_h3oplus  = 0
     761  i_hcoplus  = 0
     762  i_hco2plus = 0
     763  i_elec     = 0 
     764  ! MICROPHYSICAL TRACERS FOR CL_SCHEME=1
     765  i_h2oliq   = 0
     766  i_h2so4liq = 0 
     767
    701768  DO i=1, nqtot                                 
    702769       
    … …  
    705772       
    706773        SELECT CASE(tname(i))
     774! NEUTRAL TRACERS
    707775                CASE('co2')
    708776                i_co2=i
    709777                PRINT*,'co2',i_co2
    710                 M_tr(i_co2) = 44.0095 
     778                M_tr(i_co2) = 44.0095
     779      type_tr(i_co2) = 1 
    711780                CASE('co')
    712781                i_co=i
    713782                PRINT*,'co',i_co
    714                 M_tr(i_co)=28.0101
     783                M_tr(i_co)=28.0101
     784      type_tr(i_co) = 1
    715785                CASE('h2')
    716786                i_h2=i
    717787                PRINT*,'h2',i_h2
    718788                M_tr(i_h2)= 2.01588
     789      type_tr(i_h2) = 1
    719790                CASE('h2o')
    720791                i_h2o=i
    721792                PRINT*,'h2o',i_h2o
    722793                M_tr(i_h2o)=18.0153
     794      type_tr(i_h2o) = 1
    723795                CASE('o1d')
    724796                i_o1d=i
    725797                PRINT*,'o1d',i_o1d
    726                 M_tr(i_o1d)=15.994
     798                M_tr(i_o1d)=15.994
     799      type_tr(i_o1d) = 1
    727800                CASE('o')
    728801                i_o=i
    729802                PRINT*,'o',i_o
    730                 M_tr(i_o)=15.994
     803                M_tr(i_o)=15.994
     804      type_tr(i_o) = 1
    731805                CASE('o2')
    732806                i_o2=i
    733807                PRINT*,'o2',i_o2
    734808                M_tr(i_o2)=31.9988
     809      type_tr(i_o2) = 1
    735810                CASE('o2dg')
    736811                i_o2dg=i
    737812                PRINT*,'o2dg',i_o2dg
    738813                M_tr(i_o2dg)=31.9988
     814      type_tr(i_o2dg) = 1
    739815                CASE('o3')
    740816                i_o3=i
    741817                PRINT*,'o3',i_o3
    742                 M_tr(i_o3)= 47.9982
     818                M_tr(i_o3)= 47.9982
     819      type_tr(i_o3) = 1
    743820                CASE('h')
    744821                i_h=i
    745822                PRINT*,'h',i_h
    746                 M_tr(i_h)= 1.00794
     823                M_tr(i_h)= 1.00794
     824      type_tr(i_h) = 1 
    747825                CASE('oh')
    748826                i_oh=i
    749827                PRINT*,'oh',i_oh
    750828                M_tr(i_oh)=17.0073
     829      type_tr(i_oh) = 1
    751830                CASE('ho2')
    752831                i_ho2=i
    753832                PRINT*,'ho2',i_ho2
    754833                M_tr(i_ho2)=33.0067
     834      type_tr(i_ho2) = 1
    755835                CASE('h2o2')
    756836                i_h2o2=i
    757837                PRINT*,'h2o2',i_h2o2
    758                 M_tr(i_h2o2)=34.0147
     838                M_tr(i_h2o2)=34.0147
     839      type_tr(i_h2o2) = 1
    759840                CASE('cl')
    760841                i_cl=i
    761842                PRINT*,'cl',i_cl
    762843                M_tr(i_cl)=35.453
     844      type_tr(i_cl) = 1
    763845                CASE('clo')
    764846                i_clo=i
    765847                PRINT*,'clo',i_clo
    766848                M_tr(i_clo)=51.452
     849      type_tr(i_clo) = 1
    767850                CASE('cl2')
    768851                i_cl2=i
    769852                PRINT*,'cl2',i_cl2
    770853                M_tr(i_cl2)=70.906
     854      type_tr(i_cl2) = 1
    771855                CASE('hcl')
    772856                i_hcl=i
    773857                PRINT*,'hcl',i_hcl
    774858                M_tr(i_hcl)=36.461
     859      type_tr(i_hcl) = 1
    775860                CASE('hocl')
    776861                i_hocl=i
    777862                PRINT*,'hocl',i_hocl
    778863                M_tr(i_hocl)=52.46
    779                 CASE('clco')
     864                type_tr(i_hocl) = 1
     865      CASE('clco')
    780866                i_clco=i
    781867                PRINT*,'clco',i_clco
    782868                M_tr(i_clco)=63.463
    783                 CASE('clco3')
     869                type_tr(i_clco) = 1
     870      CASE('clco3')
    784871                i_clco3=i
    785872                PRINT*,'clco3',i_clco3
    786873                M_tr(i_clco3)=95.462
    787                 CASE('cocl2')
     874                type_tr(i_clco3) = 1
     875      CASE('cocl2')
    788876                i_cocl2=i
    789877                PRINT*,'cocl2',i_cocl2
    790878                M_tr(i_cocl2)=98.916
    791                 CASE('s')
     879                type_tr(i_cocl2) = 1
     880      CASE('s')
    792881                i_s=i
    793882                PRINT*,'s',i_s
    794883                M_tr(i_s)=32.065
    795                 CASE('so')
     884                type_tr(i_s) = 1
     885      CASE('so')
    796886                i_so=i
    797887                PRINT*,'so',i_so
    798888                M_tr(i_so)=48.0644
    799                 CASE('so2')
     889                type_tr(i_so) = 1
     890      CASE('so2')
    800891                i_so2=i
    801892                PRINT*,'so2',i_so2
    802893                M_tr(i_so2)=64.064
    803                 CASE('so3')
     894                type_tr(i_so2) = 1
     895      CASE('so3')
    804896                i_so3=i
    805897                PRINT*,'so3',i_so3
    806898                M_tr(i_so3)=80.063
    807                 CASE('s2o2')
     899                type_tr(i_so3) = 1
     900      CASE('s2o2')
    808901                i_s2o2=i
    809902                PRINT*,'s2o2',i_s2o2
    810903                M_tr(i_s2o2)= 96.1288
    811                 CASE('ocs')
     904                type_tr(i_s2o2) = 1
     905      CASE('ocs')
    812906                i_ocs=i
    813907                PRINT*,'ocs',i_ocs
    814908                M_tr(i_ocs)=60.0751
    815                 CASE('hso3')
     909                type_tr(i_ocs) = 1
     910      CASE('hso3')
    816911                i_hso3=i
    817912                PRINT*,'hso3',i_hso3
    818913                M_tr(i_hso3)=81.071
    819                 CASE('h2so4')
     914                type_tr(i_hso3) = 1
     915      CASE('h2so4')
    820916                i_h2so4=i
    821917                PRINT*,'h2so4',i_h2so4
    822918                M_tr(i_h2so4)=98.078
    823                 CASE('s2')
     919                type_tr(i_h2so4) = 1
     920      CASE('s2')
    824921                i_s2=i
    825922                PRINT*,'s2',i_s2
    826923                M_tr(i_s2)=64.13
    827                 CASE('clso2')
     924                type_tr(i_s2) = 1
     925      CASE('clso2')
    828926                i_clso2=i
    829927                PRINT*,'clso2',i_clso2
    830928                M_tr(i_clso2)=99.517
    831                 CASE('oscl')
     929                type_tr(i_clso2) = 1
     930      CASE('oscl')
    832931                i_oscl=i
    833932                PRINT*,'oscl',i_oscl
    834933                M_tr(i_oscl)=83.517
    835                 CASE('n2')
    836                 i_n2=i
    837                 M_tr(i_n2)=28.013
    838                 CASE('he')
    839                 i_he=i
    840                 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
     934                  type_tr(i_oscl) = 1
     935      CASE('n2')
     936                  i_n2=i
     937      PRINT*,'n2',i_n2
     938                  M_tr(i_n2)=28.013
     939                  type_tr(i_n2) = 1
     940      CASE('he')
     941                  i_he=i
     942      PRINT*,'he',i_he
     943                  M_tr(i_he)=4.0026
     944                  type_tr(i_he) = 1
     945      CASE('n2d')
     946                  i_n2d=i
     947      PRINT*,'n2d',i_n2d
     948                  M_tr(i_n2d)=14.0067
     949                  type_tr(i_n2d) = 1
     950      CASE('n')
     951                  i_n=i
     952      PRINT*,'n',i_n
     953                  M_tr(i_n)=14.0067
     954                  type_tr(i_n) = 1
     955      CASE('no')
     956                  i_no=i
     957      PRINT*,'no',i_no
     958                  M_tr(i_no)=30.0061   
     959                  type_tr(i_no) = 1
     960      CASE('no2')
     961                  i_no2=i
     962      PRINT*,'no2',i_no2
     963                  M_tr(i_no2)=46.0055
     964      type_tr(i_no2) = 1         
     965! ION TRACERS FOR OK_IONCHEM = TRUE
     966      CASE('co2plus')
     967      i_co2plus=i
     968      PRINT*,'co2plus',i_co2plus
     969      M_tr(i_co2plus)=44.0095
     970      type_tr(i_co2plus) = 2
     971      CASE('coplus')
     972      i_coplus=i
     973      PRINT*,'coplus',i_coplus
     974      M_tr(i_coplus)=28.0101
     975      type_tr(i_coplus) = 2
     976      CASE('oplus')
     977      i_oplus=i
     978      PRINT*,'oplus',i_oplus
     979      M_tr(i_oplus)=15.994
     980      type_tr(i_oplus) = 2
     981      CASE('o2plus')
     982      i_o2plus=i
     983      PRINT*,'o2plus',i_o2plus
     984      M_tr(i_o2plus)=31.9988
     985      type_tr(i_o2plus) = 2
     986      CASE('n2plus')
     987      i_n2plus=i
     988      PRINT*,'n2plus',i_n2plus
     989      M_tr(i_n2plus)=28.013
     990      type_tr(i_n2plus) = 2
     991      CASE('hplus')
     992      i_hplus=i
     993      PRINT*,'hplus',i_hplus
     994      M_tr(i_hplus)=1.00794
     995                  type_tr(i_hplus) = 2
     996      CASE('h2oplus')
     997      i_h2oplus=i
     998      PRINT*,'h2oplus',i_h2oplus
     999      M_tr(i_h2oplus)=18.0153
     1000                  type_tr(i_h2oplus) = 2     
     1001      CASE('nplus')
     1002                  i_nplus=i
     1003      PRINT*,'nplus',i_nplus
     1004                  M_tr(i_nplus)=14.0067
     1005                  type_tr(i_nplus) = 2
     1006      CASE('ohplus')
     1007                  i_ohplus=i
     1008      PRINT*,'ohplus',i_ohplus
     1009                  M_tr(i_ohplus)=17.0073
     1010                  type_tr(i_ohplus) = 2
     1011      CASE('cplus')
     1012                  i_cplus=i
     1013      PRINT*,'cplus',i_cplus
     1014                  M_tr(i_cplus)=12.011
     1015                  type_tr(i_cplus) = 2
     1016      CASE('noplus')
     1017                  i_noplus=i
     1018      PRINT*,'noplus',i_noplus
     1019                  M_tr(i_noplus)=30.0061
     1020                  type_tr(i_noplus) = 2
     1021      CASE('h3oplus')
     1022                  i_h3oplus=i
     1023      PRINT*,'h3oplus',i_h3oplus
     1024                  M_tr(i_h3oplus)=19.0232
     1025                  type_tr(i_h3oplus) = 2
     1026      CASE('hcoplus')
     1027                  i_hcoplus=i
     1028      PRINT*,'hcoplus',i_hcoplus
     1029                  M_tr(i_hcoplus)=29.0180
     1030                  type_tr(i_hcoplus) = 2
     1031      CASE('hco2plus')
     1032                  i_hco2plus=i
     1033      PRINT*,'hco2plus',i_hco2plus
     1034                  M_tr(i_hco2plus)=45.
     1035                  type_tr(i_hco2plus) = 2
     1036      CASE('elec')
     1037                  i_elec=i
     1038      PRINT*,'elec',i_elec
     1039                  M_tr(i_elec)=1./1822.89
     1040      type_tr(i_elec) = 2       
    8531041! MICROPHYSICAL TRACERS FOR CL_SCHEME=1
    8541042                CASE('h2oliq')
    … …  
    8561044                PRINT*,'h2oliq',i_h2oliq
    8571045                M_tr(i_h2oliq)=18.0153
     1046      type_tr(i_h2oliq) = 3
    8581047                CASE('h2so4liq')
    8591048                i_h2so4liq=i
    8601049                PRINT*,'h2so4liq',i_h2so4liq
    8611050                M_tr(i_h2so4liq)=98.078
     1051      type_tr(i_h2so4liq) = 3
    8621052! MICROPHYSICAL TRACERS FOR CL_SCHEME=2
    863                 CASE('M0_aer')
     1053    CASE('M0_aer')
    8641054                i_m0_aer=i
     1055    type_tr(i_m0_aer) = 10
    8651056                PRINT*,'M0_aer',i_m0_aer
    8661057                CASE('M3_aer')
    8671058                i_m3_aer=i
     1059    type_tr(i_m3_aer) = 10
    8681060                PRINT*,'M3_aer',i_m3_aer
    869                 CASE('M0_m1drop')
     1061    CASE('M0_m1drop')
    8701062                i_m0_mode1drop=i
     1063    type_tr(i_m0_mode1drop) = 10
    8711064                PRINT*,'M0_m1drop',i_m0_mode1drop
    872                 CASE('M0_m1ccn')
     1065    CASE('M0_m1ccn')
    8731066                i_m0_mode1ccn=i
     1067    type_tr(i_m0_mode1ccn) = 10
    8741068                PRINT*,'M0_m1ccn',i_m0_mode1ccn
    8751069                CASE('M3_m1sa')
    8761070                i_m3_mode1sa=i
     1071    type_tr(i_m3_mode1sa) = 10
    8771072                PRINT*,'M3_m1sa',i_m3_mode1sa
    8781073                CASE('M3_m1w')
    8791074                i_m3_mode1w=i
     1075    type_tr(i_m3_mode1w) = 10
    8801076                PRINT*,'M3_m1w',i_m3_mode1w
    8811077                CASE('M3_m1ccn')
    8821078                i_m3_mode1ccn=i
     1079    type_tr(i_m3_mode1ccn) = 10
    8831080                PRINT*,'M3_m1ccn',i_m3_mode1ccn
    884                 CASE('M0_m2drop')
     1081    CASE('M0_m2drop')
    8851082                i_m0_mode2drop=i
     1083    type_tr(i_m0_mode2drop) = 10
    8861084                PRINT*,'M0_m2drop',i_m0_mode2drop
    887                 CASE('M0_m2ccn')
     1085    CASE('M0_m2ccn')
    8881086                i_m0_mode2ccn=i
     1087    type_tr(i_m0_mode2ccn) = 10
    8891088                PRINT*,'M0_m2ccn',i_m0_mode2ccn
    8901089                CASE('M3_m2sa')
    8911090                i_m3_mode2sa=i
     1091    type_tr(i_m3_mode2sa) = 10
    8921092                PRINT*,'M3_m2sa',i_m3_mode2sa
    8931093                CASE('M3_m2w')
    8941094                i_m3_mode2w=i
     1095    type_tr(i_m3_mode2w) = 10
    8951096                PRINT*,'M3_m2w',i_m3_mode2w
    8961097                CASE('M3_m2ccn')
    8971098                i_m3_mode2ccn=i
     1099    type_tr(i_m3_mode2ccn) = 10
    8981100                PRINT*,'M3_m2ccn',i_m3_mode2ccn
    8991101        END SELECT

  • trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h

    r2580 r2836  
    1010       LOGICAL cycle_diurne,soil_model
    1111       LOGICAL ok_orodr,ok_orolf,ok_gw_nonoro
    12        LOGICAL ok_kzmin,tuneupperatm
     12       LOGICAL ok_kzmin,tuneupperatm, ok_ionchem, ok_jonline
    1313       LOGICAL callnlte,callnirco2,callthermos
    1414       LOGICAL ok_cloud, ok_chem, reinit_trac, ok_sedim
    … …  
    2929     &     tuneupperatm,callnlte,callnirco2,callthermos,                &
    3030     &     ok_cloud, ok_chem, reinit_trac, ok_sedim,                    &
    31      &     ok_clmain, physideal, startphy_file
     31     &     ok_clmain, physideal, startphy_file, ok_ionchem, ok_jonline
    3232
    3333       COMMON/clesphys_i/ nbapp_rad, nbapp_chim,                        &

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

    r2686 r2836  
    6565         nbq = 0 ! to count number of tracers used in this subroutine
    6666
     67!        aki values comes from Aeronomy part B by P.M. BANKS / G. KOCKARTS – page 12-20-24)
     68
     69!        Ions are not here because the sum of all ions abundace is not
     70!        above 10^-4 until 250 km and we don't have their cpi and aki.
     71
     72!        Heat capacity for H, He, N, N2, O, O2, CO2, CO:
     73!        Circular of the bureau of Standards no. 564: tables of thermal properties of gases comprising [...]
     74!        Tables of Thermodynamic and Transport Properties of Air, Argon, Carbon Dioxide, Carbon Monoxide,
     75!        Hydrogen (molecular and atomic), Nitrogen (molecular and atomic), Oxygen (molecular and atomic), and Steam
     76!        https://www.govinfo.gov/content/pkg/GOVPUB-C13-89baf9f9b4a43e5f25820bd51b0f3f11/pdf/GOVPUB-C13-89baf9f9b4a43e5f25820bd51b0f3f11.pdf
     77       
     78
    6779         if (i_co2 /= 0) then
    6880            nbq = nbq + 1
    … …  
    144156            cpi(nbq) = 1.870e3
    145157         end if
    146 c         if (i_n /= 0) then
    147 c            nbq = nbq + 1
    148 c            niq(nbq) = i_n
    149 c            aki(nbq) = 0.0
    150 c            cpi(nbq) = 0.0
    151 c         endif
    152 c         if(i_no /= 0) then
    153 c            nbq = nbq + 1
    154 c            niq(nbq) = i_no
    155 c            aki(nbq) = 0.0
    156 c            cpi(nbq) = 0.0
    157 c         endif
    158 c         if(i_no2 /= 0) then
    159 c            nbq = nbq + 1
    160 c            niq(nbq) = i_no2
    161 c            aki(nbq) = 0.0
    162 c            cpi(nbq) = 0.0
    163 c         endif
    164 c         if(i_n2d /= 0) then
    165 c            nbq = nbq + 1
    166 c            niq(nbq) = i_n2d
    167 c            aki(nbq) = 0.0
    168 c            cpi(nbq) = 0.0
    169 c         endif
     158         if (i_he /= 0) then
     159            nbq = nbq + 1
     160            niq(nbq) = i_he
     161            aki(nbq) = 29.9e-4
     162            cpi(nbq) = 5.2e3
     163         end if
     164         if (i_n /= 0) then
     165            nbq = nbq + 1
     166            niq(nbq) = i_n
     167            aki(nbq) = 0.0
     168            cpi(nbq) = 1.4844e3
     169         endif
     170         if(i_no /= 0) then
     171            nbq = nbq + 1
     172            niq(nbq) = i_no
     173            aki(nbq) = 0.0
     174            cpi(nbq) = 0.0
     175         endif
     176         if(i_no2 /= 0) then
     177            nbq = nbq + 1
     178            niq(nbq) = i_no2
     179            aki(nbq) = 0.0
     180            cpi(nbq) = 0.0
     181         endif
     182         if(i_n2d /= 0) then
     183            nbq = nbq + 1
     184            niq(nbq) = i_n2d
     185            aki(nbq) = 0.0
     186            cpi(nbq) = 1.4844e3 !?
     187         endif
    170188         if (i_he /= 0) then
    171189            nbq = nbq + 1
    … …  
    255273            end do
    256274 
    257 
    258275            cpnew(ig,l) = cpnew(ig,l)/ntot
    259276            akknew(ig,l)= akknew(ig,l)/ntot

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

    r2580 r2836  
    485485  tuneupperatm = .false.
    486486  call getin('tuneupperatm',tuneupperatm)
     487
     488!
     489!Config Key  = ok_jonline
     490!Config Desc =
     491!Config Def  = .false.
     492!Config Help =
     493!
     494  ok_jonline = .false.
     495  call getin('ok_jonline',ok_jonline)
     496
     497!
     498!Config Key  = ok_ionchem
     499!Config Desc =
     500!Config Def  = .false.
     501!Config Help =
     502!
     503  ok_ionchem = .false.
     504  call getin('ok_ionchem',ok_ionchem)
     505
     506  if ((ok_jonline.eq..false.).and.(.true..eq.ok_ionchem)) then
     507    write(*,*) "Attention, incoherence :"
     508    write(*,*) "ok_jonline=",ok_jonline," / ok_ionchem=",ok_ionchem
     509    write(*,*) "Si vous souhaitez les ions, ok_jonline==.true."
     510    write(*,*) "Si vous ne voulez pas des ions, ok_ionchem==.false."
     511    write(*,*) "Verifiez votre physiq.def"
     512    stop
     513  endif
    487514
    488515!
    … …  
    547574  write(lunout,*)' euveff = ',euveff
    548575  write(lunout,*)' tuneupperatm = ',tuneupperatm
    549 
     576  write(lunout,*)' ok_jonline = ',ok_jonline
     577  write(lunout,*)' ok_ionchem = ',ok_ionchem
     578 
    550579  end subroutine conf_phys
    551580

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

    r2464 r2836  
    7676!!! If the values are changed there, the same has to be done here  !!!
    7777     
    78       integer,parameter :: ix_co2=1
    79       integer,parameter :: ix_o=3
    80       integer,parameter :: ix_co=4
    81       integer,parameter :: ix_n2=13
    82 
     78!      integer,parameter :: ix_co2=1
     79!      integer,parameter :: ix_o=3
     80!      integer,parameter :: ix_co=4
     81!      integer,parameter :: ix_n2=13
     82
     83      integer,parameter :: ix_co2  =  1
     84      integer,parameter :: ix_co   =  2
     85      integer,parameter :: ix_o    =  3
     86      integer,parameter :: ix_o1d  =  4
     87      integer,parameter :: ix_o2   =  5
     88      integer,parameter :: ix_o3   =  6
     89      integer,parameter :: ix_h    =  7
     90      integer,parameter :: ix_h2   =  8
     91      integer,parameter :: ix_oh   =  9
     92      integer,parameter :: ix_ho2  = 10
     93      integer,parameter :: ix_h2o2 = 11
     94      integer,parameter :: ix_h2o  = 12
     95      integer,parameter :: ix_n    = 13
     96      integer,parameter :: ix_n2d  = 14
     97      integer,parameter :: ix_no   = 15
     98      integer,parameter :: ix_no2  = 16
     99      integer,parameter :: ix_n2   = 17
    83100
    84101! Tracer indexes in the GCM:
    85102      integer,save :: g_co2=0
    86103      integer,save :: g_o=0
     104      integer,save :: g_o2=0
     105      integer,save :: g_h2=0
     106      integer,save :: g_h2o2=0
     107      integer,save :: g_h2o=0
     108      integer,save :: g_o3=0
     109      integer,save :: g_n2=0
     110      integer,save :: g_n=0
     111      integer,save :: g_no=0
    87112      integer,save :: g_co=0
    88       integer,save :: g_n2=0
     113      integer,save :: g_h=0
     114      integer,save :: g_no2=0
     115      integer,save :: g_oh=0
     116      integer,save :: g_ho2=0
     117      integer,save :: g_o1d=0
     118      integer,save :: g_n2d=0
    89119     
    90120      logical,save :: firstcall=.true.
    … …  
    110140            stop
    111141         else
    112            nspeuv_vgcm=nspeuv_vgcm+1
     142            nspeuv_vgcm=nspeuv_vgcm+1
    113143         endif
    114144         g_co=i_co
    … …  
    122152         ! n2
    123153         g_n2=i_n2
    124             if (g_n2.eq.0) then
    125                write(*,*) "euvheat: Error; no N2 tracer !!!"
    126 !               write(*,*) "N2 needed if NO is in traceur.def"
    127                stop
    128             else
    129                nspeuv_vgcm=nspeuv_vgcm+1
    130             endif
    131 
    132 !         g_o2=igcm_o2
    133 !         if (g_o2.eq.0) then
    134 !            write(*,*) "euvheat: Error; no O2 tracer !!!"
     154         if (g_n2.eq.0) then
     155            write(*,*) "euvheat: Error; no N2 tracer !!!"
     156!                write(*,*) "N2 needed if NO is in traceur.def"
     157            stop
     158         else
     159            nspeuv_vgcm=nspeuv_vgcm+1
     160         endif
     161         g_o2=i_o2
     162         if (g_o2.eq.0) then
     163            write(*,*) "euvheat: Error; no O2 tracer !!!"
    135164!            write(*,*) "O2 is always needed if calleuv=.true."
    136 !            stop
    137 !         else
    138 !            nespeuv=nespeuv+1
    139 !         endif
    140 !         g_h2=igcm_h2
    141 !         if (g_h2.eq.0) then
    142 !            write(*,*) "euvheat: Error; no H2 tracer !!!"
     165            stop
     166         else
     167            nspeuv_vgcm=nspeuv_vgcm+1
     168         endif
     169         g_h2=i_h2
     170         if (g_h2.eq.0) then
     171            write(*,*) "euvheat: Error; no H2 tracer !!!"
    143172!            write(*,*) "H2 is always needed if calleuv=.true."
    144 !            stop
    145 !         else
    146 !            nespeuv=nespeuv+1
    147 !         endif
    148 !         g_oh=igcm_oh
    149 !         if (g_oh.eq.0) then
    150 !            write(*,*) "euvheat: Error; no OH tracer !!!"
     173            stop
     174         else
     175            nspeuv_vgcm=nspeuv_vgcm+1
     176         endif
     177         g_oh=i_oh
     178         if (g_oh.eq.0) then
     179            write(*,*) "euvheat: Error; no OH tracer !!!"
    151180!            write(*,*) "OH must always be present if thermochem=T"
    152 !            stop
    153 !         else
    154 !            nespeuv=nespeuv+1 
    155 !         endif
    156 !         g_ho2=igcm_ho2
    157 !         if (g_ho2.eq.0) then
    158 !            write(*,*) "euvheat: Error; no HO2 tracer !!!"
     181            stop
     182         else
     183            nspeuv_vgcm=nspeuv_vgcm+1 
     184         endif
     185         g_ho2=i_ho2
     186         if (g_ho2.eq.0) then
     187            write(*,*) "euvheat: Error; no HO2 tracer !!!"
    159188!            write(*,*) "HO2 must always be present if thermochem=T"
    160 !            stop
    161 !         else
    162 !            nespeuv=nespeuv+1 
    163 !         endif
    164 !         g_h2o2=igcm_h2o2
    165 !         if (g_h2o2.eq.0) then
    166 !            write(*,*) "euvheat: Error; no H2O2 tracer !!!"
     189            stop
     190         else
     191            nspeuv_vgcm=nspeuv_vgcm+1 
     192         endif
     193         g_h2o2=i_h2o2
     194         if (g_h2o2.eq.0) then
     195            write(*,*) "euvheat: Error; no H2O2 tracer !!!"
    167196!            write(*,*) "H2O2 is always needed if calleuv=.true."
    168 !            stop
    169 !         else
    170 !            nespeuv=nespeuv+1
    171 !         endif
    172 !         g_h2o=igcm_h2o_vap
    173 !         if (g_h2o.eq.0) then
    174 !            write(*,*) "euvheat: Error; no water vapor tracer !!!"
     197            stop
     198         else
     199            nspeuv_vgcm=nspeuv_vgcm+1
     200         endif
     201         g_h2o=i_h2o
     202         if (g_h2o.eq.0) then
     203            write(*,*) "euvheat: Error; no water vapor tracer !!!"
    175204!            write(*,*) "H2O is always needed if calleuv=.true."
    176 !            stop
    177 !         else
    178 !            nespeuv=nespeuv+1
    179 !         endif
    180 !         g_o1d=igcm_o1d
    181 !         if (g_o1d.eq.0) then
    182 !            write(*,*) "euvheat: Error; no O1D tracer !!!"
     205            stop
     206         else
     207            nspeuv_vgcm=nspeuv_vgcm+1
     208         endif
     209         g_o1d=i_o1d
     210         if (g_o1d.eq.0) then
     211            write(*,*) "euvheat: Error; no O1D tracer !!!"
    183212!            write(*,*) "O1D must always be present if thermochem=T"
    184 !            stop
    185 !         else
    186 !            nespeuv=nespeuv+1 
    187 !         endif
    188 !         g_h=igcm_h
    189 !         if (g_h.eq.0) then
    190 !            write(*,*) "euvheat: Error; no H tracer !!!"
     213            stop
     214         else
     215            nspeuv_vgcm=nspeuv_vgcm+1 
     216         endif
     217         g_h=i_h
     218         if (g_h.eq.0) then
     219            write(*,*) "euvheat: Error; no H tracer !!!"
    191220!            write(*,*) "H is always needed if calleuv=.true."
    192 !            stop
    193 !         else
    194 !            nespeuv=nespeuv+1
    195 !         endif
     221            stop
     222         else
     223            nspeuv_vgcm=nspeuv_vgcm+1
     224         endif
    196225         
    197 !         euvmod = 3            !Default: C/O/H chemistry
     226!         euvmod = 1            !Default: C/O/H chemistry
    198227!         !Check if O3 is present
    199 !         g_o3=igcm_o3
    200 !         if (g_o3.eq.0) then
    201 !            write(*,*) "euvheat: Error; no O3 tracer !!!"
    202 !            write(*,*) "O3 must be present if calleuv=.true."
    203 !            stop
    204 !         else
    205 !            nespeuv=nespeuv+1
    206 !            euvmod=1
    207 !         endif
     228         g_o3=i_o3
     229         if (g_o3.eq.0) then
     230            write(*,*) "euvheat: Error; no O3 tracer !!!"
     231            write(*,*) "O3 must be present if calleuv=.true."
     232            stop
     233         else
     234            nspeuv_vgcm=nspeuv_vgcm+1
     235            euvmod=1
     236         endif
    208237
    209238         !Nitrogen species
    210239         !NO is used to determine if N chemistry is wanted
    211240         !euvmod=2 -> N chemistry
    212 !         g_no=igcm_no
    213 !         if (g_no.eq.0) then
    214 !            write(*,*) "euvheat: no NO tracer"
    215 !            write(*,*) "No N species in UV heating"
    216 !         else if(g_no.ne.0) then
    217 !            nespeuv=nespeuv+1
    218 !            euvmod=2
    219 !         endif
     241         g_no=i_no
     242         if (g_no.eq.0) then
     243            write(*,*) "euvheat: no NO tracer"
     244            write(*,*) "No N species in UV heating"
     245         else if(g_no.ne.0) then
     246            nspeuv_vgcm=nspeuv_vgcm+1
     247            euvmod=2
     248         endif
    220249         ! N
    221 !         g_n=igcm_n
    222 !         if(euvmod == 2) then
    223 !            if (g_n.eq.0) then
    224 !               write(*,*) "euvheat: Error; no N tracer !!!"
    225 !               write(*,*) "N needed if NO is in traceur.def"
    226 !               stop
    227 !            else if(g_n.ne.0) then
    228 !               nespeuv=nespeuv+1
    229 !            endif
    230 !         else
    231 !            if(g_n /= 0) then
    232 !               write(*,*) "euvheat: Error: N present, but NO not!!!"
    233 !               write(*,*) "Both must be in traceur.def"
    234 !               stop
    235 !            endif
    236 !         endif   !Of if(euvmod==2)
    237 !         !NO2
    238 !         g_no2=igcm_no2
    239 !         if(euvmod == 2) then
    240 !            if (g_no2.eq.0) then
    241 !               write(*,*) "euvheat: Error; no NO2 tracer !!!"
    242 !               write(*,*) "NO2 needed if NO is in traceur.def"
    243 !               stop
    244 !            else if(g_no2.ne.0) then
    245 !               nespeuv=nespeuv+1
    246 !            endif
    247 !         else
    248 !            if(g_no2 /= 0) then
    249 !               write(*,*) "euvheat: Error: NO2 present, but NO not!!!"
    250 !               write(*,*) "Both must be in traceur.def"
    251 !               stop
    252 !            endif
    253 !         endif   !Of if(euvmod==2)
     250         g_n=i_n
     251         if(euvmod == 2) then
     252            if (g_n.eq.0) then
     253               write(*,*) "euvheat: Error; no N tracer !!!"
     254               write(*,*) "N needed if NO is in traceur.def"
     255               stop
     256            else if(g_n.ne.0) then
     257               nspeuv_vgcm=nspeuv_vgcm+1
     258            endif
     259         else
     260            if(g_n /= 0) then
     261               write(*,*) "euvheat: Error: N present, but NO not!!!"
     262               write(*,*) "Both must be in traceur.def"
     263               stop
     264            endif
     265         endif   !Of if(euvmod==2)
     266         !NO2
     267         g_no2=i_no2
     268         if(euvmod == 2) then
     269            if (g_no2.eq.0) then
     270               write(*,*) "euvheat: Error; no NO2 tracer !!!"
     271               write(*,*) "NO2 needed if NO is in traceur.def"
     272               stop
     273            else if(g_no2.ne.0) then
     274               nspeuv_vgcm=nspeuv_vgcm+1
     275            endif
     276         else
     277            if(g_no2 /= 0) then
     278               write(*,*) "euvheat: Error: NO2 present, but NO not!!!"
     279               write(*,*) "Both must be in traceur.def"
     280               stop
     281            endif
     282         endif   !Of if(euvmod==2)
    254283         !N2D
    255 !         g_n2d=igcm_n2d
    256 !         if(euvmod == 2) then
    257 !            if (g_n2d.eq.0) then
    258 !               write(*,*) "euvheat: Error; no N2D tracer !!!"
    259 !               write(*,*) "N2D needed if NO is in traceur.def"
    260 !               stop
    261 !            else
    262 !               nespeuv=nespeuv+1 
    263 !            endif
    264 !         else
    265 !            if(g_n2d /= 0) then
    266 !               write(*,*) "euvheat: Error: N2D present, but NO not!!!"
    267 !               write(*,*) "Both must be in traceur.def"
    268 !               stop
    269 !            endif
    270 !         endif   !Of if(euvmod==2)
     284         g_n2d=i_n2d
     285         if(euvmod == 2) then
     286            if (g_n2d.eq.0) then
     287               write(*,*) "euvheat: Error; no N2D tracer !!!"
     288               write(*,*) "N2D needed if NO is in traceur.def"
     289               stop
     290            else
     291               nspeuv_vgcm=nspeuv_vgcm+1 
     292            endif
     293         else
     294            if(g_n2d /= 0) then
     295               write(*,*) "euvheat: Error: N2D present, but NO not!!!"
     296               write(*,*) "Both must be in traceur.def"
     297               stop
     298            endif
     299         endif   !Of if(euvmod==2)
    271300
    272301         !Check if nespeuv is appropriate for the value of euvmod
    … …  
    299328
    300329
    301       !Allocate density vector
    302       allocate(rm(nlev,nespeuv))
     330         !Allocate density vector
     331         allocate(rm(nlev,nespeuv))
    303332
    304333         firstcall= .false.
    … …  
    310339!cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    311340
     341      ! euvmod selection security
     342      if(euvmod.gt.2.or.euvmod.lt.0) then
     343         write(*,*)'euvheat: bad value for euvmod. Stop'
     344         stop
     345      endif
     346
    312347      ! build local updated values of tracers (if any) and temperature
    313 
    314348      do l=1,nlev
    315349        do ig=1,nlon
    316 
    317 ! chemical species
    318           zq(ig,l,g_co2)=pq(ig,l,g_co2) 
    319           zq(ig,l,g_co)=pq(ig,l,g_co)   
    320           zq(ig,l,g_o)=pq(ig,l,g_o)     
    321           zq(ig,l,g_n2)=pq(ig,l,g_n2)   
    322 
    323 ! atmospheric temperature
    324          zt(ig,l)=pt(ig,l)
    325 
     350          ! chemical species
     351          zq(ig,l,g_co2)=pq(ig,l,g_co2)     ! CO2 
     352          zq(ig,l,g_o2)=pq(ig,l,g_o2)       ! O2
     353          zq(ig,l,g_o)=pq(ig,l,g_o)         ! O(3P)
     354          zq(ig,l,g_h2)=pq(ig,l,g_h2)       ! H2
     355          zq(ig,l,g_h2o2)=pq(ig,l,g_h2o2)   ! H2O2
     356          zq(ig,l,g_h2o)=pq(ig,l,g_h2o)     ! H2O
     357          zq(ig,l,g_n2)=pq(ig,l,g_n2)       ! N2   
     358          zq(ig,l,g_co)=pq(ig,l,g_co)       ! CO
     359          zq(ig,l,g_h)=pq(ig,l,g_h)         ! H
     360
     361          !Only if O3, N or ion chemistry requested
     362          if(euvmod.ge.1) then
     363             zq(ig,l,g_o3)=pq(ig,l,g_o3)    ! 03
     364          endif
     365          !Only if N or ion chemistry requested
     366          if(euvmod.ge.2) then
     367             zq(ig,l,g_n)=pq(ig,l,g_n)      ! N
     368             zq(ig,l,g_no)=pq(ig,l,g_no)    ! NO
     369             zq(ig,l,g_no2)=pq(ig,l,g_no2)  ! NO2
     370          endif
     371          ! atmospheric temperature
     372          zt(ig,l)=pt(ig,l)
    326373!      write(*,*),  "CHECK update densities L332 euv",   zq(ig,l,g_co2)
    327 
    328 
    329374        enddo
    330375      enddo
    331376     
    332       !Solar flux calculation
    333      
     377      !Solar flux calculation     
    334378      do ig=1,nlon
    335379         zenit=acos(mu0(ig))*180./acos(-1.)  !convers from rad to deg
    336380         
    337          do l=1,nlev
    338          
    339 !   Conversion to number density
    340                    
    341 !!!  use R specific = R/MolarMass
    342 
     381         do l=1,nlev         
     382            !Conversion to number density   
     383            !!!  use R specific = R/MolarMass
    343384            dens=pplay(ig,l)/(rnew(ig,l)*zt(ig,l)) / 1.66e-21   ! [g mol-1] [cm-3]         
    344385
    345             rm(l,ix_co2)  = zq(ig,l,g_co2) * dens / M_tr(g_co2)   ! [cm-3]
     386            rm(l,ix_co2)  = zq(ig,l,g_co2) * dens / M_tr(g_co2) ! [cm-3]
     387            rm(l,ix_o2)   = zq(ig,l,g_o2)  * dens / M_tr(g_o2)
    346388            rm(l,ix_o)    = zq(ig,l,g_o)   * dens / M_tr(g_o)
    347             rm(l,ix_co)   = zq(ig,l,g_co) * dens  / M_tr(g_co)
    348             rm(l,ix_n2)   = zq(ig,l,g_n2) * dens  / M_tr(g_n2)
     389            rm(l,ix_h2)   = zq(ig,l,g_h2)  * dens / M_tr(g_h2)
     390            rm(l,ix_h2o)  = zq(ig,l,g_h2o) * dens / M_tr(g_h2o)
     391            rm(l,ix_h2o2) = zq(ig,l,g_h2o2)* dens / M_tr(g_h2o2)
     392            rm(l,ix_co)   = zq(ig,l,g_co)  * dens / M_tr(g_co)
     393            rm(l,ix_n2)   = zq(ig,l,g_n2)  * dens / M_tr(g_n2)
     394            rm(l,ix_h)    = zq(ig,l,g_h)   * dens / M_tr(g_h)
     395
     396            !Only if O3, N or ion chemistry requested
     397            if(euvmod.ge.1) then
     398               rm(l,ix_o3)   = zq(ig,l,g_o3) * dens / M_tr(g_o3)
     399            endif
     400            !Only if N or ion chemistry requested
     401            if(euvmod.ge.2) then
     402               rm(l,ix_n)    = zq(ig,l,g_n)    * dens / M_tr(g_n)
     403               rm(l,ix_no)   = zq(ig,l,g_no)   * dens / M_tr(g_no)
     404               rm(l,ix_no2)  = zq(ig,l,g_no2)  * dens / M_tr(g_no2)
     405            endif
    349406
    350407!           write(*,*),  "CHECK n density", l, rm(l,ix_co2)
    351 
    352 
    353408         enddo
    354409
    … …  
    368423         
    369424
     425                       !euveff: UV heating efficiency. Following Fox et al. ASR 1996
     426                       !should vary between 19% and 23%. Lower values
     427                       !(i.e. 16%) can be used to compensate
     428                       !underestimation
     429                       !of 15-um cooling (see Forget et al. JGR 2009 and
     430                       !Gonzalez-Galindo et al. JGR 2009) for details
     431
    370432        !Calculates the UV heating from the total photoabsorption coefficient
    371433        do l=1,nlev

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

    r2686 r2836  
    4242!!! If the values are changed there, the same has to be done here  !!!
    4343
    44       integer,parameter :: i_co2=1
    45       integer,parameter :: i_n2=13
    46       integer,parameter :: i_n=14
    47       integer,parameter :: i_o=3
    48       integer,parameter :: i_co=4
     44!      integer,parameter :: i_co2=1
     45!      integer,parameter :: i_n2=13
     46!      integer,parameter :: i_n=14
     47!      integer,parameter :: i_o=3
     48!      integer,parameter :: i_co=4
    4949
     50      integer,parameter :: ix_co2  =  1
     51      integer,parameter :: ix_co   =  2
     52      integer,parameter :: ix_o    =  3
     53      integer,parameter :: ix_o1d  =  4
     54      integer,parameter :: ix_o2   =  5
     55      integer,parameter :: ix_o3   =  6
     56      integer,parameter :: ix_h    =  7
     57      integer,parameter :: ix_h2   =  8
     58      integer,parameter :: ix_oh   =  9
     59      integer,parameter :: ix_ho2  = 10
     60      integer,parameter :: ix_h2o2 = 11
     61      integer,parameter :: ix_h2o  = 12
     62      integer,parameter :: ix_n    = 13
     63      integer,parameter :: ix_n2d  = 14
     64      integer,parameter :: ix_no   = 15
     65      integer,parameter :: ix_no2  = 16
     66      integer,parameter :: ix_n2   = 17
    5067
    5168c*************************PROGRAM STARTS*******************************
    … …  
    7390      !
    7491      do i=1,klev
    75          xabsi(1,i)  = rm(i,i_co2)
    76          xabsi(3,i)  = rm(i,i_o)
    77          xabsi(8,i)  = rm(i,i_n2)
    78          xabsi(11,i) = rm(i,i_co)
    79 
    80 c         xabsi(6,i)  = rm(i,i_h2o2)
     92         xabsi(1,i)  = rm(i,ix_co2)    ! CO2
     93         xabsi(2,i)  = rm(i,ix_o2)     ! O2
     94         xabsi(3,i)  = rm(i,ix_o)      ! O(3P)
     95         xabsi(4,i)  = rm(i,ix_h2o)    ! H2O
     96         xabsi(5,i)  = rm(i,ix_h2)     ! H2
     97         xabsi(6,i)  = rm(i,ix_h2o2)   ! H2O2
    8198         !Only if O3, N or ion chemistry requested
    82 c         if(euvmod.ge.1) then
    83 c            xabsi(7,i)  = rm(i,i_o3)
    84 c         endif
     99         if(euvmod.ge.1) then
     100            xabsi(7,i)  = rm(i,ix_o3)  ! O3
     101         endif
     102         xabsi(8,i)  = rm(i,ix_n2)     ! N2
    85103         !Only if N or ion chemistry requested
    86 c         if(euvmod.ge.2) then
    87 c            xabsi(8,i)  = rm(i,i_n2)
    88 c            xabsi(9,i)  = rm(i,i_n)
    89 c            xabsi(10,i) = rm(i,i_no)
    90 c            xabsi(13,i) = rm(i,i_no2)
    91 c         endif
     104         if(euvmod.ge.2) then
     105            xabsi(9,i)  = rm(i,ix_n)   ! N
     106            xabsi(10,i) = rm(i,ix_no)  ! NO
     107            xabsi(13,i) = rm(i,ix_no2) ! NO2
     108         endif
     109         xabsi(11,i) = rm(i,ix_co)     ! CO
     110         xabsi(12,i) = rm(i,ix_h)      ! H
    92111      end do
    93112

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

    r2809 r2836  
    9090      real*8      limup                        !        ""
    9191
    92       !!!ATTENTION. Here i_co2 has to have the same value than in euvheat.F90
     92      !!!ATTENTION. Here ix_co2 has to have the same value than in euvheat.F90
    9393      !!!If the value is changed there, if has to be changed also here !!!!
    94       integer,parameter :: i_co2=1
     94      integer,parameter :: ix_co2=1
    9595
    9696      character*20 modname
    … …  
    123123      !Calculation of column amounts
    124124      call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
    125      $     co2colx,o3pcolx,n2colx,cocolx)
     125     $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx,
     126     $     n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
    126127
    127128      !Auxiliar column to include the temperature dependence
    … …  
    130131      do i=nlayer-1,1,-1
    131132        coltemp(i)=!coltemp(i+1)+     PQ SE ELIMINA? REVISAR
    132      $         ( rm(i,i_co2) + rm(i+1,i_co2) ) * 0.5
     133     $         ( rm(i,ix_co2) + rm(i+1,ix_co2) ) * 0.5
    133134     $         * 1e5 * (iz(i+1)-iz(i)) * abs(t2(i)-t0(i))
    134135      end do
    … …  
    11151116
    11161117      subroutine column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
    1117      $     co2colx,o3pcolx, n2colx, cocolx)
    1118 
     1118     $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx,
     1119     $     n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
     1120
     1121c     Jun 2022        AM            readapted to Venus GCM
    11191122c     mar 2014        gg            adapted to Venus GCM
    11201123c     nov 2002        fgg           first version
    … …  
    11351138      integer    ig
    11361139      integer    chemthermod
    1137       integer    nesptherm                   !# of species undergoing
    1138 chemistry, input
     1140      integer    nesptherm                  !# of species undergoing chemistry, input
    11391141      real       rm(klev,nesptherm)         !densities (cm-3), input
    11401142      real       tx(klev)                   !temperature profile, input
    … …  
    11421144      real       zenit                      !SZA, input
    11431145      real       co2colx(klev)              !column density of CO2 (cm^-2), output
     1146      real       o2colx(klev)               !column density of O2(cm^-2), output
    11441147      real       o3pcolx(klev)              !column density of O(3P)(cm^-2), output
     1148      real       h2colx(klev)               !H2 column density (cm-2), output
     1149      real       h2ocolx(klev)              !H2O column density (cm-2), output
     1150      real       h2o2colx(klev)             !column density of H2O2(cm^-2), output
     1151      real       o3colx(klev)               !O3 column density (cm-2), output
    11451152      real       n2colx(klev)               !N2 column density (cm-2), output
     1153      real       ncolx(klev)                !N column density (cm-2), output
     1154      real       nocolx(klev)               !NO column density (cm-2), output
    11461155      real       cocolx(klev)               !CO column density (cm-2), output
     1156      real       hcolx(klev)                !H column density (cm-2), output
     1157      real       no2colx(klev)              !NO2 column density (cm-2), output
     1158
    11471159
    11481160c     local variables
    … …  
    11531165      real       Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2
    11541166
     1167      ! density
    11551168      real       co2x(klev)
     1169      real       o2x(klev)
    11561170      real       o3px(klev)
    11571171      real       cox(klev)
     1172      real       hx(klev)
     1173      real       h2x(klev)
     1174      real       h2ox(klev)
     1175      real       h2o2x(klev)
     1176      real       o3x(klev)
    11581177      real       n2x(klev)
    11591178      real       nx(klev)
     1179      real       nox(klev)
     1180      real       no2x(klev)
    11601181
    11611182      integer    i,j,k,icol,indexint          !indexes
    … …  
    11791200!!! If the values are changed there, the same has to be done here  !!!
    11801201
    1181       integer,parameter :: ix_co2=1
    1182       integer,parameter :: ix_n2=13
    1183       integer,parameter :: ix_o=3
    1184       integer,parameter :: ix_co=4
     1202!      integer,parameter :: ix_co2=1
     1203!      integer,parameter :: ix_n2=13
     1204!      integer,parameter :: ix_o=3
     1205!      integer,parameter :: ix_co=4
     1206
     1207      integer,parameter :: ix_co2  =  1
     1208      integer,parameter :: ix_co   =  2
     1209      integer,parameter :: ix_o    =  3
     1210      integer,parameter :: ix_o1d  =  4
     1211      integer,parameter :: ix_o2   =  5
     1212      integer,parameter :: ix_o3   =  6
     1213      integer,parameter :: ix_h    =  7
     1214      integer,parameter :: ix_h2   =  8
     1215      integer,parameter :: ix_oh   =  9
     1216      integer,parameter :: ix_ho2  = 10
     1217      integer,parameter :: ix_h2o2 = 11
     1218      integer,parameter :: ix_h2o  = 12
     1219      integer,parameter :: ix_n    = 13
     1220      integer,parameter :: ix_n2d  = 14
     1221      integer,parameter :: ix_no   = 15
     1222      integer,parameter :: ix_no2  = 16
     1223      integer,parameter :: ix_n2   = 17
    11851224
    11861225c*************************PROGRAM STARTS*******************************
    … …  
    11951234      xx = kboltzman * tx(klev) * n_avog / grav(klev)   ! g cm mol-1
    11961235
    1197       Ho3p  = xx / mmolo
    11981236      Hco2  = xx / mmolco2
     1237      Ho2   = xx / mmolo2
     1238      Ho3p  = xx / mmolo     ! Oxygen 3P
     1239      Hh2   = xx / mmolh2
     1240      Hh2o2 = xx / mmolh2o2
     1241      Hh2o  = xx / mmolh2o     
     1242
     1243      !Only if O3 chem. required
     1244      if(chemthermod.ge.1)
     1245     $     Ho3   = xx / mmolo3
     1246      Hn2   = xx / mmoln2
     1247      !Only if N or ion chem.
     1248      if(chemthermod.ge.2) then
     1249         Hn    = xx / mmoln
     1250         Hno   = xx / mmolno
     1251         Hno2  = xx / mmolno2
     1252      endif
    11991253      Hco   = xx / mmolco
    1200       Hn2   = xx / mmoln2
    1201       Hn    = xx / mmoln
    1202 
    1203       !Only if O3 chem. required
    1204 c      if(chemthermod.ge.1)
    1205 !     $     Ho3   = xx / mmol(igcm_o3)
    1206 c     $     Ho3   = xx / mmolo3
    1207 c      !Only if N or ion chem.
    1208 c      if(chemthermod.ge.2) then
    1209 c         Hn2   = xx / mmoln2
    1210 c         Hn    = xx / mmoln
    1211 c         Hno   = xx / mmolno
    1212 c         Hno2  = xx / mmolno2
    1213 c      endif
     1254      Hh    = xx / mmolh
    12141255      ! first loop in altitude : initialisation
    12151256      do i=klev,1,-1
    12161257         !Column initialisation
    12171258         co2colx(i)  = 0.
     1259         o2colx(i)   = 0.
    12181260         o3pcolx(i)  = 0.
     1261         h2colx(i)   = 0.
     1262         h2ocolx(i)  = 0.
     1263         h2o2colx(i) = 0.
     1264         o3colx(i)   = 0.
    12191265         n2colx(i)   = 0.
     1266         ncolx(i)    = 0.
     1267         nocolx(i)   = 0.
    12201268         cocolx(i)   = 0.
    1221 
     1269         hcolx(i)    = 0.
     1270         no2colx(i)  = 0.
    12221271         !--Densities [cm-3]
    12231272         co2x(i)  = rm(i,ix_co2)
     1273         o2x(i)   = rm(i,ix_o2)
    12241274         o3px(i)  = rm(i,ix_o)
     1275         h2x(i)   = rm(i,ix_h2)
     1276         h2ox(i)  = rm(i,ix_h2o)
     1277         h2o2x(i) = rm(i,ix_h2o2)
    12251278         cox(i)   = rm(i,ix_co)
     1279         hx(i)    = rm(i,ix_h)  ! write(*,*), '--jthermcalc--', co2x(i)
     1280
     1281         !Only if O3 chem. required
     1282         if(chemthermod.ge.1)
     1283     $        o3x(i)   = rm(i,ix_o3)
    12261284         n2x(i)   = rm(i,ix_n2)
    1227 
    1228 c         write(*,*), '--jthermcalc--', co2x(i)
    1229 
    1230          !Only if O3 chem. required
    1231 c         if(chemthermod.ge.1)
    1232 c     $        o3x(i)   = rm(i,i_o3)
    1233 c         !Only if Nitrogen of ion chemistry requested
    1234 c         if(chemthermod.ge.2) then
    1235 c            n2x(i)   = rm(i,i_n2)
    1236 c            nx(i)    = rm(i,i_n)
    1237 c            nox(i)   = rm(i,i_no)
    1238 c            no2x(i)  = rm(i,i_no2)
    1239 c         endif
     1285         !Only if Nitrogen of ion chemistry requested
     1286         if(chemthermod.ge.2) then
     1287            nx(i)    = rm(i,ix_n)
     1288            nox(i)   = rm(i,ix_no)
     1289            no2x(i)  = rm(i,ix_no2)
     1290         endif
    12401291      enddo     ! end first loop
    12411292
    … …  
    12471298         if(ilayesp(nlayesp).eq.-1) then
    12481299            co2colx(i)=1.e25
     1300            o2colx(i)=1.e25
    12491301            o3pcolx(i)=1.e25
     1302            h2colx(i)=1.e25
     1303            h2ocolx(i)=1.e25
     1304            h2o2colx(i)=1.e25
     1305            o3colx(i)=1.e25
    12501306            n2colx(i)=1.e25
    12511307            cocolx(i)=1.e25
    1252 
    1253 c            o2colx(i)=1.e25
    1254 c            o3pcolx(i)=1.e25
    1255 c            h2colx(i)=1.e25
    1256 c            h2ocolx(i)=1.e25
    1257 c            h2o2colx(i)=1.e25
    1258 c            o3colx(i)=1.e25
    1259 c            ncolx(i)=1.e25
    1260 c            nocolx(i)=1.e25
    1261 
    1262 c            cocolx(i)=1.e25
    1263 c            hcolx(i)=1.e25
    1264 c            no2colx(i)=1.e25
     1308            hcolx(i)=1.e25
     1309            ncolx(i)=1.e25
     1310            nocolx(i)=1.e25
     1311            no2colx(i)=1.e25
    12651312         else
    12661313            rcmnz = ( radio + iz(klev) ) * 1.e5    ! km --> cm
    12671314            rcmmini = ( radio + zmini ) * 1.e5
    1268             !Column calculation taking into account the geometrical
    1269             !depth
     1315            !Column calculation taking into account the geometrical depth
    12701316            !calculated before
    12711317            do j=1,nlayesp
    … …  
    12781324                     co2colx(i)=co2colx(i)+co2x(klev)*Hco2*esp(j)
    12791325     $                    *1.e-5
     1326                     h2o2colx(i)=h2o2colx(i)+
     1327     $                    h2o2x(klev)*Hh2o2*esp(j)*1.e-5
     1328                     o2colx(i)=o2colx(i)+o2x(klev)*Ho2*esp(j)
     1329     $                    *1.e-5
     1330                     h2colx(i)=h2colx(i)+h2x(klev)*Hh2*esp(j)
     1331     $                    *1.e-5
     1332                     h2ocolx(i)=h2ocolx(i)+h2ox(klev)*Hh2o*esp(j)
     1333     $                    *1.e-5
     1334                     n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)
     1335     $                    *1.e-5                     
    12801336                     cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j)
    12811337     $                    *1.e-5
    1282                      n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)
     1338                     hcolx(i)=hcolx(i)+hx(klev)*Hh*esp(j)
    12831339     $                    *1.e-5
    12841340
    1285 c                     h2o2colx(i)=h2o2colx(i)+
    1286 c     $                    h2o2x(klev)*Hh2o2*esp(j)*1.e-5
    1287 c                     o2colx(i)=o2colx(i)+o2x(klev)*Ho2*esp(j)
    1288 c     $                    *1.e-5
    1289 c                     h2colx(i)=h2colx(i)+h2x(klev)*Hh2*esp(j)
    1290 c     $                    *1.e-5
    1291 c                     h2ocolx(i)=h2ocolx(i)+h2ox(klev)*Hh2o*esp(j)
    1292 c     $                    *1.e-5                     
    1293 c                     cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j)
    1294 c     $                    *1.e-5
    1295 c                     hcolx(i)=hcolx(i)+hx(klev)*Hh*esp(j)
    1296 c     $                    *1.e-5
    12971341                     !Only if O3 chemistry required
    1298 c                     if(chemthermod.ge.1) o3colx(i)=
    1299 c     $                    o3colx(i)+o3x(klev)*Ho3*esp(j)
    1300 c     $                    *1.e-5
     1342                     if(chemthermod.ge.1) o3colx(i)=
     1343     $                    o3colx(i)+o3x(klev)*Ho3*esp(j)
     1344     $                    *1.e-5
    13011345                     !Only if N or ion chemistry requested
    1302 c                     if(chemthermod.ge.2) then
    1303 c                        n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)
    1304 c     $                    *1.e-5
    1305 
    1306 c                     endif
     1346                     if(chemthermod.ge.2) then
     1347                        ncolx(i)=ncolx(i)+nx(klev)*Hn*esp(j)
     1348     $                       *1.e-5
     1349                        nocolx(i)=nocolx(i)+nox(klev)*Hno*esp(j)
     1350     $                       *1.e-5
     1351                        no2colx(i)=no2colx(i)+no2x(klev)*Hno2*esp(j)
     1352     $                       *1.e-5
     1353                     endif
    13071354                  else if(zenit.gt.60.) then
    13081355                     espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j)
     1356                     espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
    13091357                     espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j)
     1358                     esph2  = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j)
     1359                     esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j)
     1360                     esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j)
    13101361                     espco  = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j)
    13111362                     espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
    1312                      espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
    1313 
    1314 c                     espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
    1315 c                     esph2  = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j)
    1316 c                     esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j)
    1317 c                     esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j)
    1318 c                     esph   = sqrt((rcmnz+Hh)**2 -rcmmini**2)  - esp(j)
     1363                     esph   = sqrt((rcmnz+Hh)**2 -rcmmini**2)  - esp(j)
     1364
    13191365                     !Only if O3 chemistry required
    1320 c                     if(chemthermod.ge.1)                     
    1321 c     $                   espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j)
    1322 c                     !Only if N or ion chemistry requested
    1323 c                     if(chemthermod.ge.2) then
    1324 c                        espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
    1325 c                        espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
    1326 c                        espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j)
    1327 c                        espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j)
    1328 c                     endif
     1366                     if(chemthermod.ge.1)                     
     1367     $                   espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j)
     1368                     !Only if N or ion chemistry requested
     1369                     if(chemthermod.ge.2) then
     1370                        espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
     1371                        espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j)
     1372                        espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j)
     1373                     endif
    13291374
    13301375                     co2colx(i) = co2colx(i) + espco2*co2x(klev)
     1376                     o2colx(i)  = o2colx(i)  + espo2*o2x(klev)
    13311377                     o3pcolx(i) = o3pcolx(i) + espo3p*o3px(klev)
     1378                     n2colx(i)  = n2colx(i)  + espn2*n2x(klev)
     1379                     h2colx(i)  = h2colx(i)  + esph2*h2x(klev)
     1380                     h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(klev)
     1381                     h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(klev)
    13321382                     cocolx(i)  = cocolx(i)  + espco*cox(klev)
    1333                      n2colx(i)  = n2colx(i)  + espn2*n2x(klev)
    1334 
    1335 c                     o2colx(i)  = o2colx(i)  + espo2*o2x(klev)
    1336 c                     h2colx(i)  = h2colx(i)  + esph2*h2x(klev)
    1337 c                     h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(klev)
    1338 c                     h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(klev)
    1339 c                     cocolx(i)  = cocolx(i)  + espco*cox(klev)
    1340 c                     hcolx(i)   = hcolx(i)   + esph*hx(klev)
     1383                     hcolx(i)   = hcolx(i)   + esph*hx(klev)
     1384
    13411385                     !Only if O3 chemistry required
    1342 c                     if(chemthermod.ge.1)                     
    1343 c     $                  o3colx(i) = o3colx(i)  + espo3*o3x(klev)
    1344 c                     !Only if N or ion chemistry requested
    1345 c                     if(chemthermod.ge.2) then
    1346 c                        n2colx(i)  = n2colx(i)  + espn2*n2x(klev)
    1347 c                        ncolx(i)   = ncolx(i)   + espn*nx(klev)
    1348 c                        nocolx(i)  = nocolx(i)  + espno*nox(klev)
    1349 c                        no2colx(i) = no2colx(i) + espno2*no2x(klev)
    1350 c                     endif
     1386                     if(chemthermod.ge.1)                     
     1387     $                  o3colx(i) = o3colx(i)  + espo3*o3x(klev)
     1388                     !Only if N or ion chemistry requested
     1389                     if(chemthermod.ge.2) then
     1390                        ncolx(i)   = ncolx(i)   + espn*nx(klev)
     1391                        nocolx(i)  = nocolx(i)  + espno*nox(klev)
     1392                        no2colx(i) = no2colx(i) + espno2*no2x(klev)
     1393                     endif
    13511394                  endif !Of if zenit.lt.60
    13521395               !Other layers
    13531396               else
    1354                   co2colx(i)  = co2colx(i) +
     1397                  co2colx(i)  = co2colx(i) + 
    13551398     $                 esp(j) * (co2x(jj)+co2x(jj+1)) / 2.
    1356                   o3pcolx(i)  = o3pcolx(i) +
     1399                  o2colx(i)   = o2colx(i) +
     1400     $                 esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
     1401                  o3pcolx(i)  = o3pcolx(i) +
    13571402     $                 esp(j) * (o3px(jj)+o3px(jj+1)) / 2.
    1358                   cocolx(i)   = cocolx(i) +
    1359      $                 esp(j) * (cox(jj)+cox(jj+1)) / 2.
     1403                  h2colx(i)   = h2colx(i) +
     1404     $                 esp(j) * (h2x(jj)+h2x(jj+1)) / 2.
     1405                  h2ocolx(i)  = h2ocolx(i) +
     1406     $                 esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2.
     1407                  h2o2colx(i) = h2o2colx(i) +
     1408     $                 esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2.
    13601409                  n2colx(i)   = n2colx(i) +
    13611410     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
    1362 
    1363 c
    1364 c                  o2colx(i)   = o2colx(i) +
    1365 c     $                 esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
    1366 c                  h2colx(i)   = h2colx(i) +
    1367 c     $                 esp(j) * (h2x(jj)+h2x(jj+1)) / 2.
    1368 c                  h2ocolx(i)  = h2ocolx(i) +
    1369 c     $                 esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2.
    1370 c                  h2o2colx(i) = h2o2colx(i) +
    1371 c     $                 esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2.
    1372 c                  hcolx(i)    = hcolx(i) +
    1373 c     $                 esp(j) * (hx(jj)+hx(jj+1)) / 2.
     1411                  cocolx(i)   = cocolx(i) +
     1412     $                 esp(j) * (cox(jj)+cox(jj+1)) / 2.
     1413                  hcolx(i)    = hcolx(i) +
     1414     $                 esp(j) * (hx(jj)+hx(jj+1)) / 2.
     1415
    13741416                  !Only if O3 chemistry required
    1375 c                  if(chemthermod.ge.1)
    1376 c     $                 o3colx(i) = o3colx(i) +
    1377 c     $                 esp(j) * (o3x(jj)+o3x(jj+1)) / 2.
    1378 c                  !Only if N or ion chemistry requested
    1379 c                  if(chemthermod.ge.2) then
    1380 c                     n2colx(i)   = n2colx(i) +
    1381 c     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
    1382 c                     ncolx(i)    = ncolx(i) +
    1383 c     $                    esp(j) * (nx(jj)+nx(jj+1)) / 2.
    1384 c                     nocolx(i)   = nocolx(i) +
    1385 c     $                    esp(j) * (nox(jj)+nox(jj+1)) / 2.
    1386 c                     no2colx(i)  = no2colx(i) +
    1387 c     $                    esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
    1388 c                  endif
    1389 
     1417                  if(chemthermod.ge.1)
     1418     $                 o3colx(i) = o3colx(i) +
     1419     $                 esp(j) * (o3x(jj)+o3x(jj+1)) / 2.
     1420                  !Only if N or ion chemistry requested
     1421                  if(chemthermod.ge.2) then
     1422                     ncolx(i)    = ncolx(i) +
     1423     $                    esp(j) * (nx(jj)+nx(jj+1)) / 2.
     1424                     nocolx(i)   = nocolx(i) +
     1425     $                    esp(j) * (nox(jj)+nox(jj+1)) / 2.
     1426                     no2colx(i)  = no2colx(i) +
     1427     $                    esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
     1428                  endif
    13901429               endif  !Of if jj.eq.klev
    13911430            end do    !Of do j=1,nlayesp
    … …  
    14061445C escout(nlayer) on pressure grid p(nlayer).
    14071446C
    1408       real*8 wm(nlayer),wp(nlayer),p(nlayer)
    1409       integer nm(nlayer)
    1410       real*8 pin(nl)
    1411       real*8 limup,limdown
    1412       integer nl,nlayer,n1,n,np,nini
     1447      real*8,intent(out) :: wm(nlayer),wp(nlayer) ! interpolation weights
     1448      integer,intent(out) :: nm(nlayer) ! index of nearest point
     1449      real*8,intent(in) :: pin(nl),p(nlayer)
     1450      real*8,intent(in) :: limup,limdown
     1451      integer,intent(in) :: nl,nlayer
     1452      integer :: n1,n,np,nini
    14131453      nini=1
    14141454      do n1=1,nlayer
    … …  
    14511491      real        szadeg                ! I. SZA [rad]
    14521492      real        z                     ! I. altitude of interest [km]
    1453       integer     nz3,ig                   ! I. dimension of esp, ylayesp, etc...
    1454                                         !  (=2*klev= max# of layers in
    1455                                         !  ray path)
    1456       real     iz(klev+1)              ! I. Altitude of each layer
     1493      integer     nz3,ig                ! I. dimension of esp, ylayesp, etc...
     1494                                        !  (=2*klev= max# of layers in ray path)
     1495      real     iz(klev+1)                ! I. Altitude of each layer
    14571496      real*8        esp(nz3)            ! O. layer widths after geometrically
    1458                                         !    amplified; in [cm] except
    1459                                         !    at TOA
    1460                                         !    where an auxiliary value is
    1461                                         !    used
     1497                                        !    amplified; in [cm] except at TOA
     1498                                        !    where an auxiliary value is  used
    14621499      real*8        ilayesp(nz3)        ! O. Indexes of layers along ray path
    14631500      real*8        szalayesp(nz3)      ! O. SZA [deg]    "     "     "
    14641501      integer       nlayesp
    1465 !      real*8        nlayesp             ! O. # layers along ray path at
    1466 !      this z
     1502!      real*8        nlayesp             ! O. # layers along ray path at this z
    14671503      real*8        zmini               ! O. Minimum altitud of ray path [km]
    14681504
    … …  
    15011537        ! The optical thickness will be given by  1/cos(sza)
    15021538        ! We deal with 2 different regions:
    1503         !   1: First, all layers between z and ztop ("upper part of
    1504         !   ray")
     1539        !   1: First, all layers between z and ztop ("upper part of ray")
    15051540        !   2: Second, the layer at ztop
    15061541        if(szadeg.lt.60.d0) then
    … …  
    15291564        ! The optical thickness is evaluated.
    15301565        !    (the magnitude of the effect of not using cos(sza) is 3.e-5
    1531         !     for z=60km & sza=30, and 5e-4 for z=60km & sza=60,
    1532         !     approximately)
     1566        !     for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately)
    15331567        ! We deal with 2 different regions:
    1534         !   1: First, all layers between z and ztop ("upper part of
    1535         !   ray")
     1568        !   1: First, all layers between z and ztop ("upper part of ray")
    15361569        !   2: Second, the layer at ztop ("uppermost layer")
    15371570        else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
    … …  
    15531586              rkmj = radio+diz(j)
    15541587              szalayesp(nlayesp) = asin( rkmmini/rkmj )             ! [rad]
    1555               szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592
    1556 ! [deg]
     1588              szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
    15571589           end do
    15581590 1470      continue
    … …  
    16031635                rkmj = radio+diz(j)
    16041636                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
    1605                 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592
    1606 ! [deg]
     1637                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592! [deg]
    16071638              end do
    16081639
    … …  
    16291660                rkmj = radio+diz(j+1)
    16301661                szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
    1631                 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592
    1632 ! [deg]
     1662                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592! [deg]
    16331663              endif
    16341664
    … …  
    16441674                rkmj = radio+diz(j)
    16451675                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
    1646                 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592
    1647 ! [deg]
     1676                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592! [deg]
    16481677              end do
    16491678
    … …  
    16551684                esp(nlayesp) =
    16561685     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
    1657      $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )
    1658 ! [km]
    1659                 esp(nlayesp) = esp(nlayesp) * 1.d5
    1660 ! [cm]
     1686     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )! [km]
     1687                esp(nlayesp) = esp(nlayesp) * 1.d5! [cm]
    16611688                rkmj = radio+diz(j)
    1662                 szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )
    1663 ! [rad]
    1664                 szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592
    1665 ! [deg]
     1689                szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )! [rad]
     1690                szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592! [deg]
    16661691              end do
    16671692

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

    r2795 r2836  
    421421! DIffusion coefficient
    422422        CALL DCOEFF(nn,dij,Praf,Traf,Nraf,Nrafk,Draf,nlraf,ncompdiff)
     423        !Draf(:) = max(100.,Draf(:))
     424        !Draf(:) = 0.01*Draf(:)
    423425        Drafmol(:,nn)=Draf(:)
    424426! Scale height of the density of the specie

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

    r2686 r2836  
    4242
    4343     USE ioipsl_getin_p_mod, ONLY : getin_p
    44      use geometry_mod, only: cell_area
     44     USE geometry_mod, ONLY: cell_area, latitude_deg
    4545!#ifdef CPP_XIOS
    4646!     use xios_output_mod, only: send_xios_field
    … …  
    107107     REAL, parameter:: cmin = 1.    ! Min horizontal absolute phase velocity
    108108
    109 !!! Nominal as Gilli2017
    110      REAL, parameter:: epflux_max = 0.005  ! Max EP flux value at launching altitude (previous name: RUWMAX)
     109!----------------------------------------
     110! Nominal as Gilli2017
     111!     REAL, parameter:: epflux_max = 0.005  ! Max EP flux value at launching altitude (previous name: RUWMAX)
     112!----------------------------------------
     113! VCD 2.1 tuning
     114     REAL, parameter:: epflux_max = 0.001  ! Max EP flux value at launching altitude (previous name: RUWMAX)
     115!----------------------------------------
    111116
    112117     REAL cpha                      ! absolute phase velocity frequency
    … …  
    251256     DO jw = 1, nw
    252257        DO ii = 1, ngrid
    253            
     258            
    254259           ran_num_1 = mod(tt(ii, jw) * 10., 1.)
    255260           ran_num_2 = mod(tt(ii, jw) * 100., 1.)
    256261
    257262!! OPTIONS GENERIC DIFF VENUS !!
    258            ! angle (random)
    259            zp(jw, ii) = (sign(1., 0.5 - ran_num_1) &
     263           ! angle (random) - reference
     264             zp(jw, ii) = (sign(1., 0.5 - ran_num_1) &
    260265                        + 1.) * RPI / 2.
    261            ! angle (0 or RPI)
    262            !zp(jw, ii) = RPI*mod(jw,2)
     266           ! --------- TRY 00 -----------------------
     267           !IF((abs(latitude_deg(ii)) .le. 55.)) THEN
     268           !  zp(jw, ii) = 0. + 2. * ( ran_num_1 - 0.5 ) * RPI /180. &
     269           !                       * (10.) ! +/- 10 deg par rapport equateur
     270           !ELSE IF ((abs(latitude_deg(ii)) .le. 75.)) THEN
     271           !  zp(jw, ii) = 0. + 2. * ( ran_num_1 - 0.5 ) * RPI /180. &
     272           !                       * (10. +(90.-10.)                 &
     273           !                       * (latitude_deg(ii)-55.)/20.)
     274           !ELSE
     275           !  zp(jw, ii) = ran_num_1 * RPI
     276           !ENDIF
     277           !zp(jw, ii) = mod(zp(jw, ii),2.*RPI)
     278           ! ------ TRY 01-------------------------------
     279           !IF((abs(latitude_deg(ii)) .le. 55)) THEN
     280           !  zp(jw, ii) = (sign(1., 0.5 - ran_num_1) &
     281           !             + 1.) * RPI / 2.
     282           !ELSE
     283           !  zp(jw, ii) = ran_num_1 * RPI
     284           !ENDIF
     285           ! ---- angle (0 or RPI) -----
     286           !!zp(jw, ii) = RPI*mod(jw,2)
    263287
    264288           ! horizontal wavenumber amplitude

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

    r2464 r2836  
    4747       write(*,*)'It should be in :', trim(datafile),'/'
    4848       write(*,*)'1) You can change this directory address in '
    49        write(*,*)'   callphys.def with datadir=/path/to/dir'
     49       write(*,*)'   callphys.def with datafile=/path/to/dir'
    5050       write(*,*)'2) If necessary, EUVDAT (and other datafiles)'
    5151       write(*,*)'   can be obtained online on:'
    … …  
    279279c     dissociation and ionization efficiencies
    280280
    281 !      do inter=1,ninter
    282 !         efdisco2(inter)=0.
    283 !         efdiso2(inter)=0.
    284 !         efdish2(inter)=0.
    285 !         efdish2o(inter)=0.
    286 !         efdish2o2(inter)=0.
    287 !         efdiso3(inter)=0.
    288 !         efdisco(inter)=0.
    289 !         efdisn2(inter)=0.
    290 !         efdisno(inter)=0.
    291 !         efdisno2(inter)=0.
    292 !         efionco2(inter,1)=0.
    293 !         efionco2(inter,2)=0.
    294 !         efionco2(inter,3)=0.
    295 !         efionco2(inter,4)=0.
    296 !         efiono3p(inter)=0.
    297 !         efionn2(inter,1)=0.
    298 !         efionn2(inter,2)=0.
    299 !         efionco(inter,1)=0.
    300 !         efionco(inter,2)=0.
    301 !         efionco(inter,3)=0.
    302 !         efionn(inter)=0.
    303 !         efionh(inter)=0.
    304 !         efionno(inter)=0.
    305 !      enddo
     281      do inter=1,ninter
     282         efdisco2(inter)=0.
     283         efdiso2(inter)=0.
     284         efdish2(inter)=0.
     285         efdish2o(inter)=0.
     286         efdish2o2(inter)=0.
     287         efdiso3(inter)=0.
     288         efdisco(inter)=0.
     289         efdisn2(inter)=0.
     290         efdisno(inter)=0.
     291         efdisno2(inter)=0.
     292         efionco2(inter,1)=0.
     293         efionco2(inter,2)=0.
     294         efionco2(inter,3)=0.
     295         efionco2(inter,4)=0.
     296         efiono2(inter,1)=0.
     297         efiono2(inter,2)=0.
     298         efiono3p(inter)=0.
     299         efionn2(inter,1)=0.
     300         efionn2(inter,2)=0.
     301         efionco(inter,1)=0.
     302         efionco(inter,2)=0.
     303         efionco(inter,3)=0.
     304         efionn(inter)=0.
     305         efionh(inter)=0.
     306         efionno(inter)=0.
     307      enddo
    306308
    307309
    … …  
    320322!         efdisno(inter)=1.
    321323!      enddo
     324!      close(120)
    322325
    323326
    324327c     N2
    325328
    326 !      efdisn2(15)=0.1
    327 !      do inter=16,ninter
    328 !         efdisn2(inter)=1.
    329 !      enddo
     329      efdisn2(15)=0.1
     330      do inter=16,ninter
     331         efdisn2(inter)=1.
     332      enddo
    330333
    331334
    332335c     CO
    333336
    334 !      efdisco(16)=0.5
    335 !      do inter=17,ninter
    336 !         efdisco(inter)=1.
    337 !      enddo
     337      efdisco(16)=0.5
     338      do inter=17,ninter
     339         efdisco(inter)=1.
     340      enddo
    338341
    339342     
    340343c     O, N, H
    341344
    342 !      do inter=1,ninter
    343 !         efdiso(inter)=0.
    344 !         efdisn(inter)=0.
    345 !         efdish(inter)=0.
    346 !      enddo
     345      do inter=1,ninter
     346         efdiso(inter)=0.
     347         efdisn(inter)=0.
     348         efdish(inter)=0.
     349      enddo
    347350
    348351
    349352c     H2O, H2O2, O3, NO2
    350353
    351 !      do inter=25,31
    352 !         efdish2o(inter)=1.
    353 !      enddo
    354 !      do inter=25,35
    355 !         efdish2o2(inter)=1.
    356 !      enddo
    357 !      do inter=34,36
    358 !         efdiso3(inter)=1.
    359 !      enddo
    360 !      do inter=27,36
    361 !         efdisno2(inter)=1.
    362 !      enddo
    363 !      do inter=1,15
    364 !         efdish2(inter)=1.
    365 !      enddo
     354      do inter=25,31
     355         efdish2o(inter)=1.
     356      enddo
     357      do inter=25,35
     358         efdish2o2(inter)=1.
     359      enddo
     360      do inter=34,36
     361         efdiso3(inter)=1.
     362      enddo
     363      do inter=27,36
     364         efdisno2(inter)=1.
     365      enddo
     366      do inter=1,15
     367         efdish2(inter)=1.
     368      enddo
    366369         
    367370      !4 possible channels for CO2 ionization
     371      open(130,file=trim(datafile)//'/EUVDAT'//
     372     $     '/co2ion_branchingratio_schunkandnagy2000_param.dat')
     373      do inter=1,16
     374         read(130,*)i,nada,efionco2(inter,1),efionco2(inter,2),
     375     $        efionco2(inter,3),efionco2(inter,4)
     376         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
     377         efdisco2(inter)=1.-nada
     378         efionco2(inter,1)=(1.-efdisco2(inter))*efionco2(inter,1)
     379         efionco2(inter,2)=(1.-efdisco2(inter))*efionco2(inter,2)
     380         efionco2(inter,3)=(1.-efdisco2(inter))*efionco2(inter,3)
     381         efionco2(inter,4)=(1.-efdisco2(inter))*efionco2(inter,4)
     382      enddo
     383      close(130)
     384
     385      do inter=17,36
     386         efdisco2(inter)=1.
     387      enddo
     388
    368389!      do inter=14,16
    369390!         efionco2(inter,1)=1.-efdisco2(inter)
    … …  
    380401!      enddo
    381402
     403      !2 possible channels for O2 ionization
     404      open(131,file=trim(datafile)//'/EUVDAT'//
     405     $     '/o2ion_branchingratio_schunkandnagy2000_param.dat')
     406      do inter=1,23
     407         read(131,*)i,nada,efiono2(inter,1),efiono2(inter,2)
     408         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
     409         efdiso2(inter)=1.-nada
     410         efiono2(inter,1)=(1.-efdiso2(inter))*efiono2(inter,1)
     411         efiono2(inter,2)=(1.-efdiso2(inter))*efiono2(inter,2)
     412      enddo
     413      close(131)
     414
     415      do inter=24,36
     416         efdiso2(inter)=1.
     417      enddo
     418
     419
    382420      !For O(3p) total ionization under 91.1 nm
    383 !      do inter=1,16
    384 !         efiono3p(inter)=1.d0
    385 !      enddo
     421      do inter=1,16
     422         efiono3p(inter)=1.d0
     423      enddo
    386424
    387425      !2 channels for N2 ionization
     426      open(132,file=trim(datafile)//'/EUVDAT'//
     427     $     '/n2ion_branchingratio_schunkandnagy2000_param.dat')
     428      do inter=1,15
     429         read(132,*)i,nada,efionn2(inter,1),efionn2(inter,2)
     430         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
     431         efdisn2(inter)=1.-nada
     432         efionn2(inter,1)=(1.-efdisn2(inter))*efionn2(inter,1)
     433         efionn2(inter,2)=(1.-efdisn2(inter))*efionn2(inter,2)
     434      enddo
     435      close(132)
     436     
     437      do inter=16,36
     438         efdisn2(inter)=1.
     439      enddo
     440
    388441!      do inter=9,15
    389442!         efionn2(inter,1)=1.-efdisn2(inter)
    … …  
    394447     
    395448      !3 channels for CO ionization
     449       open(133,file=trim(datafile)//'/EUVDAT'//
     450     $     '/coion_branchingratio_schunkandnagy2000_param.dat')
     451      do inter=1,16
     452         read(133,*)i,nada,efionco(inter,1),efionco(inter,2),
     453     $        efionco(inter,3)
     454         !Multiply the relative efficiency of each channel by the total ionization efficiency (second column)
     455         efdisco(inter)=1.-nada
     456         efionco(inter,1)=(1.-efdisco(inter))*efionco(inter,1)
     457         efionco(inter,2)=(1.-efdisco(inter))*efionco(inter,2)
     458         efionco(inter,3)=(1.-efdisco(inter))*efionco(inter,3)
     459      enddo
     460      close(133)
     461     
     462     
     463     
     464      do inter=17,36
     465         efdisco(inter)=1.
     466      enddo
     467     
    396468!      do inter=11,16
    397469!         efionco(inter,1)=1.-efdisco(inter)
    … …  
    408480!      enddo
    409481
     482
    410483      !Total ionization under 85 nm for N
    411 !      do inter=1,16
    412 !         efionn(inter)=1.
    413 !      enddo
     484      do inter=1,16
     485         efionn(inter)=1.
     486      enddo
    414487
    415488      !NO
    416 !      do inter=2,28
    417 !         efionno(inter)=1.-efdisno(inter)
    418 !      enddo
     489      do inter=2,28
     490         efionno(inter)=1.-efdisno(inter)
     491      enddo
    419492
    420493      !Total ionization under 90 nm for H
    421 !      do inter=3,16
    422 !         efionh(inter)=1.
    423 !      enddo
     494      do inter=3,16
     495         efionh(inter)=1.
     496      enddo
    424497
    425498

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

    r2795 r2836  
    1 subroutine 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)
     1!****************************************************************
     2!
     3!     Photochemical routine
     4!
     5!     Authors: Franck Lefevre, Francisco Gonzalez-Galindo
     6!     -------
     7!
     8!     Version: 14/11/2020
     9!
     10!     ASIS scheme : for details on the method see
     11!     Cariolle et al., Geosci. Model Dev., 10, 1467-1485, 2017.
     12!
     13!     -------
     14!
     15!     2022/09/15: adding ion chemistry by Antoine Martinez
     16!
     17!*****************************************************************
     18
     19subroutine photochemistry_venus(nz, n_lon, zlocal, ptimestep,                 &
     20                                    ok_jonline,ok_ionchem,tuneupperatm,       &
     21                                    nb_reaction_3_max,nb_reaction_4_max,      &
     22                                    nb_phot_max,nphotion,ig,                  &
     23                                    p, t, t_elect, tr, vmr_dens_euv, mumean,  &
     24                                    sza_input, lon, lat, nesp, nespeuv, iter, &
     25                                    prod_tr, loss_tr, em_no, em_o2) 
    226
    327use chemparam_mod
    428use photolysis_mod
     29use param_v4_h, only: jion
     30use iono_h, only: phdisrate
     31!nphot ! number of photodissociations, include in "use photolysis_mod" 
    532     
    633implicit none
    … …  
    1037!===================================================================
    1138
    12 integer, intent(in) :: nz         ! number of atmospheric layers
    13 integer, intent(in) :: nesp       ! number of tracers in traceur.def
    14 
    15 real, dimension(nz) :: p          ! pressure (hpa)
    16 real, dimension(nz) :: t          ! temperature (k)
    17 real, dimension(nz) :: zlocal     ! altitude (km)
    18 real, dimension(nz) :: mumean     ! mean molecular mass (g/mol)
    19 real :: ptimestep                 ! physics timestep (s)
    20 real :: sza_input                 ! solar zenith angle (degrees)
    21 real :: lon, lat                  ! longitude and latitude (degrees)
    22 
    23 integer :: n_lon                  ! for 1D test
     39integer, intent(in) :: nz          ! number of atmospheric layers
     40integer, intent(in) :: nesp        ! number of tracers in traceur.def
     41integer, intent(in) :: nespeuv     ! number of tracers for jthermcal_e107.F
     42
     43integer, intent(in) :: nb_reaction_3_max   ! number of quadratic reactions
     44integer, intent(in) :: nb_reaction_4_max   ! number of bimolecular reactions
     45integer, intent(in) :: nb_phot_max         ! total number of photolysis+photoionizations+quenching reactions
     46integer, intent(in) :: nphotion            ! number of photoionizations
     47
     48logical, intent(in) :: ok_ionchem  ! switch for ion chemistry
     49logical, intent(in) :: ok_jonline  ! switch for on-line calculation of photolysis rates
     50logical, intent(in) :: tuneupperatm! upper atmosphere chemistry tuning
     51
     52real, dimension(nz) :: p                 ! pressure (hpa)
     53real, dimension(nz) :: t                 ! temperature (k)
     54real, dimension(nz) :: t_elect           ! electronic temperature (k)
     55real, dimension(nz) :: zlocal            ! altitude (km)
     56real, dimension(nz) :: mumean            ! mean molecular mass (g/mol)
     57real, dimension(nz,nespeuv) :: vmr_dens_euv ! tracer mixing ratio for jthermcal_e107 routine
     58real :: ptimestep                        ! physics timestep (s)
     59real :: sza_input                        ! solar zenith angle (degrees)
     60real :: lon, lat                         ! longitude and latitude (degrees)
     61
     62integer :: ig                            ! grid point index
     63
     64integer :: n_lon                         ! for 1D test
    2465
    2566!===================================================================
    … …  
    3071real, dimension(nz,nesp) :: prod_tr ! production (cm-3.s-1)
    3172real, dimension(nz,nesp) :: loss_tr ! loss       (cm-3.s-1)
    32 real, dimension (nz) :: em_no       ! volume emission rate of no
    33 real, dimension (nz) :: em_o2       ! volume emission rate of o2(deltag)
     73real, dimension (nz)     :: em_no   ! volume emission rate of no
     74real, dimension (nz)     :: em_o2   ! volume emission rate of o2(deltag)
    3475
    3576!===================================================================
    … …  
    4384!===================================================================
    4485
    45 ! jonline
     86! ok_ jonline: see physiq.def
    4687! true : on-line calculation of photodissociation rates ! false : lookup table
    47 
    48 logical, save :: jonline = .true.
    4988
    5089logical, save :: firstcall = .true.
    … …  
    6099
    61100integer, parameter :: nj = 22, nztable = 281, nsza = 27, nso2 = 13
    62 real, dimension(nso2,nsza,nztable,nj), save :: jphot
     101!integer, parameter :: nztable = 281, nsza = 27, nso2 = 13
     102real, dimension(nso2,nsza,nztable,nj), save :: jphot ! nj must be equal to nphot
    63103real, dimension(nztable), save :: table_colair
    64104real, dimension(nso2,nztable), save :: table_colso2
    … …  
    68108! number densities
    69109
    70 real, dimension(nesp)    :: cold  ! number densities at previous timestep (molecule.cm-3)
    71 real, dimension(nz,nesp) :: c     ! number densities at current timestep (molecule.cm-3)
    72 real, dimension(nesp)    :: cnew  ! number densities at next timestep (molecule.cm-3)
    73      
     110real (kind = 8), dimension(nesp)        :: cold ! number densities at previous timestep (molecule.cm-3)
     111real (kind = 8), dimension(nz,nesp)     :: c    ! number densities at current timestep (molecule.cm-3)
     112real (kind = 8), dimension(nz,nespeuv)  :: c_euv! number densities for jthermcal_e107 at current timestep (molecule.cm-3)
     113real (kind = 8), dimension(nesp)        :: cnew ! number densities at next timestep (molecule.cm-3)
     114
    74115! timesteps
    75116
    … …  
    80121integer :: phychemrat
    81122
     123!Tracer indexes for photionization coeffs
     124
     125integer,parameter :: induv_co2 = 1
     126integer,parameter :: induv_o2  = 2
     127integer,parameter :: induv_o   = 3
     128integer,parameter :: induv_h2o = 4
     129integer,parameter :: induv_h2  = 5
     130integer,parameter :: induv_h2o2= 6
     131integer,parameter :: induv_o3  = 7
     132integer,parameter :: induv_n2  = 8
     133integer,parameter :: induv_n   = 9
     134integer,parameter :: induv_no  = 10
     135integer,parameter :: induv_co  = 11
     136integer,parameter :: induv_h   = 12
     137integer,parameter :: induv_no2 = 13
     138
    82139! reaction rates
    83140
    84 integer, parameter :: nb_phot_max       = 35
    85 integer, parameter :: nb_phot           = 22
    86 integer, parameter :: nb_reaction_3_max = 12
    87 integer, parameter :: nb_reaction_4_max = 98
    88 
    89 real, dimension(nz,nb_phot_max)      :: v_phot
    90 real, dimension(nz,nb_reaction_3_max) :: v_3
    91 real, dimension(nz,nb_reaction_4_max) :: v_4
     141!integer, parameter :: nb_phot_max       = 35
     142!integer, parameter :: nb_phot           = 22
     143!integer, parameter :: nb_reaction_3_max = 12
     144!integer, parameter :: nb_reaction_4_max = 98
     145
     146real (kind = 8), dimension(nz,      nb_phot_max) :: v_phot
     147real (kind = 8), dimension(nz,nb_reaction_3_max) :: v_3
     148real (kind = 8), dimension(nz,nb_reaction_4_max) :: v_4
    92149
    93150logical,parameter :: hetero_ice  = .false.
    … …  
    96153! matrix
    97154
    98 real, dimension(nesp,nesp) :: mat, mat1
    99 integer, dimension(nesp)   :: indx
    100 integer                    :: code
     155real (kind = 8), dimension(nesp,nesp) :: mat, mat1
     156integer, dimension(nesp)              :: indx
     157integer                               :: code
    101158
    102159! production and loss terms (for first-guess solution only)
    103160
    104 real, dimension(nesp) :: prod, loss, lossconc
     161real (kind = 8), dimension(nesp) :: prod, loss, lossconc
    105162
    106163! indexes
    … …  
    113170!===================================================================
    114171
    115    if (jonline) then
     172   ! ok jonline
     173   ! true : on-line calculation of photodissociation rates ! false : lookup table
     174   if (ok_jonline) then
    116175      print*, 'photochemistry: Read UV absorption cross-sections:'
    117176      call init_photolysis
    118177   else
    119178      print*, 'photochemistry: Read photolysis lookup table:'
    120       call init_chimie(nj, nztable, nsza, nso2, jphot, table_colair, table_colso2, table_sza)
     179      call init_chimie(nphot, nztable, nsza, nso2, jphot, table_colair, table_colso2, table_sza)
    121180   end if
    122181
    … …  
    125184!===================================================================
    126185
    127    call indice(nb_phot_max, nb_reaction_3_max, nb_reaction_4_max)
     186   call indice(ok_ionchem, nb_phot_max, nb_reaction_3_max, nb_reaction_4_max)
    128187
    129188   firstcall = .false.
    … …  
    140199
    141200do iz = 1,nz
    142    conc(iz) = p(iz)/(1.38E-19*t(iz))
    143    c(iz,:) = tr(iz,:)*conc(iz)
     201   conc(iz)    = p(iz)/(1.38E-19*t(iz))
     202   c(iz,:)     = tr(iz,:)*conc(iz)
    144203end do
    145204     
    … …  
    152211dist_sol = 0.72333
    153212
    154 if (jonline) then
    155    if (sza_input <= 95.) then ! day at 30 km
     213if (ok_jonline) then
     214   if (sza_input <= 95.) then ! day at 300 km
    156215      call photolysis_online(nz, nb_phot_max, zlocal, p,                 &
    157216                             t, mumean, i_co2,i_co, i_o, i_o1d,          &
    … …  
    161220                             i_no2, i_no, i_n2, i_n2d,                   &
    162221                             nesp, tr, sza_input, dist_sol, v_phot)
     222
     223      !Calculation of photoionization rates, if needed
     224      if (ok_ionchem) then
     225         do iz = 1,nz
     226           c_euv(iz,:) = vmr_dens_euv(iz,:)*conc(iz)
     227         end do
     228      !! FAIRE ULTRA ATTENTION LE c_euv et vmr_dens_euv NE SONT PAS DANS LE MEME REGIME
     229      !! DE i_ESPECE QUE C, CNEW, COLD et TR
     230         call jthermcalc_e107(ig,nz,2,c_euv,nespeuv,t,zlocal,sza_input)
     231         do iz=1,nz
     232            call phdisrate(ig,nz,2,sza_input,iz)
     233         end do
     234         !CO2 photoionization
     235         v_phot(:,nphot+ 1) = jion(induv_co2,:,1)
     236         v_phot(:,nphot+ 2) = jion(induv_co2,:,2)
     237         v_phot(:,nphot+ 3) = jion(induv_co2,:,2)
     238         v_phot(:,nphot+ 4) = jion(induv_co2,:,3)
     239         v_phot(:,nphot+ 5) = jion(induv_co2,:,3)
     240         v_phot(:,nphot+ 6) = jion(induv_co2,:,4)
     241         v_phot(:,nphot+ 7) = jion(induv_co2,:,4)
     242         !O2 photoionization
     243         v_phot(:,nphot+ 8) = jion(induv_o2,:,1)
     244         !O photoionization
     245         v_phot(:,nphot+ 9) = jion(induv_o,:,1)
     246         !NO photoionization
     247         v_phot(:,nphot+10) = jion(induv_no,:,1)
     248         !CO photoionization
     249         v_phot(:,nphot+11) = jion(induv_co,:,1)
     250         v_phot(:,nphot+12) = jion(induv_co,:,2)
     251         v_phot(:,nphot+13) = jion(induv_co,:,2)
     252         !N2 photoionization
     253         v_phot(:,nphot+14) = jion(induv_n2,:,1)
     254         v_phot(:,nphot+15) = jion(induv_n2,:,2)
     255         v_phot(:,nphot+16) = jion(induv_n2,:,2)
     256         !N photoionization
     257         v_phot(:,nphot+17) = jion(induv_n,:,1)
     258         !H photoionization
     259         v_phot(:,nphot+18) = jion(induv_h,:,1)
     260      end if
    163261   else ! night
    164262      v_phot(:,:) = 0.
    165263   end if
    166264else
    167    call phot(nj, nztable, nsza, nso2, sza_input, dist_sol, mumean, tr(:,i_co2), tr(:,i_so2),         &
     265   call phot(   nj, nztable, nsza, nso2, sza_input, dist_sol, mumean, tr(:,i_co2), tr(:,i_so2),         &
    168266             jphot, table_colair, table_colso2, table_sza, nz, nb_phot_max, t, p, v_phot)
    169267end if
    … …  
    173271!===================================================================
    174272                   
    175 call 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)
     273call krates(hetero_ice,hetero_dust,ok_ionchem, nphotion, nz, nesp, c, conc, t, t_elect, p, nb_phot_max, nb_reaction_3_max, &
     274            nb_reaction_4_max, tuneupperatm, v_3, v_4, v_phot, sza_input, ind_norec, ind_orec)
    177275
    178276!===================================================================
    … …  
    248346   c(iz,:) = cnew(:)
    249347   cnew(:) = 0.
     348
     349!  force charge neutrality (mod fgg, july 2019)
     350
     351   if (ok_ionchem) then
     352      if(c(iz,i_elec).ne.c(iz,i_co2plus)+c(iz,i_oplus)+c(iz,i_o2plus)+  &
     353           c(iz,i_noplus)+c(iz,i_coplus)+c(iz,i_cplus)+c(iz,i_n2plus)+  &
     354           c(iz,i_nplus)+c(iz,i_hplus)+c(iz,i_hco2plus)+                &
     355           c(iz,i_hcoplus)+c(iz,i_h2oplus)+c(iz,i_h3oplus)+             &
     356           c(iz,i_ohplus)) then
     357         c(iz,i_elec) = c(iz,i_co2plus)+c(iz,i_oplus)+c(iz,i_o2plus)+   &
     358              c(iz,i_noplus)+c(iz,i_coplus)+c(iz,i_cplus)+              &
     359              c(iz,i_n2plus)+c(iz,i_nplus)+c(iz,i_hplus)+               &
     360              c(iz,i_hco2plus)+c(iz,i_hcoplus)+c(iz,i_h2oplus)+         &
     361              c(iz,i_h3oplus)+c(iz,i_ohplus)
     362         !      write(*,*)'photochemistry/359'
     363         !      write(*,*)'Forcing charge neutrality at ilev,',ilev,' ig=',ig
     364      end if
     365   end if
    250366
    251367!  increment internal time
    … …  
    323439!======================================================================
    324440
    325  subroutine indice(nb_phot_max, nb_reaction_3_max, nb_reaction_4_max)
     441 subroutine indice(ok_ionchem, nb_phot_max, nb_reaction_3_max, nb_reaction_4_max)
    326442
    327443!================================================================
    … …  
    344460! input
    345461
    346 integer :: nb_phot_max, nb_reaction_3_max, nb_reaction_4_max
     462integer, intent(in) :: nb_reaction_3_max   ! number of quadratic reactions
     463integer, intent(in) :: nb_reaction_4_max   ! number of bimolecular reactions
     464integer, intent(in) :: nb_phot_max         ! number of processes treated numerically as photodissociations
     465logical, intent(in) :: ok_ionchem          ! True: Ion reaction
    347466
    348467! local
    … …  
    537656indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n, 1.0, i_n2d)
    538657
    539 !===========================================================
    540 !      a001 : O + O2 + CO2 -> O3 + CO2
     658!Only if ion chemistry included
     659if (ok_ionchem) then
     660
     661!===========================================================
     662!      CO2 + hv -> CO2+ + e-
     663!===========================================================
     664
     665   nb_phot = nb_phot + 1
     666
     667   indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co2plus, 1.0, i_elec)
     668
     669!===========================================================
     670!      CO2 + hv -> O+ + CO + e-
     671!===========================================================
     672!We divide this reaction in two
     673
     674!0.5 CO2 + hv -> CO
     675   nb_phot = nb_phot + 1
     676
     677   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_co, 0.0, i_dummy)
     678
     679!0.5 CO2 + hv -> O+ + e-
     680   nb_phot = nb_phot + 1
     681
     682   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_oplus, 1.0, i_elec)
     683
     684!===========================================================
     685!      CO2 + hv -> CO+ + O + e-
     686!===========================================================
     687!We divide this reaction in two
     688
     689!0.5 CO2 + hv -> O
     690   nb_phot = nb_phot + 1
     691
     692   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o, 0.0, i_dummy)
     693
     694!0.5 CO2 + hv -> CO+ + e-
     695   nb_phot = nb_phot + 1
     696
     697   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_coplus, 1.0, i_elec)
     698
     699!===========================================================
     700!      CO2 + hv -> C+ + O2 + e-
     701!===========================================================
     702!We divide this reaction in two
     703
     704!0.5 CO2 + hv -> O2
     705   nb_phot = nb_phot + 1
     706
     707   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o2, 0.0, i_dummy)
     708
     709!0.5 CO2 + hv -> C+ + e-
     710   nb_phot = nb_phot + 1
     711
     712   indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_cplus, 1.0, i_elec)
     713
     714!===========================================================
     715!      O2 + hv -> O2+ + e-
     716!===========================================================
     717
     718   nb_phot = nb_phot + 1
     719
     720   indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o2plus, 1.0, i_elec)
     721
     722!===========================================================
     723!      O + hv -> O+ + e-
     724!===========================================================
     725
     726   nb_phot = nb_phot + 1
     727
     728   indice_phot(nb_phot) = z3spec(1.0, i_o, 1.0, i_oplus, 1.0, i_elec)
     729
     730!===========================================================
     731!      NO + hv -> NO+ + e-
     732!===========================================================
     733
     734   nb_phot = nb_phot + 1
     735
     736   indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_noplus, 1.0, i_elec)
     737
     738!===========================================================
     739!      CO + hv -> CO+ + e-
     740!===========================================================
     741
     742   nb_phot = nb_phot + 1
     743
     744   indice_phot(nb_phot) = z3spec(1.0, i_co, 1.0, i_coplus, 1.0, i_elec)
     745
     746!===========================================================
     747!      CO + hv -> C+ + O + e-
     748!===========================================================
     749!We divide this reaction in two
     750
     751!0.5 CO + hv -> O
     752   nb_phot = nb_phot + 1
     753
     754   indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_o, 0.0, i_dummy)
     755
     756!0.5 CO + hv -> C+ + e-
     757   nb_phot = nb_phot + 1
     758
     759   indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_cplus, 1.0, i_elec)
     760
     761!===========================================================
     762!      N2 + hv -> N2+ + e-
     763!===========================================================
     764
     765   nb_phot = nb_phot + 1
     766
     767   indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2plus, 1.0, i_elec)
     768
     769!===========================================================
     770!      N2 + hv -> N+ + N + e-
     771!===========================================================
     772!We divide this reaction in two
     773
     774!0.5 N2 + hv -> N
     775   nb_phot = nb_phot + 1
     776
     777   indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_n, 0.0, i_dummy)
     778
     779!0.5 N2 + hv -> N+ + e-
     780   nb_phot = nb_phot + 1
     781
     782   indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_nplus, 1.0, i_elec)
     783
     784!===========================================================
     785!      N + hv -> N+ + e-
     786!===========================================================
     787
     788   nb_phot = nb_phot + 1
     789
     790   indice_phot(nb_phot) = z3spec(1.0, i_n, 1.0, i_nplus, 1.0, i_elec)
     791
     792!===========================================================
     793!      H + hv -> H+ + e-
     794!===========================================================
     795
     796   nb_phot = nb_phot + 1
     797
     798   indice_phot(nb_phot) = z3spec(1.0, i_h, 1.0, i_hplus, 1.0, i_elec)
     799
     800end if   !ok_ionchem
     801
     802!===========================================================
     803!      a001 : O + O2 + (CO2 or M) -> O3 + (CO2 or M)
    541804!===========================================================
    542805
    … …  
    546809
    547810!===========================================================
    548 !      a002 : O + O + CO2 -> O2(Dg) + CO2
     811!      a002 : O + O + (CO2 or M) -> O2(Dg) + (CO2 or M)
    549812!===========================================================
    550813
    … …  
    15021765indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2)
    15031766
    1504 !===========================================================
    1505 !      i001: O2(Dg) + CO2 -> O2 + CO2
     1767!Only if ion chemistry
     1768if (ok_ionchem) then
     1769
     1770!===========================================================
     1771!      i001 : CO2+ + O2 -> O2+ + CO2
     1772!===========================================================
     1773
     1774   nb_reaction_4 = nb_reaction_4 + 1
     1775
     1776   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_co2)
     1777
     1778!===========================================================
     1779!      i002 : CO2+ + O -> O+ + CO2
     1780!===========================================================
     1781
     1782   nb_reaction_4 = nb_reaction_4 + 1
     1783
     1784   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co2)
     1785
     1786!===========================================================
     1787!      i003 : CO2+ + O -> O2+ + CO
     1788!===========================================================
     1789
     1790   nb_reaction_4 = nb_reaction_4 + 1
     1791
     1792   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_co)
     1793
     1794!===========================================================
     1795!      i004 : O2+ + e- -> O + O
     1796!===========================================================
     1797
     1798   nb_reaction_4 = nb_reaction_4 + 1
     1799
     1800   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_elec, 2.0, i_o, 0.0, i_dummy)
     1801
     1802!===========================================================
     1803!      i005 : O+ + CO2 -> O2+ + CO
     1804!===========================================================
     1805
     1806   nb_reaction_4 = nb_reaction_4 + 1
     1807
     1808   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_co2, 1.0, i_o2plus, 1.0, i_co)
     1809
     1810!===========================================================
     1811!      i006 : CO2+ + e -> CO + O
     1812!===========================================================
     1813
     1814   nb_reaction_4 = nb_reaction_4 + 1
     1815
     1816   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_elec, 1.0, i_co, 1.0, i_o)
     1817
     1818!===========================================================
     1819!      i007 : CO2+ + NO -> NO+ + CO2
     1820!===========================================================
     1821
     1822   nb_reaction_4 = nb_reaction_4 + 1
     1823
     1824   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_co2)
     1825
     1826!===========================================================
     1827!      i008 : O2+ + NO -> NO+ + O2
     1828!===========================================================
     1829
     1830   nb_reaction_4 = nb_reaction_4 + 1
     1831
     1832   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_o2)
     1833
     1834!===========================================================
     1835!      i009 : O2+ + N2 -> NO+ + NO
     1836!===========================================================
     1837
     1838   nb_reaction_4 = nb_reaction_4 + 1
     1839
     1840   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_no)
     1841
     1842!===========================================================
     1843!      i010 : O2+ + N -> NO+ + O
     1844!===========================================================
     1845
     1846   nb_reaction_4 = nb_reaction_4 + 1
     1847
     1848   indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n, 1.0, i_noplus, 1.0, i_o)
     1849
     1850!===========================================================
     1851!      i011 : O+ + N2 -> NO+ + N
     1852!===========================================================
     1853
     1854   nb_reaction_4 = nb_reaction_4 + 1
     1855
     1856   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_n)
     1857
     1858!===========================================================
     1859!      i012 : NO+ + e -> N + O
     1860!===========================================================
     1861
     1862   nb_reaction_4 = nb_reaction_4 + 1
     1863
     1864   indice_4(nb_reaction_4) = z4spec(1.0, i_noplus, 1.0, i_elec, 1.0, i_n, 1.0, i_o)
     1865
     1866!===========================================================
     1867!      i013 : CO+ + CO2 -> CO2+ + CO
     1868!===========================================================
     1869
     1870   nb_reaction_4 = nb_reaction_4 + 1
     1871
     1872   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_co)
     1873
     1874!===========================================================
     1875!      i014 : CO+ + O -> O+ + CO
     1876!===========================================================
     1877
     1878   nb_reaction_4 = nb_reaction_4 + 1
     1879
     1880   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co)
     1881
     1882!===========================================================
     1883!      i015 : C+ + CO2 -> CO+ + CO
     1884!===========================================================
     1885
     1886   nb_reaction_4 = nb_reaction_4 + 1
     1887
     1888   indice_4(nb_reaction_4) = z4spec(1.0, i_cplus, 1.0, i_co2, 1.0, i_coplus, 1.0, i_co)
     1889
     1890!===========================================================
     1891!      i016 : N2+ + CO2 -> CO2+ + N2
     1892!===========================================================
     1893
     1894   nb_reaction_4 = nb_reaction_4 + 1
     1895
     1896   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n2)
     1897
     1898!===========================================================
     1899!      i017 : N2+ + O -> NO+ + N
     1900!===========================================================
     1901
     1902   nb_reaction_4 = nb_reaction_4 + 1
     1903
     1904   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_noplus, 1.0, i_n)
     1905
     1906!===========================================================
     1907!      i018 : N2+ + CO -> CO+ + N2
     1908!===========================================================
     1909
     1910   nb_reaction_4 = nb_reaction_4 + 1
     1911
     1912   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co, 1.0, i_coplus, 1.0, i_n2)
     1913
     1914!===========================================================
     1915!      i019 : N2+ + e -> N + N
     1916!===========================================================
     1917
     1918   nb_reaction_4 = nb_reaction_4 + 1
     1919
     1920   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_elec, 2.0, i_n, 0.0, i_dummy)
     1921
     1922!===========================================================
     1923!      i020 : N2+ + O -> O+ + N2
     1924!===========================================================
     1925
     1926   nb_reaction_4 = nb_reaction_4 + 1
     1927
     1928   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_n2)
     1929
     1930!===========================================================
     1931!      i021 : N+ + CO2 -> CO2+ + N
     1932!===========================================================
     1933
     1934   nb_reaction_4 = nb_reaction_4 + 1
     1935
     1936   indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n)
     1937
     1938!===========================================================
     1939!      i022 : CO+ + H -> H+ + CO
     1940!===========================================================
     1941
     1942   nb_reaction_4 = nb_reaction_4 + 1
     1943
     1944   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_co)
     1945
     1946!===========================================================
     1947!      i023 : O+ + H -> H+ + O
     1948!===========================================================
     1949
     1950   nb_reaction_4 = nb_reaction_4 + 1
     1951
     1952   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_o)
     1953
     1954!===========================================================
     1955!      i024 : H+ + O -> O+ + H
     1956!===========================================================
     1957
     1958   nb_reaction_4 = nb_reaction_4 + 1
     1959
     1960   indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_h)
     1961
     1962!===========================================================
     1963!      i025 : CO2+ + H2 -> HCO2+ + H
     1964!===========================================================
     1965
     1966   nb_reaction_4 = nb_reaction_4 + 1
     1967
     1968   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2, 1.0, i_hco2plus, 1.0, i_h)
     1969
     1970!===========================================================
     1971!      i026 : HCO2+ + e -> H + CO2
     1972!===========================================================
     1973
     1974   nb_reaction_4 = nb_reaction_4 + 1
     1975
     1976   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
     1977
     1978!===========================================================
     1979!      i027 : HCO2+ + e -> H + O + CO
     1980!===========================================================
     1981!We divide this reaction in two
     1982
     1983!0.5HCO2+ + 0.5e -> H
     1984
     1985   nb_reaction_4 = nb_reaction_4 + 1
     1986
     1987   indice_4(nb_reaction_4) = z4spec(.5, i_hco2plus, 0.5, i_elec, 1.0, i_h, 0.0, i_dummy)
     1988
     1989!0.5 HCO2+ + 0.5 e -> O + CO
     1990
     1991   nb_reaction_4 = nb_reaction_4 + 1
     1992
     1993   indice_4(nb_reaction_4) = z4spec(0.5, i_hco2plus, 0.5, i_elec, 1.0, i_o, 1.0, i_co)
     1994
     1995!===========================================================
     1996!      i029 : HCO2+ + e -> OH + CO
     1997!===========================================================
     1998
     1999   nb_reaction_4 = nb_reaction_4 + 1
     2000
     2001   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_oh, 1.0, i_co)
     2002
     2003
     2004!===========================================================
     2005!      i030 : HCO2+ + e -> H + CO2
     2006!===========================================================
     2007
     2008   nb_reaction_4 = nb_reaction_4 + 1
     2009
     2010   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2)
     2011
     2012
     2013!===========================================================
     2014!      i031 : HCO2+ + O -> HCO+ + O2
     2015!===========================================================
     2016
     2017   nb_reaction_4 = nb_reaction_4 + 1
     2018
     2019   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_o, 1.0, i_hcoplus, 1.0, i_o2)
     2020
     2021
     2022!===========================================================
     2023!      i032 : HCO2+ + CO -> HCO+ + CO2
     2024!===========================================================
     2025
     2026   nb_reaction_4 = nb_reaction_4 + 1
     2027   indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_co2)
     2028
     2029
     2030!===========================================================
     2031!      i033 : H+ + CO2 -> HCO+ + O
     2032!===========================================================
     2033
     2034   nb_reaction_4 = nb_reaction_4 + 1
     2035   indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_co2, 1.0, i_hcoplus, 1.0, i_o)
     2036
     2037
     2038!===========================================================
     2039!      i034 : CO2+ + H -> HCO+ + O
     2040!===========================================================
     2041
     2042   nb_reaction_4 = nb_reaction_4 + 1
     2043   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h, 1.0, i_hcoplus, 1.0, i_o)
     2044
     2045
     2046!===========================================================
     2047!      i035 : CO+ + H2 -> HCO+ + H
     2048!===========================================================
     2049
     2050   nb_reaction_4 = nb_reaction_4 + 1
     2051   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h2, 1.0, i_hcoplus, 1.0, i_h)
     2052
     2053
     2054!===========================================================
     2055!      i036 : HCO+ + e- -> CO + H
     2056!===========================================================
     2057
     2058   nb_reaction_4 = nb_reaction_4 + 1
     2059   indice_4(nb_reaction_4) = z4spec(1.0, i_hcoplus, 1.0, i_elec, 1.0, i_co, 1.0, i_h)
     2060
     2061!===========================================================
     2062!      i037 : CO2+ + H2O -> H2O+ + CO2
     2063!===========================================================
     2064
     2065   nb_reaction_4 = nb_reaction_4 + 1
     2066   indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_co2)
     2067
     2068!===========================================================
     2069!      i038 : CO+ + H2O -> H2O+ + CO
     2070!===========================================================
     2071
     2072   nb_reaction_4 = nb_reaction_4 + 1
     2073   indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_co)
     2074
     2075!===========================================================
     2076!      i039 : O+ + H2O -> H2O+ + O
     2077!===========================================================
     2078
     2079   nb_reaction_4 = nb_reaction_4 + 1
     2080   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_o)
     2081
     2082!===========================================================
     2083!      i040 : N2+ + H2O -> H2O+ + N2
     2084!===========================================================
     2085
     2086   nb_reaction_4 = nb_reaction_4 + 1
     2087   indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_n2)
     2088
     2089!===========================================================
     2090!      i041 : N+ + H2O -> H2O+ + N
     2091!===========================================================
     2092
     2093   nb_reaction_4 = nb_reaction_4 + 1
     2094   indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_n)
     2095
     2096!===========================================================
     2097!      i042 : H+ + H2O -> H2O+ + H
     2098!===========================================================
     2099
     2100   nb_reaction_4 = nb_reaction_4 + 1
     2101   indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_h)
     2102
     2103!===========================================================
     2104!      i043 : H2O+ + O2 -> O2+ + H2O
     2105!===========================================================
     2106
     2107   nb_reaction_4 = nb_reaction_4 + 1
     2108   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_h2o)
     2109
     2110!===========================================================
     2111!      i044 : H2O+ + CO -> HCO+ + OH
     2112!===========================================================
     2113
     2114   nb_reaction_4 = nb_reaction_4 + 1
     2115   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_oh)
     2116
     2117!===========================================================
     2118!      i045 : H2O+ + O -> O2+ + H2
     2119!===========================================================
     2120
     2121   nb_reaction_4 = nb_reaction_4 + 1
     2122   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_h2)
     2123
     2124!===========================================================
     2125!      i046 : H2O+ + NO -> NO+ + H2O
     2126!===========================================================
     2127
     2128   nb_reaction_4 = nb_reaction_4 + 1
     2129   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_no, 1.0, i_noplus, 1.0, i_h2o)
     2130
     2131!===========================================================
     2132!      i047 : H2O+ + e- -> H + H + O
     2133!===========================================================
     2134
     2135   nb_reaction_4 = nb_reaction_4 + 1
     2136   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 2.0, i_h, 1.0, i_o)
     2137
     2138!===========================================================
     2139!      i048 : H2O+ + e- -> H + OH
     2140!===========================================================
     2141
     2142   nb_reaction_4 = nb_reaction_4 + 1
     2143   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 1.0, i_h, 1.0, i_oh)
     2144
     2145!===========================================================
     2146!      i049 : H2O+ + e- -> H2 + O
     2147!===========================================================
     2148
     2149   nb_reaction_4 = nb_reaction_4 + 1
     2150   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 1.0, i_h2, 1.0, i_o)
     2151
     2152!===========================================================
     2153!      i050 : H2O+ + H2O -> H3O+ + OH
     2154!===========================================================
     2155
     2156   nb_reaction_4 = nb_reaction_4 + 1
     2157   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_oh)
     2158
     2159!===========================================================
     2160!      i051 : H2O+ + H2 -> H3O+ + H
     2161!===========================================================
     2162
     2163   nb_reaction_4 = nb_reaction_4 + 1
     2164   indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2, 1.0, i_h3oplus, 1.0, i_h)
     2165
     2166!===========================================================
     2167!      i052 : HCO+ + H2O -> H3O+ + CO
     2168!===========================================================
     2169
     2170   nb_reaction_4 = nb_reaction_4 + 1
     2171   indice_4(nb_reaction_4) = z4spec(1.0, i_hcoplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_co)
     2172
     2173!===========================================================
     2174!      i053: H3O+ + e -> OH + H + H
     2175!===========================================================
     2176
     2177   nb_reaction_4 = nb_reaction_4 + 1
     2178   indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_oh, 2.0, i_h)
     2179
     2180!===========================================================
     2181!      i054: H3O+ + e -> H2O + H
     2182!===========================================================
     2183
     2184   nb_reaction_4 = nb_reaction_4 + 1
     2185   indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_h2o, 1.0, i_h)
     2186
     2187!===========================================================
     2188!      i055: H3O+ + e -> HO + H2
     2189!===========================================================
     2190
     2191   nb_reaction_4 = nb_reaction_4 + 1
     2192   indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_oh, 1.0, i_h2)
     2193
     2194!===========================================================
     2195!      i056: H3O+ + e -> O + H2 + H
     2196!===========================================================
     2197!We divide this reaction in two
     2198
     2199!0.5H3O+ + 0.5e -> O
     2200
     2201   nb_reaction_4 = nb_reaction_4 + 1
     2202   indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_o, 0.0, i_dummy)
     2203
     2204!0.5H3O+ + 0.5e -> H2 + H
     2205
     2206   nb_reaction_4 = nb_reaction_4 + 1
     2207   indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_h2, 1.0, i_h)
     2208
     2209!===========================================================
     2210!      i057: O+ + H2 -> OH+ + H
     2211!===========================================================
     2212
     2213   nb_reaction_4 = nb_reaction_4 + 1
     2214   indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h2, 1.0, i_ohplus, 1.0, i_h)
     2215
     2216!===========================================================
     2217!      i058: OH+ + O -> O2+ + H
     2218!===========================================================
     2219
     2220   nb_reaction_4 = nb_reaction_4 + 1
     2221   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_h)
     2222
     2223!===========================================================
     2224!      i059: OH+ + CO2 -> HCO2+ + O
     2225!===========================================================
     2226
     2227   nb_reaction_4 = nb_reaction_4 + 1
     2228   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_co2, 1.0, i_hco2plus, 1.0, i_o)
     2229
     2230!===========================================================
     2231!      i060: OH+ + CO -> HCO+ + O
     2232!===========================================================
     2233
     2234   nb_reaction_4 = nb_reaction_4 + 1
     2235   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_o)
     2236
     2237!===========================================================
     2238!      i061: OH+ + NO -> NO+ + OH
     2239!===========================================================
     2240
     2241   nb_reaction_4 = nb_reaction_4 + 1
     2242   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_no, 1.0, i_noplus, 1.0, i_oh)
     2243
     2244!===========================================================
     2245!      i062: OH+ + H2 -> H2O+ + H
     2246!===========================================================
     2247
     2248   nb_reaction_4 = nb_reaction_4 + 1
     2249   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_h2, 1.0, i_h2oplus, 1.0, i_h)
     2250
     2251!===========================================================
     2252!      i063: OH+ + O2 -> O2+ + OH
     2253!===========================================================
     2254
     2255   nb_reaction_4 = nb_reaction_4 + 1
     2256   indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_oh)
     2257
     2258end if    !ok_ionchem
     2259
     2260!===========================================================
     2261!      j001: O2(Dg) + (CO2 or O) -> O2 + (CO2 or O)
    15062262!===========================================================
    15072263
    … …  
    15112267
    15122268!===========================================================
    1513 !      i002: O2(Dg) -> O2 + hv
     2269!      j002: O2(Dg) -> O2 + hv
    15142270!===========================================================
    15152271
    … …  
    17322488!-- TuneE
    17332489! VCD 2.0 tuning
    1734     v_phot(65:90,4) = v_phot(65:90,4)*10.
     2490    v_phot(65:90,4) = v_phot(65:90,4)*10. ! CO2 + hv ==> O(1D) + CO
    17352491!--
    17362492!   v_phot(:,4) = v_phot(:,4)*10.
    … …  
    17532509!======================================================================
    17542510
    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)
     2511 subroutine krates(hetero_ice,hetero_dust,ok_ionchem, nphotion,        &
     2512                   nz, nesp, c, conc, t, t_elect, p,                   &
     2513                   nb_phot_max, nb_reaction_3_max, nb_reaction_4_max,  &
     2514                   tuneupperatm, v_3, v_4, v_phot,sza_input,           &
     2515                   ind_norec, ind_orec)
    17562516 
    17572517!================================================================
    … …  
    17682528
    17692529 USE chemparam_mod
     2530 USE photolysis_mod, only : nphot
    17702531implicit none
    17712532
    1772 real, INTENT(IN), dimension(nz)  :: t, p, conc
    1773 integer, INTENT(IN) :: nesp, nj, nz
    1774 real, INTENT(IN) :: sza_input
     2533!----------------------------------------------------------------------
     2534!     input
     2535!----------------------------------------------------------------------
     2536
     2537integer, INTENT(IN) :: nesp, nz
     2538real, INTENT(IN)    :: sza_input  ! [degree]
     2539real, INTENT(IN), dimension(nz)  :: t_elect, t, p, conc
     2540logical, INTENT(IN) :: hetero_ice, hetero_dust, tuneupperatm
     2541
     2542real, INTENT(IN), dimension(nz,nesp) :: c
     2543
     2544integer, intent(in) :: nb_reaction_3_max ! number of quadratic reactions
     2545integer, intent(in) :: nb_reaction_4_max ! number of bimolecular reactions
     2546integer, intent(in) :: nb_phot_max       ! number of reactions treated numerically as photodissociations
     2547integer, intent(in) :: nphotion          ! number of photoionizations
     2548logical, intent(in) :: ok_ionchem        ! if .true. then ionchem reaction used
     2549
     2550!----------------------------------------------------------------------
     2551!     output
     2552!----------------------------------------------------------------------
     2553
     2554real (kind = 8), dimension(nz,       nb_phot_max) :: v_phot
     2555real (kind = 8), dimension(nz, nb_reaction_3_max) :: v_3
     2556real (kind = 8), dimension(nz, nb_reaction_4_max) :: v_4
     2557integer :: ind_norec
     2558integer :: ind_orec
     2559
     2560!----------------------------------------------------------------------
     2561!     local
     2562!----------------------------------------------------------------------
     2563
    17752564integer :: iz
    1776 logical, INTENT(IN) :: hetero_ice, hetero_dust
    1777 integer :: ind_norec, ind_orec
    17782565real    :: ak0, ak1, xpo, rate, rate1, rate2, pi, gam
    17792566real    :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y
     2567integer :: nb_phot, nb_reaction_3, nb_reaction_4
    17802568real, dimension(nz)  :: surfice1d, surfdust1d
    17812569real, dimension(nz) :: deq
    … …  
    18072595                       g029, g030, g031, g032, g033,               &
    18082596                       h001, h002, h003,                           &
    1809                        i001, i002
    1810 
    1811 real, INTENT(IN), dimension(nz,nesp) :: c
    1812 
    1813 integer :: nb_phot_max, nb_reaction_3_max, nb_reaction_4_max, nb_reaction_3, nb_reaction_4, nb_phot
    1814 
    1815 real, dimension(nz,nb_phot_max) :: v_phot
    1816 real, dimension(nz,nb_reaction_3_max) :: v_3
    1817 real, dimension(nz,nb_reaction_4_max) :: v_4
    1818 
    1819 pi = acos(-1.)
    1820 
    1821 nb_phot       = nj
    1822 nb_reaction_3 = 0
    1823 nb_reaction_4 = 0
     2597                       i001, i002, i003, i004, i005, i006,         &
     2598                       i007, i008, i009, i010, i011, i012,         &
     2599                       i013, i014, i015, i016, i017, i018, i019,   &
     2600                       i020, i021, i022, i023, i024, i025, i026,   &
     2601                       i027, i028, i029, i030, i031, i032, i033,   &
     2602                       i034, i035, i036, i037, i038, i039, i040,   &
     2603                       i041, i042, i043, i044, i045, i046, i047,   &
     2604                       i048, i049, i050, i051, i052, i053, i054,   &
     2605                       i055, i056, i057, i058, i059, i060, i061,   &
     2606                       i062, i063,                                 &
     2607                       j001, j002
     2608
     2609!----------------------------------------------------------------------
     2610!     initialisation
     2611!----------------------------------------------------------------------
     2612
     2613      pi = acos(-1.)
     2614     
     2615      nb_phot       = nphot + nphotion ! initialised to the number of photolysis + number of photoionization rates
     2616      nb_reaction_3 = 0
     2617      nb_reaction_4 = 0
    18242618
    18252619!----------------------------------------------------------------------
    … …  
    18272621!----------------------------------------------------------------------
    18282622
    1829 !---  a001: o + o2 + co2 -> o3 + co2
     2623!---  a001: o + o2 + (co2 or M) -> o3 + (co2 or M)
    18302624
    18312625!     jpl 2003
    1832 
    1833       a001(:) = 2.5*6.0E-34*(t(:)/300.)**(-2.4)*conc(:)
     2626!     ------ BEFORE VCD 2.1 -------
     2627!      a001(:) = 2.5*6.0E-34*(t(:)/300.)**(-2.4)*conc(:)
     2628!     -----------------------------
     2629
     2630!    ! nist expression + take account of the N2 and CO2 factor
     2631!     ------ VCD 2.1 update -------
     2632      a001(:) = 6.0E-34*(t(:)/300.)**(-2.4) *    &
     2633              ( 2.5 * c(:,i_co2) +               &
     2634                1.0 * (conc(:)-c(:,i_co2))  )
     2635!     -----------------------------
    18342636
    18352637      nb_reaction_4 = nb_reaction_4 + 1
    18362638      v_4(:,nb_reaction_4) = a001(:)
    18372639
    1838 !---  a002: o + o + co2 -> o2 + co2
     2640!---  a002: o + o + (co2 or M) -> o2(Dg) + (co2 or M)
    18392641
    18402642!     Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986
    … …  
    18432645
    18442646!     Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973
    1845 
    1846       a002(:) = 2.5*9.46E-34*exp(485./t(:))*conc(:) ! nist expression
    1847 
     2647!     ------ BEFORE VCD 2.1 -------
     2648!      a002(:) = 2.5*9.46E-34*exp(485./t(:))*conc(:) ! nist expression
     2649!     -----------------------------
     2650
     2651!     ! nist expression + take account of the N2 and CO2 factor
     2652!       a002(:) = 9.46E-34*exp(485./t(:)) *     &
     2653!              ( 2.5 * c(:,i_co2) +            &
     2654!                1.0 * (conc(:)-c(:,i_co2)) )
     2655
     2656!     Baulch et al., 1976 (confirmed by smith and robertson, 2008)
     2657!     ------ VCD 2.1 update -------
     2658
     2659       a002(:) = 2.76E-34*exp(720./t(:)) *    &
     2660              ( 2.5 * c(:,i_co2) +            &
     2661                1.0 * (conc(:)-c(:,i_co2)) )
     2662!     -----------------------------
     2663     
    18482664      nb_reaction_3 = nb_reaction_3 + 1
    18492665      v_3(:,nb_reaction_3) = a002(:)
    … …  
    31263942
    31273943!----------------------------------------------------------------------
     3944!     ionospheric reactions
     3945!     only if ok_ionchem=true
     3946!----------------------------------------------------------------------
     3947
     3948      if (ok_ionchem) then
     3949
     3950!---     i001: co2+ + o2 -> o2+ + co2
     3951
     3952!        aninich, j. phys. chem. ref. data 1993
     3953
     3954         i001(:) = 5.5e-11*(300./t_elect(:))**0.5
     3955
     3956         nb_reaction_4 = nb_reaction_4 + 1
     3957         v_4(:,nb_reaction_4) = i001(:)
     3958
     3959!---     i002: co2+ + o -> o+ + co2
     3960
     3961!        UMIST database
     3962
     3963         i002(:) = 9.6e-11
     3964     
     3965         nb_reaction_4 = nb_reaction_4 + 1
     3966         v_4(:,nb_reaction_4) = i002(:)
     3967
     3968!---     i003: co2+ + o -> o2+ + co
     3969
     3970!        UMIST database
     3971
     3972         i003(:) = 1.64e-10
     3973
     3974         nb_reaction_4 = nb_reaction_4 + 1
     3975         v_4(:,nb_reaction_4) = i003(:)
     3976
     3977!---     i004: o2+ + e- -> o + o
     3978
     3979!        Alge et al., J. Phys. B, At. Mol. Phys. 1983
     3980         !i004(:) = 2.0e-7*(300./t_elect(:))**0.7
     3981         
     3982         do iz = 1,nz
     3983           if (t_elect(iz)<1200.) then
     3984             i004(iz) = 2.0e-7*(300./t_elect(iz))**0.7
     3985           else
     3986             i004(iz) = 7.4e-8*(1200./t_elect(iz))**0.56
     3987           end if
     3988         end do
     3989
     3990         nb_reaction_4 = nb_reaction_4 + 1
     3991         v_4(:,nb_reaction_4) = i004(:)
     3992
     3993!---     i005: o+ + co2 -> o2+ + co
     3994
     3995!        UMIST database
     3996
     3997         i005(:) = 9.4e-10
     3998
     3999         nb_reaction_4 = nb_reaction_4 + 1
     4000         v_4(:,nb_reaction_4) = i005(:)
     4001
     4002
     4003!---     i006: co2+ + e- -> co + o
     4004
     4005!        UMIST database
     4006
     4007         i006(:) = 3.8e-7*(300./t_elect(:))**0.5
     4008
     4009         nb_reaction_4 = nb_reaction_4 + 1
     4010         v_4(:,nb_reaction_4) = i006(:)
     4011
     4012
     4013!---     i007: co2+ + no -> no+ + co2
     4014
     4015!        UMIST database
     4016
     4017         i007(:) = 1.2e-10
     4018
     4019         nb_reaction_4 = nb_reaction_4 + 1
     4020         v_4(:,nb_reaction_4) = i007(:)
     4021
     4022!---     i008: o2+ + no -> no+ + o2
     4023
     4024!        UMIST database
     4025
     4026         i008(:) = 4.6e-10
     4027
     4028         nb_reaction_4 = nb_reaction_4 + 1
     4029         v_4(:,nb_reaction_4) = i008(:)
     4030
     4031!---     i009: o2+ + n2 -> no+ + no
     4032     
     4033!        Fox & Sung 2001
     4034
     4035         i009(:) = 1.0e-15
     4036     
     4037         nb_reaction_4 = nb_reaction_4 + 1
     4038         v_4(:,nb_reaction_4) = i009(:)
     4039
     4040!---     i010: o2+ + n -> no+ + o
     4041
     4042!        Fox & Sung 2001
     4043
     4044         i010(:) = 1.0e-10
     4045
     4046         nb_reaction_4 = nb_reaction_4 + 1
     4047         v_4(:,nb_reaction_4) = i010(:)
     4048
     4049!---     i011: o+ + n2 -> no+ + n
     4050
     4051!        Fox & Sung 2001
     4052
     4053         i011(:) = 1.2e-12 * (300./t_elect(:))**0.45
     4054
     4055         nb_reaction_4 = nb_reaction_4 + 1
     4056         v_4(:,nb_reaction_4) = i011(:)
     4057
     4058!---     i012: no+ + e -> n + o
     4059
     4060!        UMIST database
     4061
     4062         i012(:) = 4.3e-7*(300./t_elect(:))**0.37
     4063
     4064         nb_reaction_4 = nb_reaction_4 + 1
     4065         v_4(:,nb_reaction_4) = i012(:)
     4066
     4067
     4068!---     i013: co+ + co2 -> co2+ + co
     4069
     4070!        UMIST database
     4071
     4072         i013(:) = 1.0e-9
     4073
     4074         nb_reaction_4 = nb_reaction_4 + 1
     4075         v_4(:,nb_reaction_4) = i013(:)
     4076
     4077
     4078!---     i014: co+ + o -> o+ + co
     4079
     4080!        UMIST database
     4081
     4082         i014(:) = 1.4e-10
     4083
     4084         nb_reaction_4 = nb_reaction_4 + 1
     4085         v_4(:,nb_reaction_4) = i014(:)
     4086
     4087!---     i015: c+ + co2 -> co+ + co
     4088
     4089!        UMIST database
     4090
     4091         i015(:) = 1.1e-9
     4092
     4093         nb_reaction_4 = nb_reaction_4 + 1
     4094         v_4(:,nb_reaction_4) = i015(:)
     4095
     4096
     4097!---     i016: N2+ + co2 -> co2+ + N2
     4098
     4099!        Fox & Song 2001
     4100
     4101         i016(:) = 9.0e-10*(300./t_elect(:))**0.23
     4102
     4103         nb_reaction_4 = nb_reaction_4 + 1
     4104         v_4(:,nb_reaction_4) = i016(:)
     4105
     4106
     4107!---     i017: N2+ + o -> no+ + N
     4108
     4109!        Fox & Song 2001
     4110
     4111         i017(:) = 1.33e-10*(300./t_elect(:))**0.44
     4112
     4113         nb_reaction_4 = nb_reaction_4 + 1
     4114         v_4(:,nb_reaction_4) = i017(:)
     4115
     4116!---     i018: N2+ + co -> co+ + N2
     4117
     4118!        UMIST
     4119
     4120         i018(:) = 7.4e-11
     4121
     4122         nb_reaction_4 = nb_reaction_4 + 1
     4123         v_4(:,nb_reaction_4) = i018(:)
     4124
     4125!---     i019: N2+ + e -> N + N
     4126
     4127!        UMIST
     4128
     4129         i019(:) = 7.7e-7*(300./t_elect(:))**0.3
     4130
     4131         nb_reaction_4 = nb_reaction_4 + 1
     4132         v_4(:,nb_reaction_4) = i016(:)
     4133
     4134!---     i020: N2+ + o -> o+ + N2
     4135
     4136!        Fox & Song 2001
     4137
     4138         i020(:) = 7.0e-12*(300./t_elect(:))**0.23
     4139
     4140         nb_reaction_4 = nb_reaction_4 + 1
     4141         v_4(:,nb_reaction_4) = i020(:)
     4142
     4143!---     i021: N+ + co2 -> co2+ + N
     4144
     4145!        UMIST
     4146
     4147         i021(:) = 7.5e-10
     4148
     4149         nb_reaction_4 = nb_reaction_4 + 1
     4150         v_4(:,nb_reaction_4) = i021(:)
     4151
     4152!---     i022: CO+ + H -> H+ + CO
     4153
     4154!        Fox & Sung 2001
     4155
     4156         i022(:) = 4.0e-10
     4157
     4158         nb_reaction_4 = nb_reaction_4 + 1
     4159         v_4(:,nb_reaction_4) = i022(:)
     4160
     4161!---     i023: O+ + H -> H+ + O
     4162
     4163!        UMIST
     4164
     4165         i023(:) = 5.66e-10*((t_elect(:)/300.)**0.36)*exp(8.6/t_elect(:))
     4166
     4167         nb_reaction_4 = nb_reaction_4 + 1
     4168         v_4(:,nb_reaction_4) = i023(:)
     4169
     4170!---     i024: H+ + O -> O+ + H
     4171
     4172!        UMIST
     4173
     4174         i024(:) = 6.86e-10*((t_elect(:)/300.)**0.26)*exp(-224.3/t_elect(:))
     4175
     4176         nb_reaction_4 = nb_reaction_4 + 1
     4177         v_4(:,nb_reaction_4) = i024(:)
     4178
     4179!---     i025: CO+ + H2 -> HCO2+ + H
     4180
     4181!        UMIST
     4182
     4183         i025(:) = 9.5e-10
     4184
     4185         nb_reaction_4 = nb_reaction_4 + 1
     4186         v_4(:,nb_reaction_4) = i025(:)
     4187
     4188!---     i026: HCO2+ + e -> H + CO2
     4189
     4190!        UMIST
     4191
     4192         i026(:) = 1.75e-8*((300./t_elect(:))**0.5)
     4193
     4194         nb_reaction_4 = nb_reaction_4 + 1
     4195         v_4(:,nb_reaction_4) = i026(:)
     4196
     4197!---     i027+i028: HCO2+ + e -> H + O + CO
     4198
     4199!        UMIST
     4200         !Reaction splitted in 2: i027: 0.5 (HCO2+ + e-) -> H
     4201         !i028: 0.5 (HCO2+ + e-) -> O + CO
     4202
     4203         i027(:) = 8.1e-7*((300./t_elect(:))**0.64)
     4204
     4205         nb_reaction_4 = nb_reaction_4 + 1
     4206         v_4(:,nb_reaction_4) = i027(:)
     4207
     4208         nb_reaction_4 = nb_reaction_4 + 1
     4209         v_4(:,nb_reaction_4) = i027(:)
     4210
     4211!---     i029: HCO2+ + e -> OH + CO
     4212
     4213!        UMIST
     4214
     4215         i029(:) = 3.2e-7*((300./t_elect(:))**0.64)
     4216
     4217         nb_reaction_4 = nb_reaction_4 + 1
     4218         v_4(:,nb_reaction_4) = i029(:)
     4219
     4220!---     i030: HCO2+ + e -> H + CO2
     4221
     4222         i030(:) = 6.0e-8*((300./t_elect(:))**0.64)
     4223         nb_reaction_4 = nb_reaction_4 + 1
     4224         v_4(:,nb_reaction_4) = i030(:)
     4225
     4226!---     i031: HCO2+ + O -> HCO+ + O2
     4227
     4228!        UMIST
     4229
     4230         i031(:) = 1.e-9
     4231         nb_reaction_4 = nb_reaction_4 + 1
     4232         v_4(:,nb_reaction_4) = i031(:)
     4233
     4234!---     i032: HCO2+ + CO -> HCO+ + CO2
     4235
     4236!        UMIST, from Prassad & Huntress 1980
     4237
     4238         i032(:) = 7.8e-10
     4239         nb_reaction_4 = nb_reaction_4 + 1
     4240         v_4(:,nb_reaction_4) = i032(:)
     4241
     4242!---     i033: H+ + CO2 -> HCO+ + O
     4243
     4244!        UMIST, from Smith et al., Int. J. Mass Spectrom. Ion Proc., 117, 457-473(1992)
     4245
     4246         i033(:) = 3.5e-9
     4247         nb_reaction_4 = nb_reaction_4 + 1
     4248         v_4(:,nb_reaction_4) = i033(:)
     4249
     4250
     4251!---     i034: CO2+ + H -> HCO+ + O
     4252
     4253!        Seen in Fox 2015, from Borodi et al., Int. J. Mass Spectrom. 280, 218-225, 2009
     4254
     4255         i034(:) = 4.5e-10
     4256         nb_reaction_4 = nb_reaction_4 + 1
     4257         v_4(:,nb_reaction_4) = i034(:)
     4258
     4259!---     i035: CO+ + H2 -> HCO+ + H
     4260
     4261         !UMIST, from Scott et al., J. Chem. Phys., 106, 3982-3987(1997)
     4262
     4263         i035(:) = 7.5e-10
     4264         nb_reaction_4 = nb_reaction_4 + 1
     4265         v_4(:,nb_reaction_4) = i035(:)
     4266
     4267!---     i036: HCO+ + e- -> CO + H
     4268
     4269         !UMIST, from Mitchell, Phys. Rep., 186, 215 (1990)
     4270
     4271         i036(:) = 2.4e-7 *((300./t_elect(:))**0.69)
     4272         nb_reaction_4 = nb_reaction_4 + 1
     4273         v_4(:,nb_reaction_4) = i036(:)
     4274
     4275!---     i037: CO2+ + H2O -> H2O+ + CO2
     4276
     4277         !UMIST, from Karpas, Z., Anicich, V.G., and Huntress, W.T., Chem. Phys. Lett., 59, 84 (1978)
     4278
     4279         i037(:) = 2.04e-9 *((300./t_elect(:))**0.5)
     4280         nb_reaction_4 = nb_reaction_4 + 1
     4281         v_4(:,nb_reaction_4) = i037(:)
     4282
     4283!---     i038: CO+ + H2O -> H2O+ + CO
     4284
     4285         !UMIST, from Huntress, W.T., McEwan, M.J., Karpas, Z., and Anicich, V.G., Astrophys. J. Supp. Series, 44, 481 (1980)
     4286
     4287         i038(:) = 1.72e-9*((300./t_elect(:))**0.5)
     4288         nb_reaction_4 = nb_reaction_4 + 1
     4289         v_4(:,nb_reaction_4) = i038(:)
     4290
     4291!---     i039: O+ + H2O -> H2O+ + O
     4292
     4293         !UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978)
     4294
     4295         i039(:) = 3.2e-9*((300./t_elect(:))**0.5)
     4296         nb_reaction_4 = nb_reaction_4 + 1
     4297         v_4(:,nb_reaction_4) = i039(:)
     4298
     4299!---     i040: N2+ + H2O -> H2O+ + N2
     4300
     4301         !UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978)
     4302
     4303         i040(:) = 2.3e-9*((300./t_elect(:))**0.5)
     4304         nb_reaction_4 = nb_reaction_4 + 1
     4305         v_4(:,nb_reaction_4) = i040(:)
     4306
     4307!---     i041: N+ + H2O -> H2O+ + N
     4308
     4309         !UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978)
     4310
     4311         i041(:) = 2.8e-9*((300./t_elect(:))**0.5)
     4312         nb_reaction_4 = nb_reaction_4 + 1
     4313         v_4(:,nb_reaction_4) = i041(:)
     4314
     4315
     4316!---     i042: H+ + H2O -> H2O+ + H
     4317
     4318         !UMIST, from D. Smith, P. Spanel and C. A. Mayhew, Int. J. Mass Spectrom. Ion Proc., 117, 457-473(1992)
     4319
     4320         i042(:) = 6.9e-9*((300./t_elect(:))**0.5)
     4321         nb_reaction_4 = nb_reaction_4 + 1
     4322         v_4(:,nb_reaction_4) = i042(:)
     4323
     4324!---     i043: H2O+ + O2 -> O2+ + H2O
     4325
     4326         !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
     4327
     4328         i043(:) = 4.6e-10
     4329         nb_reaction_4 = nb_reaction_4 + 1
     4330         v_4(:,nb_reaction_4) = i043(:)
     4331
     4332!---     i044: H2O+ + CO -> HCO+ + OH
     4333
     4334         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     4335
     4336         i044(:) = 5.0e-10
     4337         nb_reaction_4 = nb_reaction_4 + 1
     4338         v_4(:,nb_reaction_4) = i044(:)
     4339
     4340!---     i045: H2O+ + O -> O2+ + H2
     4341
     4342         !UMIST, from Viggiano, A.A, Howarka, F., Albritton, D.L., Fehsenfeld, F.C., Adams, N.G., and Smith, D., Astrophys. J., 236, 492 (1980)
     4343
     4344         i045(:) = 4.0e-11
     4345         nb_reaction_4 = nb_reaction_4 + 1
     4346         v_4(:,nb_reaction_4) = i045(:)
     4347
     4348!---     i046: H2O+ + NO -> NO+ + H2O
     4349
     4350         !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
     4351
     4352         i046(:) = 2.7e-10
     4353         nb_reaction_4 = nb_reaction_4 + 1
     4354         v_4(:,nb_reaction_4) = i046(:)
     4355
     4356!---     i047: H2O+ + e- -> H + H + O
     4357
     4358         !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
     4359         
     4360         i047(:) = 3.05e-7*((300./t_elect(:))**0.5)
     4361         nb_reaction_4 = nb_reaction_4 + 1
     4362         v_4(:,nb_reaction_4) = i047(:)
     4363
     4364!---     i048: H2O+ + e- -> H + OH
     4365         
     4366         !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
     4367
     4368         i048(:) = 8.6e-8*((300./t_elect(:))**0.5)
     4369         nb_reaction_4 = nb_reaction_4 + 1
     4370         v_4(:,nb_reaction_4) = i048(:)
     4371
     4372!---     i049: H2O+ + e- -> O + H2
     4373
     4374         !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295
     4375
     4376         i049(:) = 3.9e-8*((300./t_elect(:))**0.5)
     4377         nb_reaction_4 = nb_reaction_4 + 1
     4378         v_4(:,nb_reaction_4) = i049(:)
     4379
     4380!---     i050: H2O+ + H2O -> H3O+ + OH
     4381
     4382         !UMIST, from Huntress, W.T. and Pinizzotto, R.F., J. Chem. Phys., 59, 4742 (1973)
     4383
     4384         i050(:) = 2.1e-9*((300./t_elect(:))**0.5)
     4385         nb_reaction_4 = nb_reaction_4 + 1
     4386         v_4(:,nb_reaction_4) = i050(:)
     4387
     4388
     4389!---     i051: H2O+ + H2 -> H3O+ + H
     4390
     4391         !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980)
     4392
     4393         i051(:) = 6.4e-10
     4394         nb_reaction_4 = nb_reaction_4 + 1
     4395         v_4(:,nb_reaction_4) = i051(:)
     4396
     4397!---     i052: HCO+ + H2O -> H3O+ + CO
     4398
     4399         !UMIST, from Adams, N.G., Smith, D., and Grief, D., Int. J. Mass Spectrom. Ion Phys., 26, 405 (1978)
     4400
     4401         i052(:) = 2.5e-9*((300./t_elect(:))**0.5)
     4402         nb_reaction_4 = nb_reaction_4 + 1
     4403         v_4(:,nb_reaction_4) = i052(:)
     4404
     4405!---     i053: H3O+ + e -> OH + H + H
     4406
     4407         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
     4408
     4409         i053(:) = 3.05e-7*((300./t_elect(:))**0.5)
     4410         nb_reaction_4 = nb_reaction_4 + 1
     4411         v_4(:,nb_reaction_4) = i053(:)
     4412
     4413!---     i054: H3O+ + e -> H2O + H
     4414
     4415         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
     4416         
     4417         i054(:) = 7.09e-8*((300./t_elect(:))**0.5)
     4418         nb_reaction_4 = nb_reaction_4 + 1
     4419         v_4(:,nb_reaction_4) = i054(:)
     4420
     4421!---     i055: H3O+ + e -> OH + H2
     4422
     4423         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
     4424
     4425         i055(:) = 5.37e-8*((300./t_elect(:))**0.5)
     4426         nb_reaction_4 = nb_reaction_4 + 1
     4427         v_4(:,nb_reaction_4) = i055(:)
     4428
     4429!---     i056: H3O+ + e -> O + H2 + H
     4430
     4431         !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874
     4432
     4433         i056(:) = 5.6e-9*((300./t_elect(:))**0.5)
     4434         nb_reaction_4 = nb_reaction_4 + 1
     4435         v_4(:,nb_reaction_4) = i056(:)
     4436
     4437         nb_reaction_4 = nb_reaction_4 + 1
     4438         v_4(:,nb_reaction_4) = i056(:)
     4439
     4440!---     i057: O+ + H2 -> OH+ + H
     4441
     4442         !UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978)
     4443
     4444         i057(:) = 1.7e-9
     4445         nb_reaction_4 = nb_reaction_4 + 1
     4446         v_4(:,nb_reaction_4) = i057(:)
     4447
     4448!---     i058: OH+ + O -> O2+ + H
     4449
     4450         !UMIST, from Prasad & Huntress, 1980, ApJS, 43, 1
     4451
     4452         i058(:) = 7.1e-10
     4453         nb_reaction_4 = nb_reaction_4 + 1
     4454         v_4(:,nb_reaction_4) = i058(:)
     4455
     4456!---     i059: OH+ + CO2 -> HCO2+ + O
     4457
     4458         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     4459
     4460         i059(:) = 1.44e-9
     4461         nb_reaction_4 = nb_reaction_4 + 1
     4462         v_4(:,nb_reaction_4) = i059(:)
     4463
     4464!---     i060: OH+ + CO -> HCO+ + O
     4465
     4466         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     4467
     4468         i060(:) = 1.05e-9
     4469         nb_reaction_4 = nb_reaction_4 + 1
     4470         v_4(:,nb_reaction_4) = i060(:)
     4471
     4472!---     i061: OH+ + NO -> NO+ + OH (tasa de reacción UMIST 3.59e-10)
     4473
     4474         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     4475
     4476         i061(:) = 3.59e-10
     4477         nb_reaction_4 = nb_reaction_4 + 1
     4478         v_4(:,nb_reaction_4) = i061(:)
     4479
     4480!---     i062: OH+ + H2 -> H2O+ + H (tasa de reacción UMIST 1.01e-9,
     4481
     4482         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     4483
     4484         i062(:) = 1.01e-9
     4485         nb_reaction_4 = nb_reaction_4 + 1
     4486         v_4(:,nb_reaction_4) = i062(:)
     4487
     4488!---     i063: OH+ + O2 -> O2+ + OH (tasa de reacción UMIST 5.9e-10
     4489
     4490         !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981)
     4491
     4492         i063(:) = 5.9e-10
     4493         nb_reaction_4 = nb_reaction_4 + 1
     4494         v_4(:,nb_reaction_4) = i063(:)
     4495
     4496      end if   ! ok_ionchem
     4497
     4498!----------------------------------------------------------------------
    31284499!     reactions avec 02(Dg)
    31294500!----------------------------------------------------------------------
    31304501
    3131 !---     i001: O2(Dg) + CO2 -> O2 + CO2 + hv
    3132 
    3133 !        Krasnopolsky (2010a)
    3134 
    3135       i001(:) = 1.E-20
     4502!---     j001: O2(Dg) + (CO2 and O) -> O2 + (CO2 and O) + hv
     4503
     4504!        Krasnopolsky (2010a) for CO2 & Clark and Wayne, 1969 for O (JPL)
     4505
     4506      j001(:) = 1.E-20 * c(:,i_co2) !+ 2.E-16 * c(:,i_o)
    31364507
    31374508      nb_phot = nb_phot + 1
    3138       v_phot(:,nb_phot) = i001(:)*c(:,i_co2)
    3139 
    3140 !---     i002: O2(Dg) -> O2 + hv
     4509      v_phot(:,nb_phot) = j001(:)
     4510
     4511!---     j002: O2(Dg) -> O2 + hv
    31414512
    31424513!        Lafferty et al; (1998)
    31434514
    3144       i002(:) = 2.2E-4
     4515      j002(:) = 2.2E-4
    31454516
    31464517      nb_phot = nb_phot + 1
    3147       v_phot(:,nb_phot) = i002(:)
     4518      v_phot(:,nb_phot) = j002(:)
     4519
     4520      !!! TEST: artificial increase of CO2 photodissociation
     4521      if (tuneupperatm) then
     4522      !-- TuneA
     4523      !   v_phot(65:78,4) = v_phot(65:78,4)*10.
     4524      !   v_phot(60:64,4) = v_phot(60:64,4)*3.
     4525      !   v_phot(55:59,3) = v_phot(55:59,3)*2.
     4526      !--
     4527      !-- TuneB
     4528      !   v_phot(65:78,4) = v_phot(65:78,4)*10.
     4529      !   v_phot(55:59,3) = v_phot(55:59,3)*5.
     4530      !--
     4531      !-- TuneC
     4532      ! VCD 1.1 tuning
     4533      !   v_phot(65:78,4) = v_phot(65:78,4)*10.
     4534      !   v_phot(52:59,3) = v_phot(52:59,3)*5.
     4535      !--
     4536      !-- TuneE
     4537      ! VCD 2.0 tuning
     4538      !  v_phot(65:90,4) = v_phot(65:90,4)*10. ! CO2 + hv ==> O(1D) + CO
     4539      !--
     4540      ! VCD 2.1 tuning
     4541          v_phot(65:90,4) = v_phot(65:90,4)*8. ! CO2 + hv ==> O(1D) + CO
     4542      !--
     4543      !   v_phot(:,4) = v_phot(:,4)*10.
     4544      !do ij=3,4
     4545      !   v_phot(:,ij) = v_phot(:,ij)*10.
     4546      !end do
     4547      endif
     4548      !!!!!!!!!!!!!!!
    31484549
    31494550return
    … …  
    31764577                               ! number of processes treated numerically as photodissociations
    31774578
    3178 real, dimension(nlayer,nesp)              :: c    ! number densities
    3179 real, dimension(nlayer,      nb_phot_max) :: v_phot
    3180 real, dimension(nlayer,nb_reaction_3_max) :: v_3
    3181 real, dimension(nlayer,nb_reaction_4_max) :: v_4
     4579real (kind = 8), dimension(nlayer,nesp)              :: c    ! number densities
     4580real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
     4581real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
     4582real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
    31824583
    31834584! output
    31844585
    3185 real, dimension(nesp,nesp), intent(out) :: mat  ! matrix
    3186 real, dimension(nesp), intent(out)      :: prod, loss, lossconc
     4586real (kind = 8), dimension(nesp,nesp), intent(out) :: mat  ! matrix
     4587real (kind = 8), dimension(nesp), intent(out)      :: prod, loss, lossconc
    31874588
    31884589! local
    … …  
    31944595integer :: iphot,i3,i4
    31954596
    3196 real :: eps, eps_4  ! implicit/explicit coefficient
     4597real(kind = 8) :: eps, eps_4  ! implicit/explicit coefficient
    31974598
    31984599! initialisations
    … …  
    32994700
    33004701real :: dtold, ctimestep
    3301 real, dimension(nesp)      :: cold, ccur
    3302 real, dimension(nesp,nesp) :: mat1
    3303 real, dimension(nesp)      :: prod, loss
     4702real (kind = 8), dimension(nesp)      :: cold, ccur
     4703real (kind = 8), dimension(nesp,nesp) :: mat1
     4704real (kind = 8), dimension(nesp)      :: prod, loss
    33044705real                       :: dens
    33054706real                       :: lon, lat
    … …  
    33114712! local
    33124713
    3313 real, dimension(nesp)      :: cnew
    3314 real, dimension(nesp,nesp) :: mat
    3315 real :: atol, ratio, e, es, coef
     4714real (kind = 8), dimension(nesp)      :: cnew
     4715real (kind = 8), dimension(nesp,nesp) :: mat
     4716real (kind = 8) :: atol, ratio, e, es, coef
    33164717
    33174718integer                  :: code, iesp, iter
    … …  
    33234724! parameters
    33244725
    3325 real, parameter    :: dtmin   = 10.      ! minimum time step (s)
    3326 real, parameter    :: vmrtol  = 1.e-11   ! absolute tolerance on vmr
    3327 real, parameter    :: rtol    = 0.05     ! rtol recommended value : 0.1-0.02
    3328 integer, parameter :: niter   = 3        ! number of iterations
    3329 real, parameter    :: coefmax = 2.
    3330 real, parameter    :: coefmin = 0.1
    3331 logical            :: fast_guess = .true.
     4726real (kind = 8), parameter :: dtmin   = 10.      ! minimum time step (s)
     4727real (kind = 8), parameter :: vmrtol  = 1.e-11   ! absolute tolerance on vmr
     4728real (kind = 8), parameter :: rtol    = 0.05     ! rtol recommended value : 0.1-0.02
     4729integer,         parameter :: niter   = 3        ! number of iterations
     4730real (kind = 8), parameter :: coefmax = 2.
     4731real (kind = 8), parameter :: coefmin = 0.1
     4732logical                    :: fast_guess = .true.
     4733
    33324734
    33334735dttest = dtold   ! dttest = dtold = dt_guess
    … …  
    34034805
    34044806!======================================================================
     4807!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     4808!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     4809! CE CODE EST OBSOLETE !! A NE PAS UTILISER !!!!!!!!!!!!!!!!!!!!!!!!!!!
     4810!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     4811!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    34054812
    34064813      SUBROUTINE  rate_save(            &

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

    r2799 r2836  
    569569      implicit none
    570570
     571#include "clesphys.h" ! fixed_euv_value
     572
    571573!     input
    572574 
    … …  
    582584      integer, parameter :: kdata = 20000    ! max dimension of input solar flux
    583585      integer :: msun   ! choice of solar flux
    584       integer :: iw, nhead, ihead, n, i, ierr, kin
     586      integer :: iw, nhead, ihead, n, i, ierr, kin1, kin2
    585587
    586588      real, parameter :: deltax = 1.e-4
    587       real, dimension(kdata) :: x1, y1      ! input solar flux
    588       real, dimension(nw)    :: yg1         ! gridded solar flux
    589 
    590       character(len=100) :: fil
    591 
    592       kin = 10    ! input logical unit
    593 
     589      !  atlas1_thuillier_tuv.txt
     590      real, dimension(kdata) :: x1, y1      ! input solar flux - HIGH SOLAR ACTIVITY
     591      real :: E107_max = 196
     592      !  atlas3_thuillier_tuv.txt
     593      real, dimension(kdata) :: x2, y2      ! input solar flux - "LOW" SOLAR ACTIVITY
     594      real :: E107_min = 97
     595      ! be careful, x2 need to be equal to x1
     596      real, dimension(kdata) :: x0, y0      ! input solar flux - interpolated solar activity
     597      real, dimension(nw)    :: yg0         ! gridded solar flux
     598      real            :: factor_interp
     599           
     600      character(len=100) :: fil1, fil2
     601
     602      kin1 = 10    ! input logical unit
     603      kin2 = 11    ! input logical unit
     604     
    594605! select desired extra-terrestrial solar irradiance, using msun:
    595606
    596 ! 18 = atlas1_thuillier_tuv.txt  0-900 nm  March 1992
     607! 18 = atlas1_thuillier_tuv.txt  0-900 nm  March 29, 1992
     608!      Article: F10.7 =  192 s.f.u | <F10.7>81d =  171 s.f.u | <Rz> = 121 (sunsport number)
    597609!      Thuillier et al., Adv. Space. Res., 34, 256-261, 2004
    598  
    599       msun = 18
    600 
    601       if (msun == 18) THEN
    602 
    603          fil = 'solar_fluxes/atlas1_thuillier_tuv.txt'
    604          print*, 'solar flux : ', fil
    605 
    606       if(is_master) then
    607 
    608          open(kin, file=fil, status='old', iostat=ierr)
    609 
    610          if (ierr /= 0) THEN
    611             write(*,*)'cant find solar flux : ', fil
    612             write(*,*)'It should be in : INPUT/solar_fluxes'
    613             write(*,*)'1) You can change this directory address in '
    614             write(*,*)'   callphys.def with datadir=/path/to/dir'
    615             write(*,*)'2) If necessary, /solar fluxes (and other datafiles)'
    616             write(*,*)'   can be obtained online on:'
    617             write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
    618             stop
    619          end if
    620 
    621          nhead = 9
    622          n = 19193
    623          DO ihead = 1, nhead
    624             READ(kin,*)
    625          ENDDO
    626          DO i = 1, n
    627             READ(kin,*) x1(i), y1(i)
    628             y1(i) = y1(i)*1.e-3    ! mw -> w
    629          ENDDO
    630          CLOSE (kin)
    631          CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
    632          CALL addpnt(x1,y1,kdata,n,          0.,0.)
    633          CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
    634          CALL addpnt(x1,y1,kdata,n,      1.e+38,0.)
    635          CALL inter2(nw,wl,yg1,n,x1,y1,ierr)
    636 
    637          IF (ierr .NE. 0) THEN
    638             WRITE(*,*) ierr, fil
    639             STOP
    640          ENDIF
     610!      Model SOLAR2000 version 2015/12: E10.7 = 195.8 s.f.u | <E10.7>81d = 196.9 s.f.u
     611!      Choix de la limite haute: E10.7 = 196 s.f.u
     612
     613! 20 = atlas3_thuillier_tuv.txt  0-900 nm  November 11, 1994
     614!      Article: F10.7 = 77.5 s.f.u | <F10.7>81d = 83.5 s.f.u | <Rz> = 20 (sunsport number)   
     615!      Thuillier et al., Adv. Space. Res., 34, 256-261, 2004
     616!      Model SOLAR2000 version 2015/12: E10.7 =  96.9 s.f.u | <E10.7>81d = 100.0 s.f.u
     617!      Choix de la limite basse: E10.7 =  97 s.f.u
     618
     619      if (fixed_euv_value .ge. 196) THEN
     620        msun = 18
     621        print*, 'Atlas1 spectrum Thuiller chosen'
     622      else
     623        msun = 19
     624        print*, 'interpolated Spectrum case'
     625      end if
     626
     627      IF (is_master) THEN
     628     
     629        fil1 = 'solar_fluxes/atlas1_thuillier_tuv.txt'
     630        print*, 'solar flux : ', fil1
     631       
     632        open(kin1, file=fil1, status='old', iostat=ierr)
     633       
     634        if (ierr /= 0) THEN
     635          write(*,*)'cant find solar flux : ', fil1
     636          write(*,*)'It should be in : INPUT/solar_fluxes'
     637          write(*,*)'1) You can change this directory address in '
     638          write(*,*)'   callphys.def with datadir=/path/to/dir'
     639          write(*,*)'2) If necessary, /solar fluxes (and other datafiles)'
     640          write(*,*)'   can be obtained online on:'
     641          write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
     642          stop
     643        end if
     644         
     645        nhead = 9
     646        n = 19193
     647        DO ihead = 1, nhead
     648          READ(kin1,*)
     649        ENDDO
     650        DO i = 1, n
     651          READ(kin1,*) x1(i), y1(i)
     652          y1(i) = y1(i)*1.e-3    ! mw -> w
     653        ENDDO
     654        CLOSE (kin1)
     655         
     656        fil2 = 'solar_fluxes/atlas3_thuillier_tuv.txt'
     657        print*, 'solar flux : ', fil2
     658       
     659        open(kin2, file=fil2, status='old', iostat=ierr)
     660       
     661        if (ierr /= 0) THEN
     662          write(*,*)'cant find solar flux : ', fil2
     663          write(*,*)'It should be in : INPUT/solar_fluxes'
     664          write(*,*)'1) You can change this directory address in '
     665          write(*,*)'   callphys.def with datadir=/path/to/dir'
     666        write(*,*)'2) If necessary, /solar fluxes (and other datafiles)'
     667          write(*,*)'   can be obtained online on:'
     668        write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
     669          stop
     670        end if
     671
     672        nhead = 9
     673        n = 19193
     674        DO ihead = 1, nhead
     675          READ(kin2,*)
     676        ENDDO
     677        DO i = 1, n
     678          READ(kin2,*) x2(i), y2(i)
     679          y2(i) = y2(i)*1.e-3    ! mw -> w
     680        ENDDO
     681        CLOSE (kin2)
     682
     683        IF (msun .eq. 18) THEN
     684          DO i = 1, n
     685            x0(i) = x1(i)
     686            y0(i) = y1(i)
     687          ENDDO
     688        ELSE
     689          ! INTERPOLATION BETWEEN E107_min to E107_max and extrapolation below E107_min
     690          factor_interp=(fixed_euv_value-E107_min)/(E107_max-E107_min)
     691          DO i = 1, n
     692            x0(i) = x1(i)
     693            y0(i) = y2(i) + (y1(i) - y2(i))* factor_interp
     694          ENDDO
     695        ENDIF ! msun .ne. 18
     696       
     697        CALL addpnt(x0,y0,kdata,n,x0(1)*(1.-deltax),0.)
     698        CALL addpnt(x0,y0,kdata,n,          0.,0.)
     699        CALL addpnt(x0,y0,kdata,n,x0(n)*(1.+deltax),0.)
     700        CALL addpnt(x0,y0,kdata,n,      1.e+38,0.)
     701        CALL inter2(nw,wl,yg0,n,x0,y0,ierr)
     702
     703        IF (ierr .NE. 0) THEN
     704          WRITE(*,*) ierr, fil1, fil2
     705          STOP
     706        ENDIF
    641707
    642708!     convert to photon.s-1.nm-1.cm-2
    643709!     5.039e11 = 1.e-4*1e-9/(hc = 6.62e-34*2.998e8)
    644710
    645          DO iw = 1, nw-1
    646             f(iw) = yg1(iw)*wc(iw)*5.039e11
    647          ENDDO
     711        DO iw = 1, nw-1
     712          f(iw) = yg0(iw)*wc(iw)*5.039e11
     713        ENDDO
    648714
    649715      endif !is_master
    650716
    651717      call bcast(f)
    652 
    653       end if
    654718
    655719      end subroutine rdsolarflux

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

    r2795 r2836  
    1313
    1414      implicit none
    15 
     15     
    1616!     input
    1717
    … …  
    433433   
    434434!           production of o(1d) for wavelengths shorter than 167 nm
    435 
    436             yieldco2(:) = 1
    437435
    438436            if (wc(l) >= 167.) then

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

    r2795 r2836  
    8080      use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
    8181      use sed_and_prod_mad, only: aer_sedimentation, drop_sedimentation
     82      use iono_h, only: temp_elect, temp_ion
    8283#ifdef CPP_XIOS     
    8384      use xios_output_mod, only: initialize_xios_output,
    … …  
    237238c Declaration des procedures appelees
    238239c
    239       EXTERNAL ajsec     ! ajustement sec
    240       EXTERNAL clmain    ! couche limite
    241       EXTERNAL clmain_ideal    ! couche limite simple
    242       EXTERNAL hgardfou  ! verifier les temperatures
    243 c     EXTERNAL orbite    ! calculer l'orbite
    244       EXTERNAL phyetat0  ! lire l'etat initial de la physique
    245       EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
    246       EXTERNAL radlwsw   ! rayonnements solaire et infrarouge
    247 !      EXTERNAL suphec    ! initialiser certaines constantes
    248       EXTERNAL transp    ! transport total de l'eau et de l'energie
     240      EXTERNAL ajsec         ! ajustement sec
     241      EXTERNAL clmain        ! couche limite
     242      EXTERNAL clmain_ideal  ! couche limite simple
     243      EXTERNAL hgardfou      ! verifier les temperatures
     244c     EXTERNAL orbite         ! calculer l'orbite
     245      EXTERNAL phyetat0      ! lire l'etat initial de la physique
     246      EXTERNAL phyredem      ! ecrire l'etat de redemarrage de la physique
     247      EXTERNAL radlwsw       ! rayonnements solaire et infrarouge
     248!      EXTERNAL suphec        ! initialiser certaines constantes
     249      EXTERNAL transp        ! transport total de l'eau et de l'energie
    249250      EXTERNAL printflag
    250251      EXTERNAL zenang
    … …  
    270271c
    271272CXXX PB
    272       REAL fluxt(klon,klev)   ! flux turbulent de chaleur
    273       REAL fluxu(klon,klev)   ! flux turbulent de vitesse u
    274       REAL fluxv(klon,klev)   ! flux turbulent de vitesse v
     273      REAL fluxt(klon,klev)     ! flux turbulent de chaleur
     274      REAL fluxu(klon,klev)     ! flux turbulent de vitesse u
     275      REAL fluxv(klon,klev)     ! flux turbulent de vitesse v
    275276c
    276277      REAL flux_dyn(klon,klev)  ! flux de chaleur produit par la dynamique
    … …  
    278279      REAL flux_ec(klon,klev)   ! flux de chaleur Ec
    279280c
    280       REAL tmpout(klon,klev)  ! K s-1
     281      REAL tmpout(klon,klev)    ! [K/s]
    281282c
    282283      REAL dist, rmu0(klon), fract(klon)
    … …  
    334335c Tendencies due to molecular viscosity and conduction
    335336      real d_t_conduc(klon,klev)     ! [K/s]
    336       real d_u_molvis(klon,klev)     ! (m/s) /s
    337       real d_v_molvis(klon,klev)     ! (m/s) /s
     337      real d_u_molvis(klon,klev)     ! [m/s] /s
     338      real d_v_molvis(klon,klev)     ! [m/s] /s
    338339
    339340c Tendencies due to molecular diffusion
    … …  
    820821      IF (ancien_ok) THEN
    821822         DO k = 1, klev
    822          DO i = 1, klon
     823          DO i = 1, klon
    823824            d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
    824825            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
    825          ENDDO
     826          ENDDO
    826827         ENDDO
    827828
    … …  
    860861c
    861862      DO k = 1, klev
    862       DO i = 1, klon
     863       DO i = 1, klon
    863864         zphi(i,k) = pphi(i,k) + pphis(i)
    864       ENDDO
     865       ENDDO
    865866      ENDDO
    866867
    … …  
    996997               DO ig=1,klon
    997998                   t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
    998                    u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
    999                    v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
     999                   u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! [m/s]
     1000                   v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! [m/s]
    10001001               ENDDO
    10011002            ENDDO
    … …  
    16091610              DO ig=1,klon
    16101611                 d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k)
    1611                
    16121612              ENDDO
    16131613           ENDDO
    … …  
    20722072      CALL send_xios_field("dugwo",d_u_oro)
    20732073      CALL send_xios_field("dugwno",d_u_hin)
     2074      CALL send_xios_field("dvgwno",-1.*d_v_hin)
    20742075      CALL send_xios_field("dumolvis",d_u_molvis)
    20752076c VENUS: regardee a l envers!!!!!!!!!!!!!!!
    … …  
    21022103! production and destruction rate, cm-3.s-1
    21032104! Beware of the context*.xml file !!
    2104 !         if ((tr_scheme == 3) .and. (ok_chem)) THEN
    2105 !            do iq = 1,nqmax - nmicro
    2106 !               !if ((iq == i_co) .or. (iq == i_o)) then
    2107 !               !   call send_xios_field('prod_'+tname(iq), prod_tr(:,:,iq))
    2108 !               !   call send_xios_field('loss_'+tname(iq), loss_tr(:,:,iq))
    2109 !               !end if
    2110 !               if (iq == i_o) then
    2111 !                  call send_xios_field('prod_o', prod_tr(:,:,iq))
    2112 !                  call send_xios_field('loss_o', loss_tr(:,:,iq))
    2113 !               end if
    2114 !               if (iq == i_co) then
    2115 !                  call send_xios_field('prod_co', prod_tr(:,:,iq))
    2116 !                  call send_xios_field('loss_co', loss_tr(:,:,iq))
    2117 !               end if
    2118 !            end do
    2119 !         end if
     2105         if ((tr_scheme == 3) .and. (ok_chem)) THEN
     2106            do iq = 1,nqmax - nmicro
     2107               if ((iq.eq.i_o).or.(iq.eq.i_co)) THEN
     2108                 call send_xios_field("prod_"//tname(iq),
     2109     $                                         prod_tr(:,:,iq))
     2110                 call send_xios_field("loss_"//tname(iq),
     2111     $                                         loss_tr(:,:,iq))
     2112               end if
     2113               !if (iq.eq.i_o) then
     2114               !   call send_xios_field('prod_o', prod_tr(:,:,iq))
     2115               !   call send_xios_field('loss_o', loss_tr(:,:,iq))
     2116               !end if
     2117               !if (iq.eq.i_co) then
     2118               !   call send_xios_field('prod_co', prod_tr(:,:,iq))
     2119               !   call send_xios_field('loss_co', loss_tr(:,:,iq))
     2120               !end if
     2121            end do
     2122         end if
    21202123
    21212124! tracers in gas phase, volume mixing ratio
    … …  
    21302133         do iq = 1,nqmax - nmicro
    21312134            col_dens_tr(:,iq)=0.
    2132             do k = 1, klev
    2133               col_dens_tr(:,iq) = col_dens_tr(:,iq) +
    2134      $              tr_seri(:,k,iq)* (paprs(:,k)-paprs(:,k+1)) / RG
    2135             end do
    2136             call send_xios_field("col_"//tname(iq),col_dens_tr(:,iq))
     2135            if (type_tr(iq).eq.1) THEN
     2136               do k = 1, klev
     2137                 col_dens_tr(:,iq) = col_dens_tr(:,iq) +
     2138     $               tr_seri(:,k,iq)* (paprs(:,k)-paprs(:,k+1)) / RG
     2139               end do
     2140               call send_xios_field("col_"//tname(iq),col_dens_tr(:,iq))
     2141            end if
    21372142         end do
    21382143
    … …  
    21512156! aeronomical emissions
    21522157
    2153 !        call send_xios_field("no_emission",no_emission)
    2154 !        call send_xios_field("o2_emission",o2_emission)
     2158        call send_xios_field("no_emis",no_emission)
     2159        call send_xios_field("o2_emis",o2_emission)
    21552160
    21562161! chemical iterations

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

    r2795 r2836  
    2121      use chemparam_mod
    2222      use conc, only: mmean,rnew
     23      use photolysis_mod, only: init_photolysis, nphot
     24      use iono_h, only: temp_elect
    2325#ifdef CPP_XIOS     
    2426      use xios_output_mod, only: send_xios_field
    … …  
    6365      real :: lon_sun
    6466      real :: zlocal(nlev)  ! altitude for photochem (km)
    65 
     67      real :: t_elec(nlev)  ! electron temperature [K]
     68     
     69      integer, parameter :: t_elec_origin=2
     70      !Electronic temperature. Index 1 -> Theis et al. 1980 - model data ; Index 2-> Theis et al. 1984 - model data
     71     
    6672      integer :: i, iq
    6773      integer :: ilon, ilev
    … …  
    7682      real, dimension(nlev) :: no_em
    7783      real, dimension(nlev) :: o2_em
    78      
     84
     85      integer, save :: nb_reaction_3_max     ! number of quadratic reactions
     86      integer, save :: nb_reaction_4_max     ! number of bimolecular reactions
     87      integer, save :: nquench               ! number of quenching + heterogeneous reactions
     88      integer, save :: nphotion              ! number of photoionizations
     89      integer, save :: nb_reaction_4_ion     ! quadratic reactions for ionosphere
     90!      integer, save :: nb_reaction_4_deut    ! quadratic reactions for deuterium chem
     91      integer, save :: nb_phot_max           ! total number of photolysis+photoionizations+quenching reactions
     92
     93      ! tracer indexes for the EUV heating:
     94!!! ATTENTION. These values have to be identical to those in euvheat.F90
     95!!! If the values are changed there, the same has to be done here  !!!
     96
     97!      integer,parameter :: i_co2=1
     98!      integer,parameter :: i_n2=13
     99!      integer,parameter :: i_n=14
     100!      integer,parameter :: i_o=3
     101!      integer,parameter :: i_co=4
     102
     103      integer,parameter :: ix_co2  =  1
     104      integer,parameter :: ix_co   =  2
     105      integer,parameter :: ix_o    =  3
     106      integer,parameter :: ix_o1d  =  4
     107      integer,parameter :: ix_o2   =  5
     108      integer,parameter :: ix_o3   =  6
     109      integer,parameter :: ix_h    =  7
     110      integer,parameter :: ix_h2   =  8
     111      integer,parameter :: ix_oh   =  9
     112      integer,parameter :: ix_ho2  = 10
     113      integer,parameter :: ix_h2o2 = 11
     114      integer,parameter :: ix_h2o  = 12
     115      integer,parameter :: ix_n    = 13
     116      integer,parameter :: ix_n2d  = 14
     117      integer,parameter :: ix_no   = 15
     118      integer,parameter :: ix_no2  = 16
     119      integer,parameter :: ix_n2   = 17
     120
     121      ! NEED TO BE THE SAME THAN IN EUVHEAT.F90
     122      integer,parameter :: nespeuv=17    ! Number of species considered (11, 12 or 17 (with nitrogen))
     123
     124      real :: vmr_dens_euv(nlon,nlev,nespeuv) ! local species density for EUV heating
     125         
    79126!===================================================================
    80127!     first call : initialisations
    … …  
    90137
    91138!-------------------------------------------------------------------
     139!        Determination of chemistry reaction number
     140!-------------------------------------------------------------------
     141
     142         if (ok_chem) then
     143           ! set number of reactions, depending on ion chemistry or not
     144           nb_reaction_4_ion  = 64
     145           !nb_reaction_4_deut = 35
     146   
     147           !Default numbers if no ion and no deuterium chemistry included
     148   
     149           nb_reaction_4_max = 98     ! set number of bimolecular reactions
     150           nb_reaction_3_max = 12     ! set number of quadratic reactions
     151           nquench           = 13     ! set number of quenching + heterogeneous
     152           nphotion          = 0      ! set number of photoionizations
     153   
     154           if (ok_ionchem) then
     155              nb_reaction_4_max = nb_reaction_4_max+nb_reaction_4_ion 
     156              nphotion          = 18   ! set number of photoionizations
     157           endif
     158           !if(deutchem) then
     159           !   nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_deut 
     160           !end if
     161   
     162           !nb_phot_max is the total number of processes that are treated
     163           !numerically as a photolysis:
     164   
     165           nb_phot_max = nphot + nphotion + nquench
     166   
     167!-------------------------------------------------------------------
    92168!        case of tracers re-initialisation with chemistry
    93169!-------------------------------------------------------------------
    94 
    95          if (reinit_trac .and. ok_chem) then
    96 
    97 !!! in this reinitialisation, trac is VOLUME mixing ratio
    98 ! ONLY SO2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    99 !           convert mass to volume mixing ratio
    100             do iq = 1,nqmax - nmicro
    101                trac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq)
    102             end do
    103             print*, "SO2 is re-initialised"
    104             if (i_so2 /= 0) then
    105                trac(:,1:22,i_so2) = 10.e-6
    106 
    107 ! AL TRACERS!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    108 !           print*, "Tracers are re-initialised"
    109 !           trac(:,:,:) = 1.0e-30
    110 
    111 !           if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0)
    112 !    $           .and. (i_so2 /= 0) .and. (i_h2o /= 0) .and. (i_n2/ = 0)
    113 !    $           .and. (i_co2 /= 0)) then
    114 
    115 !              trac(:,1:22,i_ocs) = 3.e-6
    116 !              trac(:,1:22,i_co)  = 25.e-6
    117 !              trac(:,:,i_hcl)    = 0.4e-6
    118 !              trac(:,1:22,i_so2) = 7.e-6
    119 !              trac(:,1:22,i_h2o) = 30.e-6
    120 !              trac(:,:,i_n2)     = 0.35e-1
    121 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    122    
    123 !          ensure that sum of mixing ratios = 1
    124 
    125                trac_sum(:,:) = 0.
    126 
    127                do iq = 1,nqmax - nmicro
    128                   if (iq /= i_co2) then
    129                      trac_sum(:,:) = trac_sum(:,:) + trac(:,:,iq)
    130                   end if
    131                end do
    132 
    133 !          initialise co2
    134 
    135                trac(:,:,i_co2) = 1. - trac_sum(:,:)
    136 
    137             else
    138                write(*,*) "at least one tracer is missing : stop"
    139                stop
    140             end if
    141        
    142 !           update mmean
    143 
    144             mmean(:,:) = 0.
    145             do iq = 1,nqmax - nmicro
    146                mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
    147             enddo
    148             rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
    149 
    150 !           convert volume to mass mixing ratio
    151 
    152             do iq = 1,nqmax - nmicro
    153                trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:)
    154             end do
    155    
    156          end if  ! reinit_trac
     170 
     171           if (reinit_trac .and. ok_chem) then
     172 
     173  !!! in this reinitialisation, trac is VOLUME mixing ratio
     174  ! ONLY SO2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     175  !           convert mass to volume mixing ratio
     176              do iq = 1,nqmax - nmicro
     177                 trac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq)
     178              end do
     179              print*, "SO2 is re-initialised"
     180              if (i_so2 /= 0) then
     181                 trac(:,1:22,i_so2) = 10.e-6
     182 
     183  ! AL TRACERS!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     184  !           print*, "Tracers are re-initialised"
     185  !           trac(:,:,:) = 1.0e-30
     186 
     187  !           if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0)
     188  !   $          .and. (i_so2 /= 0) .and. (i_h2o /= 0)
     189  !   $          .and. (i_n2 /= 0)  .and. (i_co2 /= 0)) then
     190 
     191  !              trac(:,1:22,i_ocs) = 3.e-6
     192  !              trac(:,1:22,i_co)  = 25.e-6
     193  !              trac(:,:,i_hcl)    = 0.4e-6
     194  !              trac(:,1:22,i_so2) = 7.e-6
     195  !              trac(:,1:22,i_h2o) = 30.e-6
     196  !              trac(:,:,i_n2)     = 0.35e-1
     197  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     198     
     199  !          ensure that sum of mixing ratios = 1
     200 
     201                 trac_sum(:,:) = 0.
     202 
     203                 do iq = 1,nqmax - nmicro
     204                    if (iq /= i_co2) then
     205                       trac_sum(:,:) = trac_sum(:,:) + trac(:,:,iq)
     206                    end if
     207                 end do
     208 
     209  !          initialise co2
     210 
     211                 trac(:,:,i_co2) = 1. - trac_sum(:,:)
     212 
     213              else
     214                 write(*,*) "at least one tracer is missing : stop"
     215                 stop
     216              end if
     217         
     218  !           update mmean
     219 
     220              mmean(:,:) = 0.
     221              do iq = 1,nqmax - nmicro
     222                 mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
     223              enddo
     224              rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
     225 
     226  !           convert volume to mass mixing ratio
     227 
     228              do iq = 1,nqmax - nmicro
     229                 trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:)
     230              end do
     231   
     232           end if  ! reinit_trac
     233
     234         end if  ! ok_chem
    157235
    158236!-------------------------------------------------------------------
    … …  
    286364         lon_local(:) = lon(:)*rpi/180.
    287365         lat_local(:) = lat(:)*rpi/180.
    288 
     366         
     367         if (ok_ionchem) then
     368       
     369           vmr_dens_euv(:,:,ix_co2) = ztrac(:,:,i_co2)  ! CO2
     370           vmr_dens_euv(:,:,ix_co)  = ztrac(:,:,i_co)   ! CO
     371           vmr_dens_euv(:,:,ix_o)   = ztrac(:,:,i_o)    ! O
     372           vmr_dens_euv(:,:,ix_o1d) = ztrac(:,:,i_o1d)  ! O(1D)
     373           vmr_dens_euv(:,:,ix_o2)  = ztrac(:,:,i_o2)   ! O2
     374           vmr_dens_euv(:,:,ix_o3)  = ztrac(:,:,i_o3)   ! O3
     375           vmr_dens_euv(:,:,ix_h)   = ztrac(:,:,i_h)    ! H
     376           vmr_dens_euv(:,:,ix_h2)  = ztrac(:,:,i_h2)   ! H2
     377           vmr_dens_euv(:,:,ix_oh)  = ztrac(:,:,i_oh)   ! OH
     378           vmr_dens_euv(:,:,ix_ho2) = ztrac(:,:,i_ho2)  ! HO2
     379           vmr_dens_euv(:,:,ix_h2o2)= ztrac(:,:,i_h2o2) ! H2O2
     380           vmr_dens_euv(:,:,ix_h2o) = ztrac(:,:,i_h2o)  ! H2O
     381           vmr_dens_euv(:,:,ix_n)   = ztrac(:,:,i_n)    ! N
     382           vmr_dens_euv(:,:,ix_n2d) = ztrac(:,:,i_n2d)  ! N(2D)
     383           vmr_dens_euv(:,:,ix_no)  = ztrac(:,:,i_no)   ! NO
     384           vmr_dens_euv(:,:,ix_no2) = ztrac(:,:,i_no2)  ! NO2
     385           vmr_dens_euv(:,:,ix_n2)  = ztrac(:,:,i_n2)   ! N2
     386           
     387         end if
     388         
    289389         do ilon = 1,nlon
    290390            zlocal(:)=zzlay(ilon,:)/1000.
    291 
     391                       
    292392!           solar zenith angle
    293393
    … …  
    296396     $                 *sin(lon_local(ilon))*sin(lon_sun))*180./rpi
    297397
     398!           electron temperature
     399            do ilev = 1, nlev
     400              t_elec(ilev) = temp_elect(zlocal(ilev),temp(ilon,ilev),
     401     $                                    sza_local,t_elec_origin) 
     402            end do
     403
    298404            call photochemistry_venus(nlev, nlon, zlocal, pdtphys,
    299      $                                pplay(ilon,:)/100.,
    300      $                                temp(ilon,:),
    301      $                                ztrac(ilon,:,:),
    302      $                                mmean(ilon,:),
    303      $                                sza_local,
    304      $                                lon(ilon), lat(ilon),
    305      $                                nqmax, iter(ilon,:),
    306      $                                prod_tr(ilon,:,:),
    307      $                                loss_tr(ilon,:,:),
    308      $                                no_em, o2_em)
     405     $                           ok_jonline,ok_ionchem,tuneupperatm,
     406     $                           nb_reaction_3_max,nb_reaction_4_max,
     407     $                           nb_phot_max,nphotion,ilon,         
     408     $                           pplay(ilon,:)/100.,
     409     $                           temp(ilon,:), t_elec(:),
     410     $                           ztrac(ilon,:,:),vmr_dens_euv(ilon,:,:),
     411     $                           mmean(ilon,:),
     412     $                           sza_local,
     413     $                           lon(ilon), lat(ilon),
     414     $                           nqmax, nespeuv, iter(ilon,:),
     415     $                           prod_tr(ilon,:,:),
     416     $                           loss_tr(ilon,:,:),
     417     $                           no_em, o2_em)
    309418
    310419             no_emi(ilon,:) = no_em(:)
Note: See TracChangeset for help on using the changeset viewer.